Trying to implement a variable size criteria but the boundary mesh size seems too big

[citystrawman]

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

Issue Details

I've a tunnel which is in the middle of a big rectangle box (the width and height of the box is about several times the tunnel width and height ), and My final goal is to geneate mesh for this regtangle box which has a smaller mesh size near around the tunnel and bigger mesh size near the box boudnary. To simplify this model , I firstly made a rectangle box without tunnel, and I would like to generate 3D mesh within this box. I studied the Mesh_3/mesh_implicit_sphere_variable_size.cpp example, but in my model, the the mesh size should change within box, thus I proposed a Box_sizing_field as follows:

// 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_ratio;
    FT outer_ratio;

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

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

    void SetRatio(double inner, double outer) {
        if ((inner >= 1 && inner < 0) || (outer > 1 && outer <= 0) || inner >= outer) {
            throw std::exception("invalid ratio!");
        }
        else {
            this->inner_ratio = inner;
            this->outer_ratio = outer;
        }
    }

    void SetMeshSize(double min, double 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_ratio;
        FT y_outer_dist = y_size * 0.5 * outer_ratio;
        FT z_outer_dist = z_size * 0.5 * outer_ratio;
        FT x_inner_dist = x_size * 0.5 * inner_ratio;
        FT y_inner_dist = y_size * 0.5 * inner_ratio;
        FT z_inner_dist = z_size * 0.5 * inner_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) 
        {
            double 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));
            double 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;
        }
    }
};

Then I tested this code using my model. The result works almost OK at the inner core, but the boundary meshes looks too big and the shape is stretched: I doubt I made some mistake but could not solve it myself.

The boundingbox-large.off file is attached at the end of the post.

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_ratio;
    FT outer_ratio;

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

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

    void SetRatio(double inner, double outer) {
        if ((inner >= 1 && inner < 0) || (outer > 1 && outer <= 0) || inner >= outer) {
            throw std::exception("invalid ratio!");
        }
        else {
            this->inner_ratio = inner;
            this->outer_ratio = outer;
        }
    }

    void SetMeshSize(double min, double 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_ratio;
        FT y_outer_dist = y_size * 0.5 * outer_ratio;
        FT z_outer_dist = z_size * 0.5 * outer_ratio;
        FT x_inner_dist = x_size * 0.5 * inner_ratio;
        FT y_inner_dist = y_size * 0.5 * inner_ratio;
        FT z_inner_dist = z_size * 0.5 * inner_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) 
        {
            double 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));
            double 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;
        }
    }
};

// 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 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tunnel/boundingbox-large.off");
    std::ifstream input(fname);
    Polyhedron smbounding;
    input >> smbounding;
    if (input.fail()) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }
    /// [Domain creation] (Warning: Sphere_3 constructor uses squared radius !)
    Mesh_domain domain(smbounding);
    domain.detect_features();
    /// [Domain creation]

    // 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.SetRatio(0.3, 0.8);
    double min_zize = 5; 
    double max_size = 15;
    size_1.SetMeshSize(min_zize, max_size);

    Mesh_criteria criteria(facet_angle = 25, facet_size = max_size, facet_distance = max_size*0.33,
        cell_radius_edge_ratio = 1.2, cell_size = size_1);

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

    // Output
    std::ofstream medit_file("mesh_tunnel_variable_size.mesh");
    c3t3.output_to_medit(medit_file);

    return 0;
}

Environment

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

[janetournois]

Hello this sizing artefact is due to the feature edges size which is too large in your output. To have a smoother sizing transition between sharp features, surfaces and volume, you should add the edge_size parameter to the meshing criteria. For example :

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

[citystrawman]

Hello this sizing artefact is due to the feature edges size which is too large in your output. To have a smoother sizing transition between sharp features, surfaces and volume, you should add the edge_size parameter to the meshing criteria. For example :

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

This works. Thank you so much.