Incorrect adjacencies in polyhedral mesh

[Johannes Probst]

We discovered that MOAB returns incorrect adjacent cells under certain circumstances. We have observed this only in polyhedral meshes so far.

With this issue I am supplying a minimum working example to reproduce the issue. We create one pyramid (i.e. a polyhedron with 5 faces) and one tetrahedron (a polyhedron with 4 faces) which touch at their base. The vertices are shared such that the pyramid’s base face which consists of 4 vertices also contains the 3 vertices which the base of the tetrahedron is made of.

The image below represents the situation.

‌

When we query the adjacent cells for the base face of the pyramid (= the quad), we get one cell, which is what we expected. However, if we query the adjacent cells for the base face of the tetrahedron, we get two cells: the tet and the pyramid. This is surprising because we would expect only one adjacent cell here as well.

I am able to reproduce this behavior with the current tip of master (d6e6ae2e7040f7259762a1ef6e80d4b3c8fc466e).

Please find attached a small C++ program which reproduces the issue. We compile it on Linux with

g++ -std=c++11 reproduce.cpp -I/opt/moab/include -L/opt/moab/lib -lMOAB -L/opt/hdf5-1.8/lib -lhdf5 -o reproduce

In the attachments you also find the mesh which is generated by the program written out as H5M file.

I would also like to mention that we have repeated this test with 2 tets and we get 1 adjacent cell each, which is what we expect.

Although such a mesh certainly is a bit of an edge case, it is not invalid from our point of view. The base faces where the cells touch actually represent a system boundary on each side, so they are some kind of “infinitesimal wall”. The fact that the vertices are shared on both sides is probably caused by the mesher which merges vertices once they become too close.

[Iulian Grindeanu]

Thanks for reporting this; I looked at the code, and it will be hard to get the behavior you need (or want). In the adjacency query, we use vertex adjacencies as the main driver; So a face is adjacent to a solid if all vertices of the face are adjacent to that solid; A not so immediate consequence of this is , for example, a triangle in the interior of a prism is “adjacent” , if all vertices are part of the prism, even if the triangle is not a face of the prism;
(for a prism abcdef, the triangle formed by abf is “adjacent” to the prism)
the relevant code is in src/AEntityFactory.cpp, lines 1123-1150

In your case, the triangle is part of an actual face of the pyramid, and it is reported as adjacent to the pyramid

My guess is that it will not be easy to get the behavior you want; Maybe you can afford to verify one by one the faces of the polyhedra reported as “adjacent” .

I tried your code with the true pyramid and tetrahedron, and the behavior is the same;

[Johannes Probst]

Hi Iulian

thanks a lot for taking a look and the explanation! I was indeed unaware of this behavior.

There seems to be something I don’t yet fully understand. I have modified the initial test slightly to use 2 tets (instead of 1 pyramid and 1 tet). The situation looks like in the following image:

‌

This time, MOAB gives 1 adjacent cell each. This time around, both cells contain all vertices of the shared face, so I would have expected 2 adjacent cells each time. However, the test tells me something different. So there seems to be something I’m still missing.

I’m attaching the modified code as “reproduce-2tets.cpp” to allow you to reproduce it easily.

As you suggested, a possible workaround is to check the adjacencies returned by MOAB whether or not the faces are actually part of the cell. This is also the fix we’ve implemented inside of our application. We’re treating the adjacencies as “maybe adjacent” and do another check to filter out the ones which don’t match. It comes at a slight performance penalty but we can live with it for now.

Johannes

[Johannes Probst]

Source code to reproduce the adjacency test for 2 tets

[Iulian Grindeanu]

Yes, this seems like a bug; or at least, the behavior is inconsistent ; it should probably be 2 in this case too;

If you modify the second tet to use the first face triangle (c1f1), not the second face (c2f1), you will get 2 polyhedra adjacent to c1f1 , and 0 adjacent to c2f1

The logic is inconsistent, I think;

when there are MBPOLYGON cells sitting on top of each other, that border polyhedra, we are looking for the connectivity list of faces in the polyhedra, for an exact match; when there is only one face (poly or tri, or anything) we look for adjacency of each vertex.

So what do you want to see ? I do not know the reason for having the “equiv_entities” logic in src/AEntityFactory.cpp , separate for polyhedra/polygons; in my opinion, it should be just based on vertices, for simplicity; it seems that you prefer an “exact” match in some cases; I do not like the exact match also because an mbpolygon ABC is different from BCA, or from ACB; we do not expect any orientation of faces, or order in polyhedra;

Maybe it would make sense to orient all faces in a polyhedra to point “outwards” , then we would never share a polygon between 2 polyhedra, all shared polygons would be duplicated (because of orientation?) I do not like that, again.

‌

‌

‌

[Johannes Probst]

Hi Iulian

thanks for your reply!

Our desired behavior would be that MOAB return only 1 adjacent cell in both cases (pyramid and tet, or tet and tet). Let me explain why.

In the following drawing I have sketched up the situation again (hopefully) in a clearer way.

‌

In this case there are 2 faces (F1 and F2) which are made of the same vertices (V1, V2, V3), even in the same order (forgetting about normal vectors for a second).

As mentioned earlier, the faces F1 and F2 represent a border to the mesh on either side, like a wall with zero thickness (which are sometimes called “baffles”).

This mesh is consistent with the general rule that boundary faces have exactly 1 adjacent cell. This is a very strong guarantee because it does not require us to special-case the baffles in any way. This mesh would behave just like any other.

We would expect the adjacency to F1 of dimension 3 to be only P1. Analogously, we expect the adjacency of F2 to be only P2. We would not expect that F1 and F2 are considered “equal” based on the fact that they share the same vertices. We consider 2 face entities to be equal only if the EntityHandle has the same value. I think it can be left to the developer of an application to decide whether a comparison should instead consider the vertices, or the order of the vertices (i.e. the direction of the face normal).

In case the programmer wants something different, for example “the neighboring cells which contain the same vertices as the face”, it is very easy to get with MOAB, using moab::Interface::INTERSECT in get_adjacencies (and passing the 3 face vertices instead of the face). I have modified the program for the 2 tets to do just this, and this makes MOAB return 2 cells each time, just what I was looking for (code is attached). In fact, It seems to me that this is what the INTERSECT operation was made for (but I could be wrong).

Since this workflow allows the programmer to adjust the behavior of MOAB explicitly in both ways, it is the behavior which we would prefer.

---

If we consider the alternative for a moment, where adjacencies are always such that they are based on the vertices, there would be no straightforward way to distinguish a “regular internal face” from a baffle. This forces a developer to query the face connectivity of the polyhedron and check if the face is actually part of them (which it might not be). I find this behavior relatively surprising compared to the other variant discussed above.

I hope I was able to provide some context. Happy to hear your thoughts and continue the discussion!

Johannes

[Vijay M]

@Iulian Grindeanu Has this been resolved ? If yes, close the issue.

[Iulian Grindeanu]

This one slipped under the radar, I missed it; Sorry Johannes Probst ! The code/logic inside moab is not modified at all, but I think Johannes found some workarounds, which are adding overhead;

[Tim Tautges]

The ability to represent multiple faces that have the exact same connectivity is intentional, and IMO important; some of the work I did on local hex element modifications in fact depends on it. Not all mesh databases are capable of doing this, and there was some debate in TSTT/ITAPS about whether cases like that could be called valid (IMO they are valid in some cases).

When you create the first element that uses the first face, you can do everything normally. But when you create the 2nd face that has the same connectivity of the first, you need to also create an explicit adjacency between one of the faces and the element, to tell the code which of the faces of duplicate connectivity is adjacent to that element. When you create the 2nd cell, you also need to set an explicit adjacency between the 2nd duplicate face and the 2nd cell. If you don’t set the explicit adjacencies correctly, when you ask the elements for the adjacent faces, or vice versa, you’ll probably get MB_EQUIVALENT_ENTITIES back as an error code, telling you you have duplicate-connectivity elements (faces) without enough adjacency information to disambiguate the adjacencies properly.

‌