Face and line orientation dependent dof permutations in FE_Q

[Konrad Simon]

Dear Deal.II community,

I am currently fixing an orientation issue for some (non-primitive) vector elements with the help of two other Deal.II developers/maintainers.

I have a quick question: On complicated meshes some faces can have a non-standard orientation or they can be rotated against the neighboring face. This causes inconsistencies for some vector elements.

In principle this potentially causes inconsistencies for FE_Q elements as well (for all dofs located at vertices, lines and faces). But that is being taken care of in the library, so everything seems fine.

In order to understand the issue better for other elements: where in the library is this being taken care of? I have difficulties to find that. 

Remark: I found a table in (fe.h/fe.cc) that derived elements need to implement that takes care of local dof permutations on quads (not full faces and hence not line dofs or vertex dofs). Now, lines can have different orientations but lines also permute within a cell upon face orientation/rotation changes.

Can anyone point me to the place in Deal.II where this is being taken care of for FE_Q?

Best,
Konrad

[Wolfgang Bangerth]

On 1/13/21 12:32 PM, Konrad Simon wrote:
>
>
> I have a quick question: On complicated meshes some faces can have a
> non-standard orientation or they can be rotated against the neighboring face.
> This causes inconsistencies for some vector elements.
>
> In principle this potentially causes inconsistencies for FE_Q elements as well
> (for all dofs located at vertices, lines and faces). But that is being taken
> care of in the library, so everything seems fine.
>
> In order to understand the issue better for other elements: where in the
> library is this being taken care of? I have difficulties to find that.
>
> Remark: I found a table in (fe.h/fe.cc) that derived elements need to
> implement that takes care of local dof permutations on quads (not full faces
> and hence not line dofs or vertex dofs). Now, lines can have different
> orientations but lines also permute within a cell upon face
> orientation/rotation changes.
>
> Can anyone point me to the place in Deal.II where this is being taken care of
> for FE_Q?

Konrad -- my apologies not only for not getting to this question quicker, but
in fact for not helping you along over the past few months. I saw that you are
working on this, but it's been a very long time since I understood these
issues sufficiently well to help out and I haven't had the time to really read
up on it again in enough detail to be useful. I do know that it's a topic for
which we, collectively, have probably lost our expertise and that there is
likely nobody else who knows it well either. You might be the expert on this
at the moment.

Fundamentally, the place to look is to search through the finite element
classes where they ask for face and line orientations. It is possible that you
don't actually have to do much there: Finite element classes generally work in
the local (cell) coordinate system, where you enumerate shape functions as
0...dofs_per_cell. So when, for example, you compute the local matrix, it
really is local. The interesting part happens when you translate what *global*
degrees of freedom correspond to the local ones, i.e., when you call
cell->get_dof_indices() (==DoFCellAccessor::get_dof_indices). This is where on
every cell has to collect which DoFs live on it, and that has to take into
account the orientations of lines and faces come into play. The cell DoF
indices are cached, but you can see how this works in
dof_accessor.templates.h, starting in line 902 using the process_dof_indices()
functions.

What this means is that if you have two cells that disagree about the relative
orientation of their shared face, they when you call cell->get_dof_indices()
on these two cells, you'd end up with a different ordering of the DoFs on the
face. My usual approach to check this would be to create a mesh with two cells
so that one has a non-standard face orientation, assign a Q3 element to the
DoFHandler, and check that that is really what happens. (I generally turn
these little programs into tests for the test suite, just so that what I just
verified to be correct -- or not! -- is actually checked automatically in the
future. I would very much love to see you submit such tests as well!)

So this is the process for FE_Q, where there is really not very much to do, or
in fact for all interpolatory elements. The situation becomes a bit more
difficult for things such as the Nedelec element (also the Raviart-Thomas one)
where the shape functions are tangential to the edge. Let's say you have the
lowest order Nedelec element. Then there is just one DoF on each edge, and
both adjacent cells agree about this. But they also have to agree *which
direction the vector-valued shape function should point to*, and that's
trickier. I can't remember how we do this (or whether we actually do this at
all). I *think* that the RT element gets this right in 2d, but I don't
remember about the 3d case. Again, it would probably be quite helpful to write
little programs and turn them into tests.

I know or strongly suspect that there are bugs for RT and Nedelec elements on
meshes where faces have non-standard orientation. I've been thinking for years
how I would go about debugging this, but I've never had the time to actually
do it. But I'll share what I think would be the tools to do this:

1/ Creating meshes with non-standard orientations: You can feed
Triangulation::create_triangulation() directly so that the faces have
non-standard orientations with regard to each other. (See step-14 for how to
call this function.) In 2d, this would work as follows: Think of a mesh like this:

3---4---5
| | |
0---1---2

You'd describe the two cells via their vertex indices as
0 1 3 4
1 2 4 5
but if you want a mismatched edge orientation, you'd just use
0 1 3 4 (thinks the edge is 1->4)
5 4 2 1 (thinks the edge is 4->1)
or any other rotation. A good test would probably just try all 2x4 possible
rotations of both cells.

2/ How do you check whether things are right or not? There are many aspects,
but one I've recently thought would be useful to try out is FEInterfaceValues.
Try this: assign a DoFHandler with a continuous but possibly high-order
element with multiple DoFs per edge. Create a solution vector and give its
elements random elements. This should correspond to some finite element field,
and because it is associated with a continuous FE, the jump across the common
edge should be zero.

So create an FEInterfaceValues object and check that indeed the jump of the
solution (or, more simply, of each shape function) at every quadrature point
along that face is really zero. If it is not, there's a bug somewhere.

3a/ For FE_RaviartThomas, the jump itself doesn't have to be zero, but the
normal component has to be. That's also easy to check in the same way as
above: Use FEInterfaceValues, compute the jump in shape functions or the
solution, compute the dot product with the normal to the face, and check that
it's zero. If it isn't, there's a bug.

3b/ For FE_Nedelec, it's the tangential component of course.

4/ Repeat these steps for all 2x4 possible rotations of the two cells. In 3d,
there are even more one can and should try. Little test programs would just
loop over all possible rotations -- I think that GeometryInfo has ways to
rotate cells, though I'm not entirely sure.

5/ You can also output the solution that corresponds to the solution vector as
chosen above, and visually inspect whether what you get is continuous/has a
continuous normal/tangential component.


My best guess is that with such small test programs, you'd uncover cases that
don't work and that could, because they don't involve solving whole problems
or more complicated set ups, be easier to debug and/or fix. Everything that
turns out to work as expected should then be put into the test suite as soon
as you've verified that it works -- we accept small additions to the test
suite pretty regularly, so feel free to open lots of pull requests in this regard.

As mentioned above, regrettably, we've lost collective knowledge about these
issues over the years, and that translates into long response times :-( I'm
sorry you get caught up in this. I'd love to have fixed this long ago, and
would have loved to help you more with it as well, but time availability is
not always on my side :-((

Best
Wolfgang

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

[Konrad Simon]

Dear Wolfgang, 

 

Fundamentally, the place to look is to search through the finite element
classes where they ask for face and line orientations. It is possible that you
don't actually have to do much there: Finite element classes generally work in
the local (cell) coordinate system, where you enumerate shape functions as
0...dofs_per_cell. So when, for example, you compute the local matrix, it
really is local. The interesting part happens when you translate what *global*
degrees of freedom correspond to the local ones, i.e., when you call
cell->get_dof_indices() (==DoFCellAccessor::get_dof_indices). This is where on
every cell has to collect which DoFs live on it, and that has to take into
account the orientations of lines and faces come into play. The cell DoF
indices are cached, but you can see how this works in
dof_accessor.templates.h, starting in line 902 using the process_dof_indices()
functions.

Yes, I found that function. It requires though, that two certain tables

 `adjust_quad_dof_index_for_face_orientation_table` and `adjust_line_dof_index_for_face_orientation_table`

are implemented in derived classes. This is the case for FE_Q but not for any element derived from FE_PolyTensor.

For the DG elemements this is not problematic since all dofs are interior. But for H(div) conforming elements

(FE_ABF, FE_RaviartThomas, FE_BDM) and for  H(curl), which is essentially only FE_Nedelec, it is problematic.

The issue is even worse than just the signs: Dofs are permuted for higher order and signs can switch depending on the permutation.

 (NedelecSZ btw does not inherit from FE_PolyTensor and has a dynamic way to fix the orientation issue but can not be used adaptively.)

What this means is that if you have two cells that disagree about the relative
orientation of their shared face, they when you call cell->get_dof_indices()
on these two cells, you'd end up with a different ordering of the DoFs on the
face. My usual approach to check this would be to create a mesh with two cells
so that one has a non-standard face orientation, assign a Q3 element to the
DoFHandler, and check that that is really what happens. (I generally turn
these little programs into tests for the test suite, just so that what I just
verified to be correct -- or not! -- is actually checked automatically in the
future. I would very much love to see you submit such tests as well!)

This is exactly what I was doing and exactly the way I tested it. :-) 

So this is the process for FE_Q, where there is really not very much to do, or
in fact for all interpolatory elements. The situation becomes a bit more
difficult for things such as the Nedelec element (also the Raviart-Thomas one)
where the shape functions are tangential to the edge. Let's say you have the
lowest order Nedelec element. Then there is just one DoF on each edge, and
both adjacent cells agree about this. But they also have to agree *which
direction the vector-valued shape function should point to*, and that's
trickier. I can't remember how we do this (or whether we actually do this at
all). I *think* that the RT element gets this right in 2d, but I don't
remember about the 3d case. Again, it would probably be quite helpful to write
little programs and turn them into tests.

I know or strongly suspect that there are bugs for RT and Nedelec elements on
meshes where faces have non-standard orientation. I've been thinking for years
how I would go about debugging this, but I've never had the time to actually
do it. But I'll share what I think would be the tools to do this: 

I was already working on a static fix. First, I fixed RaviartThomas. BDM and ABF can be taken care of in a similar static fashion.

For this I implement the above mentioned tables puls an additional table for the sign switches. I needed to change the implementation

of all elements but I do not break legacy behavior (apart form the fixed RaviartThomas) or the interface.

I actually got help from two developers/maintainers, see PR #11414.

 

1/ Creating meshes with non-standard orientations: You can feed
Triangulation::create_triangulation() directly so that the faces have
non-standard orientations with regard to each other. (See step-14 for how to
call this function.) In 2d, this would work as follows: Think of a mesh like this:

3---4---5
| | |
0---1---2

You'd describe the two cells via their vertex indices as
0 1 3 4
1 2 4 5
but if you want a mismatched edge orientation, you'd just use
0 1 3 4 (thinks the edge is 1->4)
5 4 2 1 (thinks the edge is 4->1)
or any other rotation. A good test would probably just try all 2x4 possible
rotations of both cells.

Yes, I submitted another pull request containing exactly this.

 

2/ How do you check whether things are right or not? There are many aspects,
but one I've recently thought would be useful to try out is FEInterfaceValues.
Try this: assign a DoFHandler with a continuous but possibly high-order
element with multiple DoFs per edge. Create a solution vector and give its
elements random elements. This should correspond to some finite element field,
and because it is associated with a continuous FE, the jump across the common
edge should be zero.

So create an FEInterfaceValues object and check that indeed the jump of the
solution (or, more simply, of each shape function) at every quadrature point
along that face is really zero. If it is not, there's a bug somewhere.

Thanks for the hint with  FEInterfaceValues. This I did not think about but probably this is a good think for the test suite. 

I used visual debugging, a modified step-20 on problematic meshes and applications with a known solution. I could reproduce

expected convergence rates for RaviartThomas.

 

My best guess is that with such small test programs, you'd uncover cases that
don't work and that could, because they don't involve solving whole problems
or more complicated set ups, be easier to debug and/or fix. Everything that
turns out to work as expected should then be put into the test suite as soon
as you've verified that it works -- we accept small additions to the test
suite pretty regularly, so feel free to open lots of pull requests in this regard.

This is how I fixed RaviartThomas. Works nicely. 

As mentioned above, regrettably, we've lost collective knowledge about these
issues over the years, and that translates into long response times :-( I'm
sorry you get caught up in this. I'd love to have fixed this long ago, and
would have loved to help you more with it as well, but time availability is
not always on my side :-((

I actually need all these vector elements so we both have an interest to fix it. 

Looking at the current implementation it seems like the developer who did it had good ideas but left the

work in the middle (there are many TODO comments and missing implementations).

There is a lot of work to do. Maybe I can help a bit here.

Thank you a lot, Wolfgang, for your elaborate answer :-)

Cheers,
Konrad

[Wolfgang Bangerth]

On 1/30/21 10:35 AM, Konrad Simon wrote:
>
> I actually need all these vector elements so we both have an interest to fix it.
> Looking at the current implementation it seems like the developer who did it
> had good ideas but left the
> work in the middle (there are many TODO comments and missing implementations).

I think that's a fair implementation. I believe that I implemented the first
version of these elements, possibly even before we could deal with meshes with
non-standard orientation. I don't know if anyone else every looked at this
area later, but it's certainly true that I've lost most of the knowledge I had
when I wrote that first version and never found the time or motivation to read
up on it again. It's been a long time since I've used RT elements myself, so
there was little cause to spend a week or two trying to understand and fix
these issues.

> There is a lot of work to do. Maybe I can help a bit here.

We would certainly love for you to take this area over. It is certainly up for
grabs :-) Let me know where you have questions. I'm not sure I can help much
with writing code, but I'm happy to help answer whatever you can't figure out
yourself!

Best
W.