Combining Taylor-Hood and Raviart Thomas spaces with hp::FECollection

[Eldar Khattatov]

Hi,

I am trying to implement a coupled Stokes-Darcy problem using Step-46 tutorial as an example. So, essentially I would want to use something like the following

FESystem<dim> fe_stokes(FE_Q<dim>(2), dim,       // Stokes pair
                        FE_Q<dim>(1), 1,
                        FE_Nothing<dim>(dim), 1, // Darcy pair
                        FE_Nothing<dim>(), 1,
                        FE_Nothing<dim>(), 1);   // Lagrange multiplier

FESystem<dim> fe_darcy(FE_Nothing<dim>(), dim,
                       FE_Nothing<dim>(), 1,
                       FE_RaviartThomas<dim>(0), 1,
                       FE_DGQ<dim>(0), 1,
                       FE_DGQ<dim>(0), 1);

hp::FECollection<dim> fe;
fe.push_back(fe_stokes);
fe.push_back(fe_darcy);

//...

dof_handler.distribute_dofs (fe);

But this results in the error message:

--------------------------------------------------------

An error occurred in line <936> of file </home/eldar/TmpDeal/source/fe/fe.cc> in function

    dealii::FiniteElementDomination::Domination dealii::FiniteElement<<anonymous>, <anonymous> >::compare_for_face_domination(const dealii::FiniteElement<<anonymous>, <anonymous> >&) const [with int dim = 2; int spacedim = 2]

The violated condition was: 

    false

Additional information: 

    You are trying to use functionality in deal.II that is currently not implemented. In many cases, this indicates that there simply didn't appear much of a need for it, or that the author of the original code did not have the time to implement a particular case. If you hit this exception, it is therefore worth the time to look into the code to find out whether you may be able to implement the missing functionality. If you do, please consider providing a patch to the deal.II development sources (see the deal.II website on how to contribute).

Stacktrace:

-----------

#0  /home/eldar/Libs/DealII_8.5.1/lib/libdeal_II.g.so.9.0.0-pre: dealii::FiniteElement<2, 2>::compare_for_face_domination(dealii::FiniteElement<2, 2> const&) const

#1  /home/eldar/Libs/DealII_8.5.1/lib/libdeal_II.g.so.9.0.0-pre: dealii::FESystem<2, 2>::compare_for_face_domination(dealii::FiniteElement<2, 2> const&) const

From what I could see the implementation of dealii::FiniteElement<2, 2>::compare_for_face_domination(dealii::FiniteElement<2, 2> const&) const  consists of Assert(false, ...) so for some reason what I am trying to do was not meant to be done this way. 

Could you please suggest what else could be done, or, maybe, there is just another way of approaching coupled problems like this one using Deal.II?

Thank you!

[Wolfgang Bangerth]

> But this results in the error message:
> |
> --------------------------------------------------------
> An error occurred in line <936> of file </home/eldar/TmpDeal/source/fe/fe.cc>
> in function
> dealii::FiniteElementDomination::Domination
> dealii::FiniteElement<<anonymous>, <anonymous>
> >::compare_for_face_domination(const dealii::FiniteElement<<anonymous>,
> <anonymous> >&) const [with int dim = 2; int spacedim = 2]
> The violated condition was:
> false

The fact that you end up in the function of the FiniteElement base class means
that one of the finite element classes on which you are calling this function
does not implement this virtual function itself.

From your backtrace, I do see that you're ending up in
FESystem::compare_for_face_domination(), but I can't see where it is trying to
go from there. If you're in a debugger, can you see which element it's
currently trying to compare?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Eldar Khattatov]

It tries to compare FE_RaviartThomas. This element indeed does not implement compare_for_face_domination() function, so I tried using FE_RaviartThomasNodal instead, as it provides an implementation for this function. Unfortunately, from what I understood from the implementation in RT_Nodal, it works only if another element in the comparison is also of the same type RT_Nodal, otherwise it has same Assert(false, ...) before return statement, so I end up with similar exception message:

An error occurred in line <537> of file <.../fe/fe_raviart_thomas_nodal.cc> in function

   virtual FiniteElementDomination::Domination dealii::FE_RaviartThomasNodal<2>::compare_for_face_domination(const FiniteElement<dim> &) const [dim = 2]

The violated condition was:

   false

[Wolfgang Bangerth]

On 08/18/2017 06:19 AM, Eldar Khattatov wrote:
> It tries to compare FE_RaviartThomas. This element indeed does not implement
> compare_for_face_domination() function,

Correct. Though you could think about implementing it.

> so I tried using FE_RaviartThomasNodal
> instead, as it provides an implementation for this function. Unfortunately,
> from what I understood from the implementation in RT_Nodal, it works only if
> another element in the comparison is also of the same type RT_Nodal, otherwise
> it has same Assert(false, ...) before return statement,

Yes, indeed. But it shouldn't be very difficult to extend it for the
particular case you are considering: that the other element is an FE_Nothing.
The FE_Nothing describes a field that is constant zero, and so the appropriate
interface condition is that the FE_RT_Nodal should also be zero at that
interface, or be left unconstrained, depending on how you want the FE_Nothing
to be interpreted. That should be easy to implement.

Take a look at how the FE_Q element does this, in one of the classes that sit
below class FiniteElement and class FE_Q in the class inheritance diagram. I
think you'll see pretty quickly how you want to implement things for your case.

[Eldar Khattatov]

Thank you for the answer, it was indeed very easy to modify the RT_Nodal according to your comments. If this patch could be of use to the library/community I can commit my changes to the repo.

Best,

Eldar

[Wolfgang Bangerth]

On 08/31/2017 08:55 AM, Eldar Khattatov wrote:
> Thank you for the answer, it was indeed very easy to modify the RT_Nodal
> according to your comments. If this patch could be of use to the
> library/community I can commit my changes to the repo.

Yes, please do open a pull request for this!