Discovering the capabilities of t8code for my needs

[CsatiZoltan]

Feel free to move this issue to the Discussions; I opened an issue as I cannot open a discussion thread.

I collected a lot of materials (theses, presentation, etc.) about t8code, and I think it would suit my needs. However, before settling with this code, I would like to make sure that it is indeed the right tool.

Goal

A high-order solution exists on a high-order mesh. We want to replace the high-order solution and the high-order cells with an approximate solution on linear cells. Hence, we want to use AMR for post-processing, not for solving.

Characteristics of the problem

• Large unstructured mesh up to a few billion elements
• Mesh is distributed on several processors
• Full quad (2D) or full hexa (3D)
• Non-isoparametric elements: high-order element geometry with a high-order solution (e.g. DGSEM)

Requirements

• Can handle large number of trees (one tree for one original high-order element, I guess)
• Non-duplicate nodes on a subdomain, i.e. a "standard" finite element mesh with a single field value associated with a given node.
• At the end of the AMR procedure, be able to create the explicit mesh (essential for post-processing) based on the initial coarse mesh and the space filling curve
• Geometrical information available to the callback function: for a given forest leaf, I need to know where it is located in the real space so that I can perform the interpolation of the high-order solution. That will serve as an indicator whether to refine the element or not.

Nice to have

• Hybrid meshes that may arise in the future
• Elimination of hanging nodes
• Hanging node elimination not by cutting all the quad/hexa cells to tri/tetra cells, but by keeping as many quad/hexa cells as possible and only including a few transition cells
• Global connectivity across subdomains (i.e. all trees are merged so there will be no duplicate nodes even on subdomain boundaries)
• Data structure that can use existing memory to avoid data copy whenever possible.

---

Based on the demands, is t8code an appropriate tool for me? If not, could you please point me to alternatives?

Thanks,
Zoltan

[jmark]

Hello @CsatiZoltan !

Welcome to t8code!

At a first glance, your given goals, problem description and requirements are a nice match with the capabilities t8code offers. For the nice-to-have section, I can say we already support hybrid meshes. Resolving hanging nodes exactly in the manner as you described is work-in-progress but will be part of t8code in the midterm future. About global connectivity across subdomains and in-place datastructures is, to my knowledge, not available as of now but, of course, nice to have.

[jmark]

Let us know if we can be of help in your project! Do not hesitate to ask if you have more questions!

[CsatiZoltan]

Hi @jmark,

Thank you for the rapid and reassuring reply.

I would certainly have questions, but those are more specific to certain steps of the AMR process. I thought about restricting this issue to the features of t8code, while opening new ones concerning the concrete workflow. Moreover, since these will not be developer's issues, you could give me right to open discussion threads. What do you think?

[jmark]

@CsatiZoltan I converted this issue to a discussion.

[CsatiZoltan]

About global connectivity across subdomains and in-place datastructures is, to my knowledge, not available as of now

But conformity within a subdomain is possible with t8code e.g. by tetrahedralization?

---

When you create an adaptive mesh for a high-order geometry, is the resulting mesh watertight? I cannot see how it can be achieved for non-uniform quad/hex meshes. For visualization, small holes or overlapping elements are not a problem, but for many post-processing filters a manifold mesh is expected.
I see overlappings and holes in the leftmost figure of page 15 in [1] and also in the bottom figures in [2].
I was thinking about it previously; the only solution I see currently is to introduce transition elements between the different refinement levels, resulting in a hybrid mesh.
How exactly does AMR work in t8code for high-order meshes?

---

Resolving hanging nodes exactly in the manner as you described is work-in-progress but will be part of t8code in the midterm future.

Until then, is it possible to turn every quad/hex into conforming tri/tetra? Another idea: what about defining a new subdivision pattern for quad/hex so that they are first decomposed into tri/tetra and then continue their subdivision as usual?

---

Could you please point me to the relevant parts of the code, which allows me executing the following algorithm:

1. give a high-order mesh (in practice, usually quadratic) to initialize the forest
2. for each adaptation step:
   1. obtain the geometry of the linear subcells, described by the tree leafs (in other words, be able to locate the subcell in the physical space).
      It should be possible according to [2]: "When the geometrical position of a fine mesh element or a point inside such an element is required for computation, we compute it using the available geometry information of the coarse mesh.". You rely on Open CASCADE as a geometry kernel. Is there an alternative; I would not like to introduce a heavy dependency? Don't you think that using the ordered set of nodes of the high-order element is enough? I mean, using those nodes with the Lagrange basis functions establishes the geometrical mapping between the high-order reference element and that element in the physical space (i.e. part of the mesh). Another issue is that unless it is you who generates the mesh, you don't have access to the geometrical description (Gmsh's .geo file, a CAD file, etc.). All you have is a mesh and a node ordering convention (Gmsh, Abaqus, CGNS, etc.).
   2. once you have the explicit representation of a subcell, be able to invoke a callback function, which interpolates the high-order solution living on the original high-order element (root of the tree).
   3. the callback function returns a True/Flag based on the interpolated value.
   4. the AMR refines the leafs that returned True.
3. in the end, merge the trees within a subdomain to obtain a mesh with non-duplicate nodes.
4. convert the tree-based description to an explicit unstructured mesh

References

[1] Constructing a Volume Geometry Map For Hexahedra With Curved Boundary Geometries - presentation
[2] Constructing a Volume Geometry Map For Hexahedra With Curved Boundary Geometries - tech report

[jmark]

About global connectivity across subdomains and in-place datastructures is, to my knowledge, not available as of now

But conformity within a subdomain is possible with t8code e.g. by tetrahedralization?

We discussed this concept within the group but there is no implementation of this in t8code. We have an experimental branch in t8code
with a scheme that resolves hanging nodes in a quadrilateral mesh by inserting conforming triangles as a transition layer between different levels of refinement. Right now we are also working on a 3D version of it. But it takes time (and money) till production readiness.

When you create an adaptive mesh for a high-order geometry, is the resulting mesh watertight? I cannot see how it can be achieved for non-uniform quad/hex meshes. For visualization, small holes or overlapping elements are not a problem, but for many post-processing filters a manifold mesh is expected. I see overlappings and holes in the leftmost figure of page 15 in [1] and also in the bottom figures in [2]. I was thinking about it previously; the only solution I see currently is to introduce transition elements between the different refinement levels, resulting in a hybrid mesh. How exactly does AMR work in t8code for high-order meshes?

The examples you referenced are mere visual glitches. When the input mesh is watertight then t8code will not introduce any leaks by its AMR operations on the given mesh. Maybe some words for clarification. t8code is not a classical mesh generator. It reads conformal input meshes, for example from Gmsh, and builds a so called "coarse mesh". In each of these coarse elements an (oct-)tree is initialized which we use for AMR, parallelization, load balancing, partitioning, etc.

Resolving hanging nodes exactly in the manner as you described is work-in-progress but will be part of t8code in the midterm future.

Until then, is it possible to turn every quad/hex into conforming tri/tetra? Another idea: what about defining a new subdivision pattern for quad/hex so that they are first decomposed into tri/tetra and then continue their subdivision as usual?

In general yes, t8code offers the modularity and the interfaces to implement such a scheme you envisioned.

Could you please point me to the relevant parts of the code, which allows me executing the following algorithm:

1. give a high-order mesh (in practice, usually quadratic) to initialize the forest

2. for each adaptation step:
   
   1. obtain the geometry of the linear subcells, described by the tree leafs (in other words, be able to locate the subcell in the physical space).
      It should be possible according to [2]: "When the geometrical position of a fine mesh element or a point inside such an element is required for computation, we compute it using the available geometry information of the coarse mesh.". You rely on Open CASCADE as a geometry kernel. Is there an alternative; I would not like to introduce a heavy dependency? Don't you think that using the ordered set of nodes of the high-order element is enough? I mean, using those nodes with the Lagrange basis functions establishes the geometrical mapping between the high-order reference element and that element in the physical space (i.e. part of the mesh). Another issue is that unless it is you who generates the mesh, you don't have access to the geometrical description (Gmsh's .geo file, a CAD file, etc.). All you have is a mesh and a node ordering convention (Gmsh, Abaqus, CGNS, etc.).
   2. once you have the explicit representation of a subcell, be able to invoke a callback function, which interpolates the high-order solution living on the original high-order element (root of the tree).
   3. the callback function returns a True/Flag based on the interpolated value.
   4. the AMR refines the leafs that returned True.

3. in the end, merge the trees within a subdomain to obtain a mesh with non-duplicate nodes.

4. convert the tree-based description to an explicit unstructured mesh

I think you can implement this by sub-classing t8_geometry. For each tree you assign geometry data encoded as coefficients of a Lagrange polynomial. I recommend to check out this header https://github.com/DLR-AMR/t8code/blob/main/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx and associated source files. It builds a squared disk mesh by hand and defines the curvilinear geometry of it.

You can try this specific example with this example program https://github.com/DLR-AMR/t8code/blob/main/example/cmesh/t8_cmesh_geometry_examples.cxx .

We're in the process of reworking the API a little bit. So I recommend to check out the two PRs as well. The second PR implements a triangulated spherical shell.
#753
#755

This would cover the first point. Merging the subdomains and converting it to an unstructured mesh ... for that I need more information what your actual goal is here. Which format should the unstructured mesh be in, for example?

[Davknapp]

Hey @CsatiZoltan , thank you for your interest in t8code! @jmark is right, we can transform a forest into a vtkUnstructuredGrid in parallel. Each process writes its own vtu-file and one process aditionally writes a .pvtu file.

[CsatiZoltan]

Hi @Davknapp, this is great, thank you. What is the process of writing a forest into a vtkUnstructuredGrid? At some point, you need to provide the geometrical information as well (since the forest does not store the geometry). It would make sense that the forest is written to the .vtu file element by element. In that case, the elements (and the fields defined on the elements or on their nodes) are independent. What I want is no duplicate nodes in the unstructured grid. Can it be done from within t8code? If not, is there a filter in VTK or ParaView or in an external library that can get rid of the duplicate nodes?

[Davknapp]

To write the coordinates of the point we evaluate the geometry at the vertices of each element. If a geometry need more points than that (for quadratic elements for example), we use more points than that.
Currently points are duplicated, but there is a vtk-filter to get rid of duplicated points. We already have an Issue (#609) for that, but haven't added it to the code yet.

[CsatiZoltan]

If a geometry need more points than that (for quadratic elements for example), we use more points than that.

What I would like is the vertices of the forest elements, i.e. the nodes of the leaf elements in the physical space. For that, I am thinking of the following:

1. Finish the AMR, i.e. refine until the desired tolerance has been reached.
2. For each leaf element
   1. Find the tree root it belongs to.
   2. The tree roots were created for the original mesh elements. Therefore, we know the original (possibly high-order) element.
   3. Find the (4 for quad, 8 for hexa) vertices of the leaf element on the reference element (square or cube).
   4. The geometrical mapping is given by the nodes of the original/coarse element (found in step 2.2) and the Lagrange basis functions. Use the geometrical mapping to map the 4 or 8 vertices of the leaf element from the reference element (found in step 2.3) to the physical space.
   5. Interpolate the high-order solution in the 4 or 8 vertices of the leaf element on the reference square/cube.
   6. Create a VTK cell of type VTK_QUAD or VTK_HEXAHEDRON and add it to the vtkUnstructuredGrid object.
3. Write the vtkUnstructuredGrid object to a .vtu file.

The delicate issue is that the AMR partitions are not the same as the partitions used by the solver due to the partitioning of the space filling curve. So it can happen that the coarse element corresponding to the tree root is on another processor than the lead element. If the initial mesh is scattered on all processes (as is done in p4est if I recall correctly), this is not an issue.

Do you think this workflow is realizable with t8code (even if I have to write some code myself)?

Once I get there, I can try resolving #609.

[Davknapp]

You can have a look at t8_forest_vtk_write_file_via_API, and you'll see what we are currently doing. You cann call this function with the forest you created in your step 1.
We then iterate over every tree and every element in the tree and construct a VTK_Cell from the element. Currently we support linear and quadratic VTK_Cells.
The points are mapped to the physical space with respect to the geometry used.
All these cells are stored in a 'vtkUnstructuredGrid', which we then write into a .vtu file.

[jmark]

   * utilities to interrogate the geometrical position of a leaf?

In t8code you can implement new geometries and have full access to geometrical position of leaf elements which allows to implement custom geometries.

   * handling curved elements?

t8code handles curved elements via its geometry classes. I posted examples how to do this in a previous comment. Also check out https://github.com/DLR-AMR/t8code/blob/main/example/geometry/t8_example_geometries.cxx.

   * readers for common mesh formats (Gmsh)?

t8code can read in Gmsh files v2 and v4.

   * advanced output (see the next question)?

As already discussed, VTK output.

2. I do not understand yet how a high-order coarse element is refined. I thought that the AMR replaces the coarse element with a collection of linear elements. But [apparently](https://github.com/DLR-AMR/t8code/wiki/Feature---Curved-meshes#exporting-curved-meshes) you obtain curved elements. What I would like as an output is a collection of linear elements that approximate the geometry of the coarse element. In other words, in the limit, you cover the coarse element. For this feature, you do not need any CAD engine, just the (high-order) input mesh. This is how the post-processing module of Gmsh performs the adaptation. Yes, I could simply extract the corner vertices of the high-order forest elements to obtain the linear elements, but that is an overkill: needless temporary storage and dependence on Open CASCADE.

t8code does not store or refine any geometry of the elements. It computes the geometry of a leaf element on demand. The exact or high order geometry information is stored in a derived custom geometry class (or within OpenCascade read in from a *.brep file.) t8code only stores the physical vertices of the coarse elements for simple (linear) VTK visualization.

Do not hesitate to follow up on this issue and ask further.

[CsatiZoltan]

@jmark To interactively discover t8code, I would like to use the Julia wrapper t8code.jl. Does it expose all the functionalities of t8code that I need? In other words, can I write e.g. the geometry extension discussed above in Julia (I do not want to rely on Open CASCADE) ?

[jmark]

Yes, T8code.jl exposes the whole API of t8code as Julia functions.