Add missing renumbering in UGGridEntity specialization

The corners have to be renumbered for the multilineargeometry.

Can you add/improve a unit test to trigger the problem?

Added 1 commit:

• 6a57b936 - [test-ug] add test for corner ordering

Compare with previous version

I added a simple test demonstrating the problem. The test checks the corners of the face ({0,0,0},{0.25,0,0},{0,0.25,0},{0.25.0.25,0}) of the first hexahedron (see hybridtestgrids.hh, line 295) by checking all cube subEntities. Without this patch, the ordering is ({0,0,0},{0,0.25,0},{0.25,0.25,0},{0.25,0,0}) which gives a nonlinear mapping from the reference element. With this patch, the ordering is ({0,0,0},{0,0.25,0},{0.25,0,0},{0.25.0.25,0}), which is still not exactly what I was expecting (but maybe it is ok for subEntities to have a different ordering than the father entity?) but at least gives a linear mapping.

The test makes assumptions on the orientation of subentities. However, there are no guarantees on this and it is in general implementation defined. Why do you expect a special vertex ordering?

I could not think of a better test, it is just to show the problem. Without the renumbering, the mapping from the reference element to the subEntity is nonlinear. I do not know if this is ok in general, but if we can make it linear with the renumbering, why not do it?

What do you mean by "the mapping is nonlinear" and which mapping exactly are you talking about?

The local-to-global mapping from the reference element to the actual element. The cube reference element in Dune here has the corners (number: local coordinate)

0: (0,0)

1: (1,0)

2: (0,1)

3: (1,1)

The face element in the test has the coordinates

0: (0,0,0)

1: (0,0.25,0)

2: (0.25,0,0)

3: (0.25,0.25,0)

So the local-to-global mapping is

(x,y) --> (0.25y, 0.25x, 0)

Without the renumbering, corners 2 and 3 are swapped, so the (multilinear) mapping becomes

(x,y) --> (0.25y, 0.25(x+y - 2xy), 0)

This is a common misconception. As I said, there is no guarantee on the orientation. It is especially allowed that the embedding is different from the reference elements. I.e. ReferenceElement::subEntity() is in general not consistent with the corners numbering of subentities in grid entities. It should only be used as a tool to enumerate the sub-subentities.

See, e.g, flyspray/FS#980 (closed) and flyspray/FS#818 (closed) for detailed discussions about this topic.

So even if the subEntities GeometryType is (GeometryType::cube, 2), it can have a reference element different from the Dune cube reference element? I did not know that. But then the numbering of the corners of the subEntity is arbitrary? Then we could as well do the renumbering :).

I am currently iterating over the codim cc entities of the UGGrid by iterating over the codim cc subentities of the codim 0 elements, which seems to work with the renumbering in this patch. Without the renumbering, there are problems with the jacobianTransposed() method due to the nonlinear mapping.

So even if the current behavior is not wrong, maybe we can still do the renumbering as it does not matter for correctness and provides an (inofficial) way to iterate over the codim cc entities of the UGGrid?

But this "inofficial" way is never guaranteed. When you switch to a different grid, it can be different... it can even be changed without prior notice at any time. You should really not rely on implementation details.

If you have the reference element r_e of an element e and the reference r_s of a subentity s of e, then r_e and r_s are both standard Dune reference elements (e.g. for cubes). But the embedding of s into e is not guaranteed to be consistent with the "canonical" embedding of r_s into r_e induced by r_e.subEntity(...) but is an implementation detail which may depends on e and s.

Even if your patch seems to imply this consistency for your example this is a special case, because there is in general no consistent numbering. That's why the implementation is free to pick what even it wants. Just think of the following "algorithm" to create such a numbering: You recreate your grid by inserting elements one-by-one. Then any new element must be consistent with the face numbering of all its already inserted neighbours. This may lead to the case that you cannot add the next element without violating the rules enforced by the reference element.

If you really need to know, which number the i-th vertex of a face s has within an element containing s you can derive this by looping over all vertices and comparing global indices obtained by indexSet.subIndex(...) (cf. flyspray/FS#818 (closed)). However, you will run into the problem that UGGrid does not implement indexSet.subIndex(...) for entities other than elements or vertices (cf. !95 (merged)).

That said, this question is completely unrelated to the question if the geometry mapping is linear or not (,even if it is still unclear, which map you are talking of). This holds true for the geometry map of both, the element and the face. If the map is nonlinear for a parallelepiped this must be a misuse of or bug in the geometry implementation.

I just discussed with @tobias.leibner and it seems that

1. Tobias had this misunderstanding about the guarantees w.r.t. numbering
2. there is really a bug, as the geometry obtained from the subentity yields a mapping which doesn't have a positive determinant of the jacobian.

We agreed that he will update the tests and check whether this is really the case.

@tobias.leibner: I'm sorry, I missed your explanation on the map above. Now I understand your problem. The opposite/neighbouring relation of the cube corners in the reference element and the entity is different. I.e. it's indeed a misuse of the geometry.

Added 1 commit:

• c8209a2d - [test-ug] make test not relying on any ordering assumptions

Compare with previous version

As discussed with @christi, here is a new test that does not rely on the ordering of the corners. It checks whether the local-to-global mapping for the codim 1 entities is bijective. As we did not find a simple solution to test all codim 1 entities the test checks only the codim 1 entities that lie in the x-y-plane.

Sorry for not intervening. Yes, there is a bug, but you really only see it for element faces with four corners. There, the UG numbering swaps to vertices when compared to the Dune numbering, and that leads to the singular jacobian that Tobias is seeing. For all other cases Carsten's comment applies: the Dune grid interface does not really require the vertex ordering of these geometries to be related to the rest of the grid at all.

I think the second part of the fix is correct, but the first part is wrong. I'll comment there directly.

Thanks for writing the test, but a test only for UGGrid is not of much use. Thinks like this may as well be wrong in other grid implementations (current or future), and a test should apply to all of them. I recognize that the issue is a but tricky to test for.

Added 4 commits:

• a78419d3 - [test] move codim 1 subEntity mapping check to gridcheck
• ac90c804 - [utility.parmetisgridpartitioner] fix compile error
• 516102e9 - [gridcheck] do not use assert
• b07c5c11 - Revert "[uggrid.uggridentity] another missing renumbering"

Compare with previous version