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 <https://github.com/dealii/dealii/pull/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 -- 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/652cf426-6b5d-4f3f-bfce-f1266e171d20n%40googlegroups.com.
