search

CGAL / cgal

The public CGAL repository, see the README below
https://github.com/CGAL/cgal#readme
Other
4.68k stars 1.35k forks source link

Segmentation fault when using make_mesh_3() with edges from corefine_and_compute_intersection() #7356

Closed f-p-b closed 1 year ago

f-p-b commented 1 year ago

Issue Details

I am trying to generate a volume mesh from a surface obtained by performing a corefine_and_compute_intersection(). The mesh should preserve the relevant edges of the original geometry as well as those resulting from the intersection. While the code works if I don’t pass the edges of the original inputs, as soon as I pass the CGAL::parameters::edge_is_constrained_map() detected with detect_sharp_edges() for one of the original geometries the code either gets stuck or I get a segmentation fault (if there are edges from the corresponding map to preserve in the final geometry) . The code bellow and linked files (https://we.tl/t-3B38kjyeYw) should reproduce the issue.

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/IO/STL.h>
#include <CGAL/boost/graph/copy_face_graph.h>

#include <CGAL/Polygon_mesh_processing/detect_features.h>

#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>

#include <CGAL/Surface_mesh/Surface_mesh.h>

// 3D volume mesh
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>

// 3D volume mesh (polyhedral domain with features)
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;

typedef CGAL::Surface_mesh<Point_3> SurfaceMesh;
typedef boost::graph_traits<SurfaceMesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;

// 3D volume meshing (from polyhedron domain with features)
typedef CGAL::Polyhedral_mesh_domain_with_features_3<Kernel> F_Polyhedron_domain;

typedef CGAL::Mesh_polyhedron_3<Kernel>::type F_Polyhedron;

typedef CGAL::Mesh_triangulation_3<F_Polyhedron_domain, CGAL::Default, CGAL::Parallel_if_available_tag>::type F_VTr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<F_VTr, F_Polyhedron_domain::Corner_index, F_Polyhedron_domain::Curve_index> F_C3t3;

typedef CGAL::Mesh_criteria_3<F_VTr> F_Mesh_criteria;

bool polehedral2volume (SurfaceMesh surfaceMesh, 
                        SurfaceMesh::Property_map<edge_descriptor, bool> ecmap, 
                        F_C3t3 &c3t3) {

    // Convert to polygon mesh
    F_Polyhedron polygonSurface;
    CGAL::copy_face_graph(surfaceMesh, polygonSurface);

    // Create domain
    F_Polyhedron_domain domain(polygonSurface);

    // Specify edges to protect
    std::array<vertex_descriptor, 2> v_endpoints;
    std::vector<std::vector<Point_3>> protected_features;

       // Create polylines
    int counter_B = 0;
    for (edge_descriptor e : edges(surfaceMesh)) {
        if (get(ecmap, e)) {
            std::vector<Point_3> polyline;

            // Insert first edge
            polyline.push_back( surfaceMesh.point( surfaceMesh.vertex(e, 0) ) );
            polyline.push_back( surfaceMesh.point( surfaceMesh.vertex(e, 1) ) );
            put(ecmap, e, false);

            // 
            v_endpoints[0] = surfaceMesh.vertex(e, 0);
            v_endpoints[1] = surfaceMesh.vertex(e, 1);

            int counter_b = 2;
            while (1) {
                counter_b += 1;
                // Grow polyline by finding connected edges
                int counter = 0;
                for (edge_descriptor e_new : edges(surfaceMesh)) {
                    counter += 1;
                    if (get(ecmap, e_new)) {

                        // 
                        if (surfaceMesh.vertex(e_new, 0) == v_endpoints[0]) {
                            polyline.insert(polyline.begin(), surfaceMesh.point(surfaceMesh.vertex(e_new, 1)) );
                            v_endpoints[0] = surfaceMesh.vertex(e_new, 1);
                            put(ecmap, e_new, false);

                            break;
                        } else if (surfaceMesh.vertex(e_new, 1) == v_endpoints[0]) {
                            polyline.insert(polyline.begin(), surfaceMesh.point(surfaceMesh.vertex(e_new, 0)) );
                            v_endpoints[0] = surfaceMesh.vertex(e_new, 0);
                            put(ecmap, e_new, false);

                            break;
                        } else if (surfaceMesh.vertex(e_new, 0) == v_endpoints[1]) {
                            polyline.push_back( surfaceMesh.point(surfaceMesh.vertex(e_new, 1)) );
                            v_endpoints[1] = surfaceMesh.vertex(e_new, 1);
                            put(ecmap, e_new, false);

                            break;
                        } else if (surfaceMesh.vertex(e_new, 1) == v_endpoints[1]) {
                            polyline.push_back( surfaceMesh.point(surfaceMesh.vertex(e_new, 0)) );
                            v_endpoints[1] = surfaceMesh.vertex(e_new, 0);
                            put(ecmap, e_new, false);

                            break;
                        }
                    }
                }

                // 
                if ( v_endpoints[0] == v_endpoints[1] ) {
                    break;
                } else if (counter >= surfaceMesh.number_of_vertices() || counter_b >= 100) {
                    break;
                }
            }
            protected_features.push_back(polyline);
        }
    }

    int count = 0;
    for (const std::vector<Point_3> polyline : protected_features) { //
        count += 1;
        std::cout << count << ": " << polyline.size() << " >> ";
        for ( auto ee : polyline)
            std::cout << ee << " " << std::flush;
        std::cout << std::endl;
    }

    domain.add_features(protected_features.begin(), protected_features.end());

    // Mesh criteria
    F_Mesh_criteria criteria(CGAL::parameters::edge_size = 0.1,
                             CGAL::parameters::facet_angle = 30.0,
                             CGAL::parameters::facet_size = 0.1,
                             CGAL::parameters::facet_distance = 0.015,
                             CGAL::parameters::cell_radius_edge_ratio = 3.0,
                             CGAL::parameters::cell_size = 0.1);

    c3t3 = CGAL::make_mesh_3<F_C3t3>(domain, criteria);

    return EXIT_SUCCESS;
}

int main(int argc, char** argv) {

    // 
    std::vector<Point_3> points_A, points_B;
    std::vector<std::vector<std::size_t>> polygons_A, polygons_B;
    SurfaceMesh cube_A, cube_B, cube_OUT;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_A = cube_A.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_B = cube_B.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_OUT = cube_OUT.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;

    std::cout << "  Opening " << "cube.stl" << std::endl;
    std::ifstream stl_file_A("cube.stl", std::ios::binary);
    CGAL::IO::read_STL("cube.stl", points_A, polygons_A);
    stl_file_A.close();

    std::cout << "  Opening " << "cube_b.stl" << std::endl;
    std::ifstream stl_file_B("cube_b.stl", std::ios::binary);
    CGAL::IO::read_STL("cube_b.stl", points_B, polygons_B);
    stl_file_B.close();

    //
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_A, polygons_A, cube_A);
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_B, polygons_B, cube_B);

    // Detect sharp edges
    CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_A, 60, ecmap_A);
    CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_B, 60, ecmap_B);

    if(CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(cube_B, cube_A, cube_OUT, CGAL::parameters::edge_is_constrained_map(ecmap_B), CGAL::parameters::default_values(), CGAL::parameters::edge_is_constrained_map(ecmap_OUT)) == false)
        std::cout << "\n  WARNING: corefine_and_compute_intersection output may have only been corefined!" << std::endl;

    std::ofstream("cube_OUT.off") << std::setprecision(17) << cube_OUT;

    F_C3t3 c3t3;
    polehedral2volume(cube_OUT, ecmap_OUT, c3t3);

    // Write output to .mesh file
    std::ofstream medit_file("cube_OUT.mesh");
    c3t3.output_to_medit(medit_file);
    medit_file.close();
}

Environment

lrineau commented 1 year ago

I can reproduce. There is an issue with the polylines you pass to domain.add_features(...): they do not comply with the expectations. I will try to find a fix for you.

Screenshot_20230330_105312 protected_features.polylines.txt

lrineau commented 1 year ago

Here is a new version of the code. It replaces one of your piece of code with a use of CGAL::split_graph_into_polylines.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/IO/STL.h>
#include <CGAL/boost/graph/copy_face_graph.h>

#include <CGAL/Polygon_mesh_processing/detect_features.h>

#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>

#include <CGAL/Surface_mesh/Surface_mesh.h>

#include <CGAL/boost/graph/split_graph_into_polylines.h>

// 3D volume mesh
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>

// 3D volume mesh (polyhedral domain with features)
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;

typedef CGAL::Surface_mesh<Point_3> SurfaceMesh;
typedef boost::graph_traits<SurfaceMesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;

// 3D volume meshing (from polyhedron domain with features)
typedef CGAL::Polyhedral_mesh_domain_with_features_3<Kernel, SurfaceMesh> F_Polyhedron_domain;

typedef CGAL::Mesh_triangulation_3<F_Polyhedron_domain, CGAL::Default, CGAL::Parallel_if_available_tag>::type F_VTr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<F_VTr, F_Polyhedron_domain::Corner_index, F_Polyhedron_domain::Curve_index> F_C3t3;

typedef CGAL::Mesh_criteria_3<F_VTr> F_Mesh_criteria;

struct Polyline_visitor
{
  typedef std::vector<Point_3> Polyline;

  std::vector<Polyline>& polylines;
  const SurfaceMesh& mesh;

  void start_new_polyline()
  {
    polylines.emplace_back();
  }

  void add_node(vertex_descriptor vd)
  {
    polylines.back().push_back(mesh.point(vd));
  }

  void end_polyline(){}
};

bool polehedral2volume (const SurfaceMesh& surfaceMesh, F_C3t3 &c3t3) {

    // Create domain
    F_Polyhedron_domain domain(surfaceMesh);

    // Specify edges to protect
    std::vector<std::vector<Point_3>> protected_features;
    Polyline_visitor visitor{protected_features, surfaceMesh};
    CGAL::split_graph_into_polylines(surfaceMesh, visitor);

    int count = 0;
    std::ofstream polylines_file("protected_features.polylines.txt");
    polylines_file.precision(17);
    for (const std::vector<Point_3>& polyline : protected_features) { //
            polylines_file << polyline.size() << "  ";
        count += 1;
        std::cout << count << ": " << polyline.size() << " >> ";
        for ( auto ee : polyline) {
            std::cout << ee << " " << std::flush;
            polylines_file << ee << " ";
        }
        polylines_file << std::endl;
        std::cout << std::endl;
    }

    // domain.add_features(protected_features.begin(), protected_features.end());
    domain.detect_features(10);

    // Mesh criteria
    F_Mesh_criteria criteria(CGAL::parameters::edge_size = 0.1,
                             CGAL::parameters::facet_angle = 30.0,
                             CGAL::parameters::facet_size = 0.1,
                             CGAL::parameters::facet_distance = 0.015,
                             CGAL::parameters::cell_radius_edge_ratio = 3.0,
                             CGAL::parameters::cell_size = 0.1);

    c3t3 = CGAL::make_mesh_3<F_C3t3>(domain, criteria);

    return EXIT_SUCCESS;
}

int main(int, char**) {

    // 
    std::vector<Point_3> points_A, points_B;
    std::vector<std::vector<std::size_t>> polygons_A, polygons_B;
    SurfaceMesh cube_A, cube_B, cube_OUT;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_A = cube_A.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_B = cube_B.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
    SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_OUT = cube_OUT.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;

    std::cout << "  Opening " << "cube.stl" << std::endl;
    std::ifstream stl_file_A("cube.stl", std::ios::binary);
    CGAL::IO::read_STL("cube.stl", points_A, polygons_A);
    stl_file_A.close();

    std::cout << "  Opening " << "cube_b.stl" << std::endl;
    std::ifstream stl_file_B("cube_b.stl", std::ios::binary);
    CGAL::IO::read_STL("cube_b.stl", points_B, polygons_B);
    stl_file_B.close();

    //
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_A, polygons_A, cube_A);
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_B, polygons_B, cube_B);

    // Detect sharp edges
    CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_A, 60, ecmap_A);
    CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_B, 60, ecmap_B);

    if(CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(cube_B, cube_A, cube_OUT, CGAL::parameters::edge_is_constrained_map(ecmap_B), CGAL::parameters::default_values(), CGAL::parameters::edge_is_constrained_map(ecmap_OUT)) == false)
        std::cout << "\n  WARNING: corefine_and_compute_intersection output may have only been corefined!" << std::endl;

    std::ofstream("cube_OUT.off") << std::setprecision(17) << cube_OUT;

    F_C3t3 c3t3;
    polehedral2volume(cube_OUT, c3t3);

    // Write output to .mesh file
    std::ofstream medit_file("cube_OUT.mesh");
    c3t3.output_to_medit(medit_file);
    medit_file.close();
}

And then the result is fine.

f-p-b commented 1 year ago

I have tested the new code and the issue seems to still be there. If I run it the code as posted it works since the line "domain.add_features(protected_features.begin(), protected_features.end())" is commented out and domain.detect_features(10) is used instead. But if I uncomment add_features() and comment the detect_features() it gets gets stuck like before (waited for 3 minutes and then killed it)

lrineau commented 1 year ago

I am sorry. You are right. Here is the second fixed version.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/IO/STL.h>
#include <CGAL/boost/graph/copy_face_graph.h>

#include <CGAL/Polygon_mesh_processing/detect_features.h>

#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/repair.h>

#include <CGAL/Surface_mesh/Surface_mesh.h>

#include <CGAL/boost/graph/split_graph_into_polylines.h>

// 3D volume mesh
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/make_mesh_3.h>

// 3D volume mesh (polyhedral domain with features)
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;

typedef CGAL::Surface_mesh<Point_3> SurfaceMesh;
typedef boost::graph_traits<SurfaceMesh>::edge_descriptor edge_descriptor;
typedef boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;

template <typename ECMap>
struct Is_constrained {
  ECMap ecmap;

 template <typename Edge>
  bool operator()(const Edge& e) const {
    return get(ecmap, e);
  }
};

template <typename VPMap>
struct Is_feature {
  VPMap vpmap;

 template <typename Vertex>
  bool operator()(const Vertex& v) const {
    return get(vpmap, v);
  }
};

// 3D volume meshing (from polyhedron domain with features)
typedef CGAL::Polyhedral_mesh_domain_with_features_3<Kernel, SurfaceMesh> F_Polyhedron_domain;

typedef CGAL::Mesh_triangulation_3<F_Polyhedron_domain, CGAL::Default, CGAL::Parallel_if_available_tag>::type F_VTr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<F_VTr, F_Polyhedron_domain::Corner_index, F_Polyhedron_domain::Curve_index> F_C3t3;

typedef CGAL::Mesh_criteria_3<F_VTr> F_Mesh_criteria;

struct Polyline_visitor
{
  typedef std::vector<Point_3> Polyline;

  std::vector<Polyline>& polylines;
  const SurfaceMesh& mesh;

  void start_new_polyline()
  {
    polylines.emplace_back();
  }

  void add_node(vertex_descriptor vd)
  {
    polylines.back().push_back(mesh.point(vd));
  }

  void end_polyline(){}
};

bool polehedral2volume (SurfaceMesh& surfaceMesh,
                        SurfaceMesh::Property_map<edge_descriptor, bool> ecmap,
                        F_C3t3 &c3t3) {

        // Create domain
        F_Polyhedron_domain domain(surfaceMesh);

    auto v_pmap = surfaceMesh.add_property_map<vertex_descriptor, bool>("v:on_feature_curve", false).first;
    for(auto e : edges(surfaceMesh)) {
        auto v1 = source(e, surfaceMesh);
        auto v2 = target(e, surfaceMesh);
        if (get(ecmap, e)) {
            put(v_pmap, v1, true);
            put(v_pmap, v2, true);
        }
    }

        typedef SurfaceMesh::Property_map<edge_descriptor, bool> ECMap;
    typedef SurfaceMesh::Property_map<vertex_descriptor, bool> VPMap;
        typedef boost::filtered_graph<SurfaceMesh,
                                  Is_constrained<ECMap>,
                      Is_feature<VPMap>> FG;

        // Specify edges to protect
        std::vector<std::vector<Point_3>> protected_features;
        Polyline_visitor visitor{protected_features, surfaceMesh};
        Is_constrained<ECMap> is_constrained{ecmap};
    Is_feature<VPMap> is_feature{v_pmap};
        FG constrained_edges{surfaceMesh, is_constrained, is_feature};
        CGAL::split_graph_into_polylines(constrained_edges, visitor);

        int count = 0;
        std::ofstream polylines_file("protected_features.polylines.txt");
        polylines_file.precision(17);
        for (const std::vector<Point_3>& polyline : protected_features) { //
                polylines_file << polyline.size() << "  ";
                count += 1;
                std::cout << count << ": " << polyline.size() << " >> ";
                for ( auto ee : polyline) {
                        std::cout << ee << " " << std::flush;
                        polylines_file << ee << " ";
                }
                polylines_file << std::endl;
                std::cout << std::endl;
        }
        polylines_file.close();

        domain.add_features(protected_features.begin(), protected_features.end());

        // Mesh criteria
        F_Mesh_criteria criteria(CGAL::parameters::edge_size = 0.1,
                                 CGAL::parameters::facet_angle = 30.0,
                                 CGAL::parameters::facet_size = 0.1,
                                 CGAL::parameters::facet_distance = 0.015,
                                 CGAL::parameters::cell_radius_edge_ratio = 3.0,
                                 CGAL::parameters::cell_size = 0.1);

        c3t3 = CGAL::make_mesh_3<F_C3t3>(domain, criteria);

        return EXIT_SUCCESS;
}

int main(int, char**) {

        // 
        std::vector<Point_3> points_A, points_B;
        std::vector<std::vector<std::size_t>> polygons_A, polygons_B;
        SurfaceMesh cube_A, cube_B, cube_OUT;
        SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_A = cube_A.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
        SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_B = cube_B.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;
        SurfaceMesh::Property_map<edge_descriptor, bool> ecmap_OUT = cube_OUT.add_property_map<edge_descriptor, bool>("e:is_constrained", false).first;

        std::cout << "  Opening " << "cube.stl" << std::endl;
        std::ifstream stl_file_A("cube.stl", std::ios::binary);
        CGAL::IO::read_STL("cube.stl", points_A, polygons_A);
        stl_file_A.close();

        std::cout << "  Opening " << "cube_b.stl" << std::endl;
        std::ifstream stl_file_B("cube_b.stl", std::ios::binary);
        CGAL::IO::read_STL("cube_b.stl", points_B, polygons_B);
        stl_file_B.close();

        //
        CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_A, polygons_A, cube_A);
        CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points_B, polygons_B, cube_B);

        // Detect sharp edges
        CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_A, 60, ecmap_A);
        CGAL::Polygon_mesh_processing::detect_sharp_edges(cube_B, 60, ecmap_B);

        if(CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(cube_B, cube_A, cube_OUT, CGAL::parameters::edge_is_constrained_map(ecmap_B), CGAL::parameters::default_values(), CGAL::parameters::edge_is_constrained_map(ecmap_OUT)) == false)
                std::cout << "\n  WARNING: corefine_and_compute_intersection output may have only been corefined!" << std::endl;
        CGAL::Polygon_mesh_processing::remove_isolated_vertices(cube_OUT);
        std::ofstream("cube_OUT.off") << std::setprecision(17) << cube_OUT;

        F_C3t3 c3t3;
        polehedral2volume(cube_OUT, ecmap_OUT, c3t3);

        // Write output to .mesh file
        std::ofstream medit_file("cube_OUT.mesh");
        c3t3.output_to_medit(medit_file);
        medit_file.close();
}
f-p-b commented 1 year ago

Thanks, it working very nicely now!

I have one follow up question. Even though its doing exactly what it is supposed to, sometimes the result is not ideal. In corners were the cut from compute intersection results in a small edge the mesher tends to use more elements than ideal, in the example bellow it is visible how it subdivided the 3 small edges in multiple edges while not subdividing them would have been preferred. Is there any way to ‘encourage’ the mesher a little towards the later?

image image

lrineau commented 1 year ago

I see. You have a small triangular face in the input, with small edges, and it is meshed with even smaller edges. Unfortunately, that behavior is mandatory to ensure that the triangular face is well represented in the weighted Delaunay mesh.

As you are probably not interested by the Delaunay property, you may want to try remesh the result using Tetrahedral Remeshing. Please share the two input data meshes corresponding to that picture, and we will experiment.

Cc: @janetournois (author of Tetrahedral Remeshing in CGAL).

f-p-b commented 1 year ago

I have shared them though the following link. However, to avoid wasting your time for no reason a bit more context since I think Tetrahedral Remeshing may not offer an ideal solution.

Although for the meshes I just shared it’s not necessary since they are a trivial example using a cube, in general I pass on sizing fields (for cell, edge and facet size) when calling make_mesh_3(). These sizing fields depend on local thickness. I do this because generally the models have a porous structure of varying wall thickness (one region may have a structure that is 5 to 10 times finer than that of another region for example). Now, if I recall correctly from when I looked into using Tetrahedral Remeshing a while back, it is limited to uniform isotropic meshes which would mean the thinnest sections would dictate the size of the whole mesh. This would lead to overly refined meshes for large portions of the model with its respective consequences when using the mesh in simulations than can already be quite time consuming. Compared to this, the current small regions with the overly refined mesh would be preferable, although not ideal.

Under this conditions, can Tetrahedral Remeshing still offer some fix or would a different strategy be requried?

janetournois commented 1 year ago

Hello, you are right about tetrahedral remeshing. It is only working for uniform remeshing, though we are planning to add a sizing field parameter in the future. The best solution for now would be to detect these small undesired features and fix them a priori with a "real" corner.