precice / openfoam-adapter

OpenFOAM-preCICE adapter
https://precice.org/adapter-openfoam-overview.html
GNU General Public License v3.0
135 stars 80 forks source link

Reimplement connectivity definition #258

Closed davidscn closed 1 year ago

davidscn commented 1 year ago

https://github.com/precice/openfoam-adapter/blob/2a51e6e24bc36776051a06566413bc66cbe181bd/Interface.C#L309-L316

See https://github.com/precice/precice/issues/1248#issuecomment-1103875457, we need to find a new solution here.

MakisH commented 1 year ago

To summarize:

Current implementation

The affected function is preciceAdapter::Interface::configureMesh in Interface.C. What this does is:

  1. Gets a list of faces and coordinates at the interface:

https://github.com/precice/openfoam-adapter/blob/7afeaede9d9aec8d48457c77e8be9c6131897f44/Interface.C#L271-L273

  1. Modifies the coordinates to account for any deformation
  2. Uses Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const to triangulate each face and get back a triangle:

https://github.com/precice/openfoam-adapter/blob/7afeaede9d9aec8d48457c77e8be9c6131897f44/Interface.C#L293

  1. Uses the preCICE getMeshVertexIDsFromPositions to populate an array of preCICE vertex IDs for each triangle (triVertIDs):

https://github.com/precice/openfoam-adapter/blob/7afeaede9d9aec8d48457c77e8be9c6131897f44/Interface.C#L309

  1. Passes each triangle edge to the preCICE setMeshTriangleWithEdges:

https://github.com/precice/openfoam-adapter/blob/7afeaede9d9aec8d48457c77e8be9c6131897f44/Interface.C#L316

What changes in preCICE v3

The discussion in https://github.com/precice/precice/issues/1248 concluded in removing getMeshVertexIDsFromPositions, since the implementation is complicated and unreliable. We need to directly give vertices to preCICE.

There has also been a series of changes, most recently https://github.com/precice/precice/pull/1322:

Respective changelog entry:

  • Removed edge IDs from the API: setMeshEdge now returns nothing, all setMeshXWithEdges variants were removed from the API, and setMeshTriangle|Quad now accepts vertex IDs instead of edge IDs.
  • Added efficient bulk functions setMeshEdges|Triangles|Quads|Tetrahedra for setting connectivity based on an array of vertex IDs.
  • Changed API to guarantee the creation of hierarchical connectivity. For example, setMeshTriangle now creates necessary edges internally.

This means that we can no longer give locations to setMeshTriangleWithEdges (+ renaming to setMeshTriangle).

Something else that has happened in the meantime, is that we now have setMeshQuads. Since we mostly work with hexahedral meshes in OpenFOAM, we could directly extract faces from each OpenFOAM patch and give the points to preCICE with setMeshQuads. This does not require triangulating the faces or keeping track of locations.

What am I overlooking here?

P.S.: See also the porting guide, which probably needs some example on this.

fsimonis commented 1 year ago

This means that we can no longer give locations to setMeshTriangleWithEdges (+ renaming to setMeshTriangle).

Not sure if I understand what you are saying. The API change for setMeshTriangleWithEdges is a rename only, hence trivial.

Something else that has happened in the meantime, is that we now have setMeshQuads.

Be aware that these quads must be planar. This requirement is also checked, so triangulation in the adapter will most likely be more efficient.

MakisH commented 1 year ago

This means that we can no longer give locations to setMeshTriangleWithEdges (+ renaming to setMeshTriangle).

Not sure if I understand what you are saying. The API change for setMeshTriangleWithEdges is a rename only, hence trivial.

Indeed, I was confused by the diff rendering in https://github.com/precice/precice/pull/1322, mixing up setMeshTriangle with setMeshEdges. What changes is that we can no longer give vertices to preCICE and let it associate them to edges, but we have to give edge IDs directly. I don't know how to get those IDs, then. What is the recommended approach? Call setMeshEdge explicitly? It is not so common to know how points are connected.

In that case, maybe we can work with whatever the faceTriangulation class of OpenFOAM offers (returned by faceTri).

For completeness, this is where we call setMeshVertices.

Something else that has happened in the meantime, is that we now have setMeshQuads.

Be aware that these quads must be planar. This requirement is also checked, so triangulation in the adapter will most likely be more efficient.

That's the comment that I was looking for, thank you. I don't know if we can assume that the quads are always planar. I think that snappyHeshMex can also create curved faces. But maybe it is an assumption that can cover most of the cases.

fsimonis commented 1 year ago

You essentially need an association (or map) from OpenFOAM points to preCICE vertex IDs.

Currently, you are using preCICE for this purpose with getMeshVertexIDsFromPositions, which may lead to surprising results.

I don't know the internals of OpenFOAM. Maybe you can label/tag points somehow, then you could use the preCICE vertex ID as a label. Later when you iterate over the faces, you could then extract the labels of all defining vertices and pass them to preCICE as vertex IDs.

Another alternative would be to create your own mapping from coordinate to preCICE vertex ID (possibly per patch/mesh). You can then query coordinates whenever you need the vertex ID of a/the point located at those coordinates. You could write your own comparator and use a std::map, use the boost::geometry::index::rtree to query the nearest-point, or try to reuse the OpenFOAM indexedOctree

MakisH commented 1 year ago

We currently do the following:

  1. Set the preCICE mesh vertices
  2. Iterate over all OpenFOAM faces and ask OpenFOAM for triangle node coordinates
  3. Get the preCICE IDs of these triangle node coordinates
  4. Create preCICE triangles

Even if at step (1) we store a map between OpenFOAM vertices and preCICE vertex IDs, we would need a new mapping after step (2).

@davidscn What if, instead, we first triangulate the faces, and give some duplicate nodes to preCICE? Essentially:

  1. Iterate over all OpenFOAM faces and ask OpenFOAM for triangle node coordinates
  2. Store these coordinates as triplets
  3. Use all of the coordinates to set preCICE mesh vertices
  4. (preCICE hopefully does not define duplicate nodes, but still gives the same vertex IDs for the same duplicate vertices - @fsimonis can I hope on that?)
  5. Store a map from the coordinates to IDs
  6. Set preCICE mesh triangles by parsing the map in triplets
fsimonis commented 1 year ago

preCICE hopefully does not define duplicate nodes, but still gives the same vertex IDs for the same duplicate vertices - @fsimonis can I hope on that?

If you tell preCICIE to define a new vertex, it defines a new vertex. :wink:

Hence the rtree section in the documentation.