Can CGAL::make_mesh_3 method keep the input surface mesh structure unchange while generating 3D mesh?

[citystrawman]

Please use the following template to help us solving your issue.

This is a continuation for issue #7387 . To sum up, I need to generate 3D mesh of a tunnel (tunnel.off) as well as its bounding box (boundingbox-large.off). The mesh criteria is that, the nearer a point is to the tunnel, the denser that a mesh should be. Therefore the mesh looks like this: and this is what I have accomplished so far.

However, when I checked the 3D mesh within the tunnel model (I import the mesh file as well as tunnel.off into FLAC and it is able to display the mesh within the tunnel), the mesh looks as it is remeshed and does not retain some features. You can see the 3D mesh and the original off surface mesh comparision:

(this is the 3D mesh within tunnel)

(this is the original off surface mesh)

What I want to achieve is that, the generated 3D mesh should as close as the original surface mesh while retaining the mesh criteria. I know it may be hard to achive this as the mesh size is different from the original surface mesh triangle edge length, but is there a best way to achieve this? Actually by using rhino griddle add-in I could get the 3D mesh within the tunnel looks like this: which looks perfect.

Describe your issue. Please be specific (compilation error, runtime error, unexpected behavior, wrong results, etc.).

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/clip.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>

#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/boost/graph/helpers.h>

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedron_3<K> Polyhedron_3;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT(Function)(const Point&);
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;

#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type Tr;

typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

// Sizing field
struct Box_sizing_field {
    typedef ::FT FT;
    typedef Point Point_3;
    typedef Mesh_domain::Index Index;

    Point_3 center;
    FT x_size;
    FT y_size;
    FT z_size;

    FT min_mesh_size;
    FT max_mesh_size;

    FT inner_x_ratio;
    FT outer_x_ratio;

    FT inner_y_ratio;
    FT outer_y_ratio;

    FT inner_z_ratio;
    FT outer_z_ratio;

    void SetCenter(Point_3 point) {
        this->center = point;
    }

    void SetSideSize(FT x, FT y, FT z) {
        this->x_size = x;
        this->y_size = y;
        this->z_size = z;
    }

    void SetXRatio(FT inner, FT outer) {
        if ((inner >= 1 && inner < 0) || (outer > 1 && outer <= 0) || inner >= outer) {
            throw std::exception("invalid ratio!");
        }
        else {
            this->inner_x_ratio = inner;
            this->outer_x_ratio = outer;
        }
    }

    void SetYRatio(FT inner, FT outer) {
        if ((inner >= 1 && inner < 0) || (outer > 1 && outer <= 0) || inner >= outer) {
            throw std::exception("invalid ratio!");
        }
        else {
            this->inner_y_ratio = inner;
            this->outer_y_ratio = outer;
        }
    }

    void SetZRatio(FT inner, FT outer) {
        if ((inner >= 1 && inner < 0) || (outer > 1 && outer <= 0) || inner >= outer) {
            throw std::exception("invalid ratio!");
        }
        else {
            this->inner_z_ratio = inner;
            this->outer_z_ratio = outer;
        }
    }

    void SetMeshSize(FT min, FT max) {
        if (min <= 0 || max <= 0 || min > max) {
            throw std::exception("invalid mesh size!");
        }
        else {
            this->min_mesh_size = min;
            this->max_mesh_size = max;
        }
    }
    //
    FT operator()(const Point_3& p, const int, const Index&) const
    {
        FT dist_x = CGAL::abs(p.x() - center.x());
        FT dist_y = CGAL::abs(p.y() - center.y());
        FT dist_z = CGAL::abs(p.z() - center.z());

        FT dist_to_center;

        FT x_outer_dist = x_size * 0.5 * outer_x_ratio;
        FT y_outer_dist = y_size * 0.5 * outer_y_ratio;
        FT z_outer_dist = z_size * 0.5 * outer_z_ratio;
        FT x_inner_dist = x_size * 0.5 * inner_x_ratio;
        FT y_inner_dist = y_size * 0.5 * inner_y_ratio;
        FT z_inner_dist = z_size * 0.5 * inner_z_ratio;

        if (dist_x <= x_inner_dist && dist_y <= y_inner_dist && dist_z <= z_inner_dist)
        {
            return min_mesh_size;
        }
        else if (dist_x <= x_outer_dist && dist_y <= y_outer_dist && dist_z <= z_outer_dist)
        {
            FT maxxy = CGAL::max((dist_x - x_inner_dist) / (x_outer_dist - x_inner_dist),
                (dist_y - y_inner_dist) / (y_outer_dist - y_inner_dist));
            FT max_dist_ratio = CGAL::max(maxxy, (dist_z - z_inner_dist) / (z_outer_dist - z_inner_dist));
            return min_mesh_size + max_dist_ratio * (max_mesh_size - min_mesh_size);
        }
        else
        {
            return max_mesh_size;
        }
    }
};

FT Box_Function(const Point& p, FT x_size, FT y_size, FT z_size, Point center)
{
    FT x = p.x();
    FT y = p.y();
    FT z = p.z();
    FT x_center = center.x();
    FT y_center = center.y();
    FT z_center = center.z();

    FT xy = CGAL::max(CGAL::abs(x - x_center) / x_size, CGAL::abs(y - y_center) / y_size);
    FT xyz = CGAL::max(xy, CGAL::abs(z - z_center) / z_size);
    return xyz = 1;
}

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
namespace PMP = CGAL::Polygon_mesh_processing;

int main(int argc, char* argv[])
{
    std::cout.precision(17);
    std::cerr.precision(17);

    //Read  file
    const std::string fname_t = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tunnel/tunnel.off");
    const std::string fname_b = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tunnel/boundingbox-large.off");
    std::ifstream input_t(fname_t);
    Polyhedron sm;
    input_t >> sm;
    if (input_t.fail()) {
        std::cerr << "Error: Cannot read file " << fname_t << std::endl;
        return EXIT_FAILURE;
    }
    std::ifstream input_b(fname_b);
    Polyhedron smbounding;
    input_b >> smbounding;
    if (input_b.fail()) {
        std::cerr << "Error: Cannot read file " << fname_b << std::endl;
        return EXIT_FAILURE;
    }

    PMP::clip(sm, smbounding);

    CGAL::Timer t;
    t.start();

    /// [Domain creation] (Warning: Sphere_3 constructor uses squared radius !)
    Mesh_domain domain(sm, smbounding);
    domain.detect_features();

    // Mesh criteria
    Box_sizing_field size_1;
    size_1.SetCenter(Point(82.74825, 321.0794, 2621.121)); //these numbers are corresponding to the boundingbox-large.off model
    size_1.SetSideSize(557.8839, 854.6447, 261.469); //these numbers are corresponding to the boundingbox-large.off model
    size_1.SetXRatio(0.3, 0.8);
    size_1.SetYRatio(0.75, 0.9);
    size_1.SetZRatio(0.3, 0.8);
    FT min_zize = 3;
    FT max_size = 15;
    size_1.SetMeshSize(min_zize, max_size);

    Mesh_criteria criteria(edge_size = max_size, facet_angle = 25, facet_size = max_size, facet_distance = max_size * 0.8,
        cell_radius_edge_ratio = 3, cell_size = size_1);

    // Mesh generation
    C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), exude());

    std::cerr << t.time() << " sec." << std::endl;
    // Output
    std::ofstream medit_file("mesh_tunnel_variable_size.mesh");
    c3t3.output_to_medit(medit_file);

    return 0;
}

Attached is the file in code data.zip

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Win10
• Compiler: visual studio 2019
• Release or debug mode:release
• Specific flags used (if any):
• CGAL version:
• Boost version:
• Other libraries versions if used (Eigen, TBB, etc.):

[citystrawman]

To make it more clear, what I need is to find a way to make_mesh_3 while keep the domain's inner surface's original geometric features. It is not necessary to keep the inner surface vertex remain same, but I need to keep its geometry features as much as possible(such as retaining its surface curvatures, boundary edges, etc.)

[lrineau]

I think you are looking for a mesh based on a 3D constrained Delaunay triangulation (CDT). That component does not exist (yet) in CGAL. Anyway Mesh_3 cannot provide what you ask for. Meshing tool based on 3D CDT exist. The software tetgen is quite popular.

In the future, GeometryFactory plan to have 3D CDT in CGAL, and meshing algorithms based on CDT. But that piece of software is not ready for release for now.