Advice on data structure for custom mesh optimization algorithm

[larsmunaf]

Dear all,

I need advice on deciding which CGAL data structure I should go with for my mesh optimization algorithm.

Problem

I have implemented a tetrahedral mesh optimization algorithm that is currently built upon Triangulation_3/Delaunay_triangulation_3:

using Vb = CGAL::Triangulation_vertex_base_with_info_3<int, K>;
using Cb = CGAL::Triangulation_cell_base_with_info_3<int, K>;
using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;
using Delaunay_triangulation_3 = CGAL::Delaunay_triangulation_3<K, Tds>;
using Triangulation_3 = CGAL::Triangulation_3<K, Tds>;

It does not add or remove vertices and runs well on random (Delaunay) triangulations of unit spheres. It needs an index for each vertex and cell in order to run efficiently. Now, I have some additional requirements that appear difficult to realize with Triangulation_3/Delaunay_triangulation_3:

• I want to be able to read tetmeshes from files, e.g. MEDIT. I did not see how to do this with Triangulation_3. According to the documentation, this is only directly possible with Mesh_complex_3_in_triangulation_3 or with a Triangulation_3 that consists "at least" of a Simplical_mesh_vertex_base_3 and a Simplical_mesh_cell_base_3. Interestingly, what somehow works is creating vertex and cell handles manually in the underlying tds() of Triangulation_3, but CGAL adds extra finite vertices and finite tets, which is not what I want. And adding points manually to Delaunay_triangulation_3 produces the convex hull of the points, which is obviously not what I want either.
• The optimization algorithm should be compatible to the 3D mesh generation package in some way, i.e., to Mesh_complex_3_in_triangulation_3. I would like to be able to generate a C3t3 using the refine() function and then use this C3t3 as a starting point for both the built-in optimization methods (ODT smoothing etc.) as well as for my optimization algorithm.
• I still want to be able to generate random Delaunay triangulations that I can optimize with my algorithm. I think there are CGAL::tetrahedron_soup_to_triangulation_3 and c3t3.set_triangulation() exactly for this use case, right?

Question

Now, with which CGAL data structure am I able to meet the above requirements?

Ideas

• I could use Mesh_triangulation_3 for my optimization algorithm, but add and extra info() field by somehow combining Triangulation_vertex_base_with_info_3 and Mesh_vertex_base_3(and analogously for the cell base)?
• I could use Mesh_triangulation_3 for my optimization algorithm as it is. Then I need to store cell and vertex indices in a separate structure, i.e., an std::map or std::vector.
• I could implement my algorithm solely based on CGAL::Triangulation_data_structure_3. But this is probably not feasible, as I would need to rewrite the optimization algorithm completely.
• I've also read about Tetrahedral_remeshing::Remeshing_triangulation_3. But I am not sure if using it is sensible in my case, as I do not intend to add or remove vertices.

Any help is greatly appreciated, thank you!

[larsmunaf]

I now chose to use Mesh_triangulation_3 with overridden versions of Mesh_vertex_base_3 and Compact_mesh_cell_base_3 that have an info() method:

using K = CGAL::Epick;
using Mesh_domain_3 = CGAL::Labeled_mesh_domain_3<K>;

// Based on https://doc.cgal.org/latest/Triangulation_3/Triangulation_3_2adding_handles_3_8cpp-example.html
template <class GT, class MD, class Vb = CGAL::Mesh_vertex_base_3<GT, MD>>
class Mesh_vertex_base_with_info_3 : public Vb {
   public:
    using Point = typename Vb::Point;

    Mesh_vertex_base_with_info_3() : Vb(), idx(-1) {}
    explicit Mesh_vertex_base_with_info_3(const Point& p) : Vb(p), idx(-1) {}

    int& info() { return idx; }
    const int& info() const { return idx; }

    template <class TDS2>
    struct Rebind_TDS {
        using Vb2 = typename Vb::template Rebind_TDS<TDS2>::Other;
        using Other = Mesh_vertex_base_with_info_3<GT, MD, Vb2>;
    };

   private:
    int idx;
};

// Based on https://doc.cgal.org/latest/Triangulation_3/Triangulation_3_2adding_handles_3_8cpp-example.html
template <class GT, class MD, class Cell_base = CGAL::Compact_mesh_cell_base_3<GT, MD>>
class Compact_mesh_cell_base_with_info_3 : public Cell_base {
   public:
    Compact_mesh_cell_base_with_info_3() : Cell_base(), idx(-1) {}

    int& info() { return idx; }
    const int& info() const { return idx; }

    template <class TDS2>
    struct Rebind_TDS {
        using Base2 = typename Cell_base::template Rebind_TDS<TDS2>::Other;
        using Other = Compact_mesh_cell_base_with_info_3<GT, MD, Base2>;
    };

   private:
    int idx;
};

using Vb = Mesh_vertex_base_with_info_3<K, Mesh_domain_3>;
using Cb = Compact_mesh_cell_base_with_info_3<K, Mesh_domain_3>;
using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;

using Mesh_triangulation_3 = CGAL::Mesh_triangulation_3<Mesh_domain_3, K, CGAL::Sequential_tag, Tds>::type;

using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Mesh_triangulation_3>;

However, when trying to create a C3t3, I get the error

In template: no type named 'Vertex_handle' in 'CGAL::Mesh_triangulation_3<CGAL::Labeled_mesh_domain_3<CGAL::Epick>, CGAL::Epick, CGAL::Sequential_tag, CGAL::Triangulation_data_structure_3<Mesh_vertex_base_with_info_3<CGAL::Epick, CGAL::Labeled_mesh_domain_3<CGAL::Epick>>, Compact_mesh_cell_base_with_info_3<CGAL::Epick, CGAL::Labeled_mesh_domain_3<CGAL::Epick>>>>'

[larsmunaf]

I could solve it by changing

using Mesh_triangulation_3 = CGAL::Mesh_triangulation_3<Mesh_domain_3, K, CGAL::Sequential_tag, Tds>::type;

to

using Mesh_triangulation_3 = CGAL::Mesh_triangulation_3<Mesh_domain_3, K, CGAL::Sequential_tag, Vb, Cb>::type;

But I am still interested if Mesh_triangulation_3 is the best choice for this or if I should go with some other data structure?

[janetournois]

Hello @larsmunaf ,

thank you for your interest in the CGAL data structures.

I think the best available data structure is the Mesh_complex_3_in_triangulation_3, with a triangulation as close as possible to the Remeshing_triangulation_3, which uses "only" the SimplicialMeshVertexBase_3 (and cell). As you already found out, you can use the function c3t3.set_triangulation() to set the triangulation in the C3T3.

If you use a simpler triangulation, with only the index information added to the vertex/cell bases, you may miss some of the additional info that is stored in the Medit file that you are using as input, including the Subdomain_index (index of subdomain in cells), or Surface_patch_index (index of surface patch in surface facets), etc.

You can use the function CGAL::IO::read_MEDIT() to automatically build a CGAL triangulation from the file, or one of the CGAL::tetrahedron_soup_to_triangulation_3() functions to build your triangulation.

I want to be able to read tetmeshes from files, e.g. MEDIT. I did not see how to do this with Triangulation_3. According to the documentation, this is only directly possible with Mesh_complex_3_in_triangulation_3 or with a Triangulation_3 that consists "at least" of a Simplical_mesh_vertex_base_3 and a Simplical_mesh_cell_base_3. Interestingly, what somehow works is creating vertex and cell handles manually in the underlying tds() of Triangulation_3, but CGAL adds extra finite vertices and finite tets, which is not what I want. And adding points manually to Delaunay_triangulation_3 produces the convex hull of the points, which is obviously not what I want either.

About that : there is no reason why CGAL would add extra finite vertices. It only adds one infinite vertex, to have a fully valid connectivity. Each and every facet, including the ones of the convex hull of a triangulation, has two valid neighboring cells (two finite, or one finite and one infinite, or two infinite). It adds finite cells to cover the convex hull of the triangulation vertices.

It is important to mention that a CGAL triangulation is meant to be a triangulation of the convex hull of its vertices. In the case where the input triangulation (read from Medit file for example) does not cover the entire convex hull, the Triangulation_data_structure will remain valid, but the location functions called outside the outer hull (not convex) may return wrong results. In your case, because you have no vertex insertions and vertex removals, and your optimization operations happen inside the outer hull, things should go smoothly, but be careful :-)

I hope it helps!

[larsmunaf]

Hello @janetournois,

thank you for the detailed insights! Then I think I am on the right way with Mesh_complex_3_in_triangulation_3 and Mesh_triangulation_3 using subclasses of Mesh_vertex_base_3 and Compact_mesh_cell_base_3 with an additional info() field. This works perfectly fine right now.

Two follow-up questions:

• Assuming I choose "only" the Remeshing_triangulation_3 with a SimplicialMeshVertexBase_3. According to the documentation, this will not work with CGAL::IO::read_MEDIT(), as the vertex and cell base class must be models of the concepts MeshVertexBase_3 and MeshCellBase_3, right? But when using the same Remeshing_triangulation_3 in a Mesh_complex_3_in_triangulation_3, I should be able to call CGAL::IO::read_MEDIT() on the Mesh_complex_3_in_triangulation_3, right?
• CGAL::IO::read_MEDIT() expects Surface_patch_index to be an int. create_implicit_mesh_domain() expects it to be an std::pair by default. Usually, I would need to create two variants of Mesh_triangulation_3, one with std::pair and the other with int as the Surface_patch_index. In order to avoid this, I pass a surface index wrapper to create_implicit_mesh_domain(), that converts std::pair to int. As I do not deal with multiple domains, I do not need the Subdomain_index and Surface_patch_index at all. Using this std::pair to int wrapper should thus be fine, right?

Regards!

[janetournois]

Hello @larsmunaf ,

• Assuming I choose "only" the Remeshing_triangulation_3 with a SimplicialMeshVertexBase_3. According to the documentation, this will not work with CGAL::IO::read_MEDIT(), as the vertex and cell base class must be models of the concepts MeshVertexBase_3 and MeshCellBase_3, right? But when using the same Remeshing_triangulation_3 in a Mesh_complex_3_in_triangulation_3, I should be able to call CGAL::IO::read_MEDIT() on the Mesh_complex_3_in_triangulation_3, right?

actually I think that is a documentation mistake and the Simplicial... should be enough to use CGAL::IO::read_MEDIT()

CGAL::IO::read_MEDIT() expects Surface_patch_index to be an int. create_implicit_mesh_domain() expects it to be an std::pair by default. Usually, I would need to create two variants of Mesh_triangulation_3, one with std::pair and the other with int as the Surface_patch_index. In order to avoid this, I pass a surface index wrapper to create_implicit_mesh_domain(), that converts std::pair to int. As I do not deal with multiple domains, I do not need the Subdomain_index and Surface_patch_index at all. Using this std::pair to int wrapper should thus be fine, right?

yes, that sounds good!