Error when make_hanging_node_constraints: FE_TraceQ element

[Anton]

Hi,

I am using static condensation to eliminate the internal degrees of freedom from my FEM problem and therefore end up with a system defined on the interfaces only.  When I try to run "make_hanging_node_constraints" after adaptively refining the mesh I end up with an error, which I do not really expect to be there.  It seems that the problem is occurring with FE_TraceQ element.  I have made a small example illustrating the problem (attached).  The error log is below.

If anyone could help me to get around this issue I would really appreciate it.

Sicnerely,

/Anton

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

An error occurred in line <826> of file <../source/fe/fe.cc> in function

    const FullMatrix<double> &dealii::FiniteElement<2, 2>::constraints(const dealii::internal::SubfaceCase<dim> &) const [dim = 2, spacedim = 2]

The violated condition was:

    (this->dofs_per_face == 0) || (interface_constraints.m() != 0)

The name and call sequence of the exception was:

    ExcConstraintsVoid()

Additional Information:

(none)

Stacktrace:

-----------

#0  2   libdeal_II.g.8.2.1.dylib            0x0000000106659cdc _ZNK6dealii13FiniteElementILi2ELi2EE11constraintsERKNS_8internal11SubfaceCaseILi2EEE + 348: 2   libdeal_II.g.8.2.1.dylib            0x0000000106659cdc _ZNK6dealii13FiniteElementILi2ELi2EE11constraintsERKNS_8internal11SubfaceCaseILi2EEE

#1  3   libdeal_II.g.8.2.1.dylib            0x00000001084e766a _ZN6dealii8DoFTools8internal38make_oldstyle_hanging_node_constraintsINS_10DoFHandlerILi2ELi2EEEEEvRKT_RNS_16ConstraintMatrixENS_8internal8int2typeILi2EEE + 3258: 3   libdeal_II.g.8.2.1.dylib            0x00000001084e766a _ZN6dealii8DoFTools8internal38make_oldstyle_hanging_node_constraintsINS_10DoFHandlerILi2ELi2EEEEEvRKT_RNS_16ConstraintMatrixENS_8internal8int2typeILi2EEE

#2  4   libdeal_II.g.8.2.1.dylib            0x00000001084e352e _ZN6dealii8DoFTools29make_hanging_node_constraintsINS_10DoFHandlerILi2ELi2EEEEEvRKT_RNS_16ConstraintMatrixE + 78: 4   libdeal_II.g.8.2.1.dylib            0x00000001084e352e _ZN6dealii8DoFTools29make_hanging_node_constraintsINS_10DoFHandlerILi2ELi2EEEEEvRKT_RNS_16ConstraintMatrixE

#3  5   demo_err                            0x0000000103a8fbeb main + 603: 5   demo_err                            0x0000000103a8fbeb main

#4  6   libdyld.dylib                       0x00007fff88e055c9 start + 1: 6   libdyld.dylib                       0x00007fff88e055c9 start

#5  7   ???                                 0x0000000000000001 0x0 + 1: 7   ???                                 0x0000000000000001 0x0

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

[Daniel Arndt]

Hey Anton,

the error that you get simply tells you that constraints are not yet implemented for FE_TraceQ.

Best
Daniel

[Martin Kronbichler]

> the error that you get simply tells you that constraints are not yet
> implemented for FE_TraceQ.

And to extend on this answer: Since FE_TraceQ represents an FE_Q as seen
from the faces, you may try to implement this quite easily: In the
constructor of fe_trace.cc, you can fill the matrix
'interface_constraints' with the respective matrix you get from FE_Q of
the same degree. Do you want to try that?

Best,
Martin

[Anton Evgrafov]

Sure I could look at that. /Anton

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/mRCm8c-zUzk/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.



--
--Anton Evgrafov

[Anton Evgrafov]

Modifying the constructor of FE_TraceQ fixes half of the problem. The
implementation of constraints_are_implemented() and constraints()
in FiniteElement are (methinks) wrong for FE_FaceQ/P. These are the
elements with dofs on faces but w/o meaningful dof constraints, are
they not? /Anton

--
--Anton Evgrafov

[Martin Kronbichler]

Anton,

> Modifying the constructor of FE_TraceQ fixes half of the problem. The
> implementation of constraints_are_implemented() and constraints()
> in FiniteElement are (methinks) wrong for FE_FaceQ/P. These are the
> elements with dofs on faces but w/o meaningful dof constraints, are
> they not?

The elements FE_FaceQ/P put all degrees of freedom on faces. When two
faces meet, the field is the same from both sides. For this reason, you
also need to have hanging node constraints for FE_FaceQ when the
refinement level from the two sides of the element differs. Therefore,
these elements have constraints just as usual continuous elements, even
though the solution space is discontinuous over vertices (and edges in
3D).

One might argue that the way we use constraints for FE_FaceQ in step-51
is suboptimal, though. Constraints in deal.II are constructed such that
one always takes the solution space from the coarse side. However, in
the HDG context such as step-51 you typically want to use the solution
space from the finer side because that is the richer solution space. We
have not really thought deeply about that, but I guess one could create
reasonable approximations for that case.

Thinking about the original problem a bit more, namely the constraints
for FE_TraceQ, lead me to the conclusion that you should also overload
get_face_interpolation_matrix, get_subface_interpolation_matrix,
hp_constraints_are_implemented, and compare_for_face_domination. See the
definition of FE_FaceQ. Only the first two methods require
implementation, for which you could again use the respective methods in
FE_Q_Base. (These methods are not called very often, so it is OK to
create a temporary FE_Q object within these methods.)

Best,
Martin

[Wolfgang Bangerth]

On 02/19/2015 04:52 AM, Martin Kronbichler wrote:
> Thinking about the original problem a bit more, namely the constraints
> for FE_TraceQ, lead me to the conclusion that you should also overload
> get_face_interpolation_matrix, get_subface_interpolation_matrix,
> hp_constraints_are_implemented, and compare_for_face_domination. See the
> definition of FE_FaceQ. Only the first two methods require
> implementation, for which you could again use the respective methods in
> FE_Q_Base. (These methods are not called very often, so it is OK to
> create a temporary FE_Q object within these methods.)

Or, maybe something like this (basically no modifications required at all to
the constructor and elsewhere):

class FE_TraceQ {
private:
// store an FE_Q element from which we can obtain face interpolation
// matrices, hanging node constraints, etc, without having to
// duplicate the code
FE_Q fe_q;
};

FE_TraceQ::interface_constraints (...) {
return fe_q.interface_constraints(...);
}

Anton, we would certainly accept such a patch!

Best
W.

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

[Anton Evgrafov]

Wolfgang: This indeed only solves a part of the problem as one then
cannot use both TraceQ and FaceQ in a FE_System, since the the latter
has constraints_are_implemented = false but
hp_constraints_are_implemented = true. Maybe I will indeed go ahead
and implement all the face/subface interpolation etc. /Anton

[Martin Kronbichler]

Dear Anton,

> Wolfgang: This indeed only solves a part of the problem as one then
> cannot use both TraceQ and FaceQ in a FE_System, since the the latter
> has constraints_are_implemented = false but
> hp_constraints_are_implemented = true. Maybe I will indeed go ahead
> and implement all the face/subface interpolation etc.

Indeed, you should go ahead with the face/subface interpolation and the
other two methods suggested below. Thank you for your effort!

As far as interface_constraints are concerned, one could replace the
code for interface_constraints (for FE_Q) from combining the possible
get_subface_interpolation_matrix results with the inverse of the
get_face_interpolation_matrix of the element with itself if I read the
code correctly; this is something we should do at some point, but one
would need to check this more closely.

Best,
Martin

[Anton Evgrafov]

Here is a minimal-code duplication approach to FE_TraceQ
implementation, suggested by Wolfgang (attached). I have added a
private copy of FE_Q and have delegated as many calls to it as
possible (though taking get_subface_interpolation_matrix from FE_FaceQ
is another option).



/Anton

[Martin Kronbichler]

Dear Anton,

Thank you very much for the patch. I will take care of it and add a test
to the test suite. If all goes well, we might merge it early next week.

Best,
Martin

[Martin Kronbichler]

Dear Anton,

I've now applied and merged your patch to the development sources of deal.II, see also https://github.com/dealii/dealii/pull/585 for documentation. Thanks once again for the patch.

Best,

Martin

[Anton Evgrafov]

Dear Martin,

Great, thank you.

Regarding the use of FE_FaceQ in the HDG tutorial with the
reconstruction of the local DOFs from the "coarse" trace rather than
from the "fine" trace. I suppose it would be too much to try to add
an alternative reconstruction, which then would allow one to omit the
hanging node constraints, to the tutorial example? This way one would
need to integrate over subfaces and keep track of "fine" element
neighbors/dofs, which seems to go a little bit against the overall
design of the library?

/Anton

[Martin Kronbichler]

Dear Anton,

> Regarding the use of FE_FaceQ in the HDG tutorial with the
> reconstruction of the local DOFs from the "coarse" trace rather than
> from the "fine" trace. I suppose it would be too much to try to add
> an alternative reconstruction, which then would allow one to omit the
> hanging node constraints, to the tutorial example?

In my opinion, adapting it to the general context would indeed make the
tutorial rather difficult to read because, as you correctly point out,
one would need to integrate over subfaces from the coarse side. This
leads to variable size in the element matrix and requires 'modified'
local_dof_indices that take the ones from the refined side. I have not
tried that myself so I cannot answer how much of an improvement this
would give; I definitely observed sub-optimal behavior with the strategy
to use the constraints from the coarse side as we currently do.

One could certainly do this, but to make the code nice requires some
more abstraction that we currently do not have. It is on my agenda, but
not a very high priority topic... :(

On the other hand, I've been curious to what would happen if instead of
constraining the fine space to the coarse space we do the opposite,
i.e., 'interpolate' the coarse information (matrix, rhs) into the dofs
of the refined side. This does not give full rank in itself, but since
we have the representation on the refined side the method should at
least be stable (singe-face HDG is stable). Hence, this would be
something interesting to try at some point.

Best,
Martin