Mesh duality: m_poly faces not planar nor perpendicular to the corresponding m_tet edge

[Mengesh]

First of all, thank you for this great project! It spared me a lot of trouble during my research.

I was investigating the duality of the meshes from the example 13, m_tet and m_poly.

I observed every edge from the m_tet mesh and its corresponding face from the m_poly mesh (shared face between edge points). I realized that:

• the face is not perpendicular to the edge and
• the face is not always planar

Is that a desired behaviour? Am I doing it wrong?

This is the code snippet I inserted into the example code for the purpose of this analysis:

    for (int i=0; i<m_tet.num_edges(); ++i)
      {
    std::vector<uint> vert_ids = m_tet.edge_vert_ids(i);
    uint fid = m_poly.poly_shared_face(vert_ids.at(0), vert_ids.at(1));

    std::vector<vec3d> edgev = m_tet.edge_verts(i);
    std::vector<vec3d> facev = m_poly.face_verts(fid);

    std::ofstream edgeout ("/tmp/cino/edge-" + std::to_string(i) + ".txt");
    for (auto v: edgev)
        edgeout << v.x() << " " << v.y() << " " << v.z() << "\n";
    std::ofstream faceout ("/tmp/cino/face-" + std::to_string(i) + ".txt");
    for (auto v: facev)
        faceout << v.x() << " " << v.y() << " " << v.z() << "\n";
      }

And some screenshots from Paraview:

[mlivesu]

Hi @Mengesh, thanks for appreciating the project!

In the example you mention the dual mesh also contains clipped (boundary) cells. Are the problems you are mentioning involving boundary faces of those clipped elements?

This is the only issue that comes to my mind, because that primal mesh is Delaunay, hence its dual should be a plain Voronoi diagram with all planar faces

[Mengesh]

Unfortunately no.

[mlivesu]

Ok. Then it may be a bug. I will investigate it and let you know. I'll keep the issue open as a personal reminder :)

[mlivesu]

I've had some time to investigate this issue. It seems that the reason why dual faces are not planar is because in the dualization code I take the centroid of each primal element, whereas to reproduce the Voronoi diagram starting from a Delaunay mesh I should take the center of the sphere passing through the four vertices of the tetrahedron (i.e. its circumcenter).

Note that dual_mesh is more general, and it's not meant to be a tool to pass from Voronoi to Dealunay, or vice-versa. As such, it does not make any assumption on the primal mesh, which is not required to be a Delaunay mesh. Not even to be a tetmesh, actually.

For your specific case, to obtain the desired result I suggest you first make a dual mesh as usual, and then update the vertices of the each dual vertex with the circumcenter of the primal element they correspond to. The first NP vertices in the dual mesh are in 1:1 relation with the polyhedra in the prima mesh. I have created a new function, called tetrahedron_circumcenter, which you can find inside <cinolib/geometry/tetrahedron_utils.h>. You can use it to compute the circumcenter.

[Mengesh]

Thank you very much! I did it like this

for (int i=0; i<m_tet.num_polys(); ++i)
  {
std::vector<vec3d> tetraverts = m_tet.poly_verts(i);
m_poly.vert(i) = tetrahedron_circumcenter(tetraverts.at(0),
                      tetraverts.at(1),
                      tetraverts.at(2),
                      tetraverts.at(3));
  };

and it seems to work for non-boundary edges. However, the results for boundary edges are wild :) The faces are not only non-planar, but sometimes also scattered, as you can see on the image below. The arrows are there just so I can check if face normal and edge are parallel.

I guess this happens because of the clipping mechanism. Do you have any suggestion how to fix this?

[mlivesu]

Correct, this is due to clipped cells, for which there are no guarantees that the dual faces will be planar. For the scattering, I am not sure I understand what you are doing/mean, but if you are restricting your processing to non clipped cells only, there may certainly be gaps (e.g. when two opposite clipped cells are adjacent to one another).

To my knowledge and understanding, the only meshing algorithm that will give you what you are looking for is VoroCrust. They have had a simple yet great idea: rather than clipping cells, they carefully position Voronoi seeds so that the cells will conform to the target geometry. In this case there will be no clipped cells at all. This ensures that the mesh natively meets the orthogonality condition, which is important in CFD. Is this your field?

[Mengesh]

Yes, partially.

I actually found about VoroCrust before cinolib, but I was unable to get in touch with Mohamed Ebeida. I see their website went through some serious update, so I'll try to contact them again. Thank you for the reminder :slightly_smiling_face: