[Dune] Non-unique coordinate mappings of Codim(1) entities?

[Christian Waluga]

Hello everybody,

I am trying to implement  a hybrid DG method using DUNE and I have a problem with the mapping feature of the intersections. In addition to my element DOFs, I have an additional function living on the trace of the grid, i.e. I have some DOFs distributed on the faces of the mesh. The intersections of the elements seem to be associated with the codim 1 entities of the mesh, but the orientation of their geometries seems to be not unique in my experiments (using ALUSimplexGrid<2,2>). 

Okay, this gets a bit hard to explain, I will make an example:

Lets say I have 2 neighboring elements e1 and e2. I get the face entity associated with the intersection by calling something like this...

GV::Codim<1>::EntityPointer e1 = is->inside()->subEntity<1>(is->indexInInside());

I do the same from the other side for comparison.

GV::Codim<1>::EntityPointer e2 = is->outside()->subEntity<1>(is->indexInOutside());

If I compare the local to global mappings of e1->geometry() with the mappings from e2->geometry() (which is the same but from another elements perspective), I observe that they don't seem to match. Lets say if the local coordinates 0.25 and 0.75 on e1 are mapped to (0.25,0) and (0.75,0) in the physical domain, then the geometry of e2 maps them to (0.75,0) and (0.25,0), which messes up the integration, since my DOFs on the face entities are in a specific order (I have integrals like \int_{\partial T} {L*V ds}, where L lives on the trace of T and V is a discontinuous test function that is H1 in each element T).

I worked around this by replacing e1 with e2, in case I already visited the face, to get consistent mappings on both sides. But this workaround is very ugly, since it needs an array to remember which faces were already visited and since I cannot work directly with the intersection.geometry() to do the mappings (right now, I do the coordinate mapping in the integration procedure by mapping from codim 1 reference coordinates to physical coordinates and back to coordinates on the codim 0 reference element to map the integration points that lie on the face to the reference element geometry).

Since I can also directly iterate over codim 1 entities I would assume, that their geometry is unique, regardless from which element I get it by calling subEntity. I would try it the other way round, i.e. iterate over the codim 1 entities and add the contributions from both intersecting elements, but I know of no elegant solution to get references to these intersecting elements.

I can write a minimal code example, that demonstrates this behavior, if needed.

Maybe I got something wrong here, so I am thankful for any feedback!

Best, Christian