Support nearest-projection mapping

[MakisH]

The adapter currently does not provide preCICE with the mesh edges, but only with the locations of the vertices. This means that currently nearest-projection mapping is not supported.

preCICE provides the following method in its src/precice/SolverInterface.hpp file:

int setMeshEdge (
     int meshID,
     int firstVertexID,
     int secondVertexID );

The mesh vertices are assigned in the adapter's Interface.C file, in the method preciceAdapter::Interface::configureMesh().

[uekerman]

Please note that, besides edges, in 3D also triangles need to be defined. setMeshTriangle

[MakisH]

OpenFOAM can use various types of cells, the most common being hexahedral cells. How should the triangles be defined in hexahedral cells, where we have four-edge faces?

[uekerman]

A possible solution would be to sub-divide the hexa faces into two triangles. We use a similar approach (successfully) for Ateles. The other option would be to use quads in preCICE (setMeshQuad). This feature, however, has (to my best knowledge) not yet been carefully tested.

[MakisH]

There are now thoughts on supporting quads on preCICE: https://github.com/precice/precice/issues/153

This could take a while, though.

[MakisH]

There is a question on whether to define the triangles using the nodes of one face (so... inside the face), or to use the face centers of neighboring cells.

(here, I should add a nice illustration in the near future)

Looking at the big picture:

• We may have hexahedral or tetrahedral cells.
  • Therefore, with quads we would restrict the implementation to hexahedral cells.
• preCICE will use the nearest mesh element to map onto.
  • This could be a triangle, an edge, or a vertex.
  • Therefore, in 2D it is enough to have edges.
  • Therefore, we should define vertices, then edges, and then triangles (from edges).
• We may have 3D or 2D cases (pseudo-2D).
  • Defining triangles inside a cell could cause strange effects on pseudo-2D cases.
  • Therefore, just define edges between face centers.
  • Using the face centers would also work in tetrahedral faces, at least in 3D.
• In practice, the Nearest-Projection mapping is mainly needed in CHT. In FSI, with our current setup, only the Structure solver would need to define connectivity.
  • In general: In a consistent mapping, the from participant needs to have connectivity. In a conservative mapping, the to.

[davidscn]

In general, it would be preferable to create the edges/triangles based on the face centers, since the FVM provides the information there. But we have a loss of connectivity information, when we decompose the mesh. Openfoam crosses our "virtual edges" in the decomposition and the processor is not aware of the neighbor face. The idea is now, to provide the connectivity information based on the real nodes and apply an interpolation from the cell centers to the nodes.

[MakisH]

Our current question is "what is the connectivity the the solver knows" and then "how to use this connectivity in the adapter (i.e. which methods to use)".

Another question is "how to define triangles". Does OpenFOAM offer some triangulation method?

If it has a triangulation method, then we should apply this first. Then, we need to extract the connectivity from the triangulated mesh, which would maybe be easier.

Therefore, I think we should focus on finding a way to triangulate surfaces in OpenFOAM. This will unlock us some doors.

[davidscn]

There is a faceTriangulation() object in OpenFoam and we have all building blocks to triangulate our mesh as desired. I would prefer to use the precice function setMeshTriangleWithEdges, which can be used to define triangles using vertices. Edges are created on the fly by precice.

Although this function is significantly slower than the mentioned function setMeshTriangles, we would circumvent the problem of the right precice own edgeID adressing, since only the vertexIDs are needed for the definition of the triangles.

[MakisH]

After triangulating (the decomposed mesh, I assume) do we still get the problem of edges belonging to different subdomains?

How would you use the faceTriangulation() in our context? If it returns points, we could still use the vertex-based method. But let's make a prototype first.

[davidscn]

Since we define triangles (and thus edges) based on the real nodes and not on the face centers, we will not lose any connectivity information. In the previous step we can interpolate from the faces to the nodes, as already described, using the faceToPointInterpolate() function in OpenFoam.

We apply therefore the triangulation and get three faces back, which are labelLists with the three vertices clamping each face.

[MakisH]

We apply therefore the triangulation and get three faces back, which are labelLists with the three vertices clamping each face.

Then we can also create triangles from edges, using (list[0], list[1]), (list[1], list[2]), (list[2], list[0]) for edges. But again, whatever is easier first!

(I assume you meant "and get two faces back")

[davidscn]

(I assume you meant "and get two faces back")

Yes, of course, a quad, two triangles.

Then we can also create triangles from edges, using (list[0], list[1]), (list[1], list[2]), (list[2], list[0]) for edges. But again, whatever is easier first!

The quick answer: Yes, there might be a solution. But we need to avoid defining edges multiple times. If we define an edge in precice, we get an ID from precice back, which we would need then to define the triangle. By iterating over all faces (e.g. going row by row on the interface patch), we would need at some point the edge ID we defined in the previously (e.g. the previous row on the triangle below). This is difficult to handle. Even if we provide here a solution: I guess the whole process won't be much faster, we just source the worrk out in the adapter.

[MakisH]

Some conclusions from our so far discussions here and in person:

1. In 3D, we triangulate the faces, creating two triangles out of each face
2. Therefore, in 3D, we always need a nodes-based mesh for the solver that writes
3. In 2D it not clear what to do, as we have two options: i. Define the interface edges on one of the two pseudo-faces (this assumption is used in the figure, for simplicity). ii. Define the interface edges between the cell centers. This causes problems when a patch is cut in two by a parallel boundary, as we cannot "bridge" the two neighboring cell centers.
4. In any case decided in (3), Solver A first interpolates the centers-stored fields (e.g. Temperature) to its nodes (therefore, its buffer is of size equal to the number of interface nodes).
5. When Solver B reads, it simply stores the received values to the cell centers, as the interface is defined on its cell centers.
6. Solver A and Solver B are simply swaped when B writes and A reads (e.g. the heat flux).

[davidscn]

In order to realize the described method., each participant needs to define two meshes: One based on the face nodes (which we use to read the data from the participant) and one for the centers (which we use to write the information to the mesh). This allows us, to skip an unnecessary interpolation step during the write procedure, since we directly use the face centers. However, this might be uncomfortable for users, since we need to define two mesh names (one for the centers and one for the nodes) in the precice xml file, which (at the moment) differ from the name specified in the yml file of the adapter (see also the example setup in the tutorial in the pull request #46).

Besides the input data concerning the meshes, the user can specify with "faceTriangles" instead of e.g. "faceCenters", that the adapter should perform in the described way, which means, triangles (connectivity) and the two meshes are defined for each participant with the respective mapping setup. If no nearest-projection mapping is desired, this should allow to skip the computational more expensive step of the triangulation.

[MakisH]

The feature in now implemented in develop, the issue will be closed automatically next time we merge to master.