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: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
--
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 the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c0c60681-4e4f-f4e3-c4ac-6e6c8978e7f0%40colostate.edu.
