Bug Report: Mesh adaptation failure when using CAD geometry

[saustinp]

CAD_mesh_debugging.zip

When adapting a forest defined using the CAD geometry, the resultant mesh is not valid. There are interior holes in the mesh in some regions and overlapping elements in others. This makes it impossible to use the CAD geometry feature in t8code when this occurs.

To Reproduce
I uploaded a minimum example to this bug report which includes everything to reproduce the bug: mesh generation script using the gmsh python API, cmake file and build script for the example, and paraview state file to view the output .VTU file.

Expected behavior
This behavior does not happen when using the linear element geometry, only when using the CAD geometry. Even when using the CAD geometry, I expect that only the nodes on the boundary would move, but it appears that all of the nodes in the element on the boundary (and even in elements away from the boundary) move in response to projection to the CAD boundary, resulting in an invalid mesh.

Environment (please complete the following information):
OS: Ubuntu 24.04
Compiler: g++/gcc 14.3.0 (Ubuntu 14.3.0-12ubuntu124ppa1)
MPI version: OpenMPI 5.0.8
t8code commit hash: f878a8e
CMake options:

-DCMAKE_BUILD_TYPE=Release
-DCMAKE_CXX_COMPILER=g++
-DCMAKE_C_COMPILER=gcc
-DCMAKE_INSTALL_PREFIX=~/.local/t8code
-DCMAKE_CXX_STANDARD=20
-DCMAKE_C_STANDARD=17
-DCMAKE_CXX_EXTENSIONS=OFF
-DT8CODE_ENABLE_VTK=ON
-DT8CODE_BUILD_AS_SHARED_LIBRARY=ON
-DMPI_HEADER_FOUND=TRUE
-DT8CODE_ENABLE_OCC=ON \

Estimated priority

"Priority: high" Should be solved as soon as possible

Additional context
Add any other context about the problem here.

[sandro-elsweijer]

Hey Sam,

thank you for your report!
The holes and overlapping elements are generated due to the fact, that the elements themselves are curved.

In t8code, we apply the curvature on the cmesh trees and not individual elements. We do this, since we would be really susceptible to produce inverted elements at fine refinement levels. On the top left in Fig 1 you can see a triangular tree, with one curved edge. This curvature is projected onto all elements inside that tree (in the real implementation the mapping looks a bit different, but this was easier to draw and has the same concept). If we would only curve the one layer of elements at the curved edge, we would get the top right tree, which consists of invalid elements. And with each further refinement this would get worse.

That is why we curve the complete tree. And if you evaluate the geometry of the curved elements at more points than the corner points (e.g. in a high-order code), then you would get the mapping of the curved elements as you can see in the Paraview Screenshots Fig 2 and 3.

The Problem arises, when the elements are considered linear and hanging nodes appear. In Fig 1 bottom left you can see a uniform curved tree with linear elements. Here, everything is fine. But if we coarsen the top element we introduce a hanging node and arrive at bottom right, the case you observed.

So the question is if your code can support curved elements, which do not have the same problem. Another solution would be to resolve the hanging nodes and lastly, a more complicated approach would be to map hanging nodes onto the edge, they are hanging from. But this would violate the design how geometrical evaluations are designed in t8code. The actual location of nodes should be independent of the actual refinement of the tree.

Does this explanation help you with your problem?

Fig 1

Fig 2

Fig 3

[saustinp]

Hi Sandro,

This is very helpful - thank you for the reply with the graphics. It makes sense to me why you would want to treat the entire cmesh element as nonlinear in the case of a curved boundary to prevent the element inversion that you show in the top right figure.

In my solver workflow, I'm "exporting" the mesh coordinates and connectivity from t8code prior to executing a solver step in order to maintain compatibility with my simulation. I am encountering this issue during the hanging node elimination process, which is a post processing function that I implemented outside of t8code. To eliminate hanging notes, I take two pieces of information: 1) the coarse element that needs to be split into transition elements, and 2) the face of that element that needs to be split (the face that contains the hanging node). I then compute where on the coarse element face the hanging node should be using a linear approximation (I place the new node at the midpoint of the coarse element face) and then search the mesh coordinate array for the matching point, which I then assign to the transition elements that I insert to split up the coarse element into fine elements.

You can see how this breaks down in the case of the overlapping elements/holes we have been discussing, as the hanging node no longer lies at the midpoint of the coarse element face but some distance away due to the nonlinear representation.

My code can certainly handle nonlinear elements, and this would be the preferred solution to overcome this challenge. However, I need your advice on a few aspects:

1. I'm trying to think of a way to export the mesh from t8code with the high order connectivity required for nonlinear elements. Currently, I am computing the mesh connectivity and coordinates by pulling the coordinates at the vertices of the elements using t8_forest_element_coordinate, and then computing the connectivity array by de-duplicating the coordinate array using a hash table. As you can see, because I rely on the coordinates of the element vertices, this is inherently a linear construction. For the curved elements, I need a high order element geometry representation. My proposed procedure for pulling the connectivity is as follows:

1. Query the cmesh element (tree) to see if it lies on a boundary. If it does, then instead of pulling the points using the linear vertices, move onto step 2: 2. Use t8_forest_element_from_ref_coords to pull the coordinates corresponding to the locations of the high order nodes.

2. With the nonlinear geometry representation for the boundary trees, I still need to perform hanging node elimination during refinement. As I mentioned, my current approach is to assume a linear element geometry and find the midpoint on the face of the coarse element that needs to be divided. I suppose that, similar to 1), instead of using the midpoint, I could find the physical coordinate of the hanging node by evaluating t8_forest_element_from_ref_coords at the midpoint of the element face in reference space, which performs the nonlinear mapping for us.

3. I think the native t8code elements only support up to 2nd order, but I don't believe this is a problem for me. If I switch to using t8_forest_element_from_ref_coords for performing all coordinate samplings with the nonlinear elements, I can sample the high order element geometry for all points required by my solver, including along the boundary with the CAD projection feature.

Does this approach make sense and sound reasonable? Please let me know if you'd like me to explain it in more detail.

Thanks,
Sam

[sandro-elsweijer]

Hey Sam!
We actually have a branch implementing 2D hanging node elimination for quads, but it is quite old and needs major updates and the design concept needs to be altered to be merge-able into main. But at least we are working on hanging node elimination in 2D :)

About your implementation: This sounds like a good plan, but I have some remarks/questions:

1. Why don't you always use t8_forest_element_from_ref_coords? With linear elements this will always be a (bi-)linear interpolation and hence probably as fast as your implementation. You would also not have to distinguish between linear and curved trees.
2. You do not have to check if an element/tree is at a boundary. Instead you can directly ask the tree if it uses the linear or curved geometry like so. During the msh import, t8code checks if a tree is linked to a cad geometry and if not, it will use the linear geometry. It will also use the linear geometry if the linked cad geometry is linear like a plane or a line.
3. What do you mean by that?

I think the native t8code elements only support up to 2nd order, but I don't believe this is a problem for me.

t8codes curved cad geometry is exact and queries the cad geometry. So you can be as exact as you want to.

[sandro-elsweijer]

Whoopsie