search

CGAL / cgal

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

Trying to make a variable size criteria for generating 3D mesh but failed #7384

Closed citystrawman closed 1 year ago

citystrawman commented 1 year ago

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 5 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:

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_outer_dist || dist_y > y_outer_dist || dist_z > z_outer_dist)
        {
            return max_mesh_size;
        }
        else if (dist_x <= x_inner_dist && dist_y <= y_inner_dist && dist_z <= z_inner_dist)
        {
            return min_mesh_size;
        }
        else {
            double maxxy = CGAL::max((dist_x - x_outer_dist) / (x_size * 0.5),
                (dist_y - y_outer_dist) / (y_size * 0.5));
            double max_dist_ratio = CGAL::max(maxxy, (dist_z - z_outer_dist) / (z_size * 0.5));
            return min_mesh_size + max_dist_ratio * (max_mesh_size-min_mesh_size);
        }
    }
};

Then I tested this code using my model, but the program keeps running without stopping. I doubt I made some mistake but could not solve it myself.

BTW I also looked the Lipschitz Sizing Field example, but in that example the mesh size is bigger inside and smaller outside, which is the opposite way. I swapped its min_size and max_size and the program also kept running without stopping.

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_outer_dist || dist_y > y_outer_dist || dist_z > z_outer_dist)
        {
            return max_mesh_size;
        }
        else if (dist_x <= x_inner_dist && dist_y <= y_inner_dist && dist_z <= z_inner_dist)
        {
            return min_mesh_size;
        }
        else {
            double maxxy = CGAL::max((dist_x - x_outer_dist) / (x_size * 0.5),
                (dist_y - y_outer_dist) / (y_size * 0.5));
            double max_dist_ratio = CGAL::max(maxxy, (dist_z - z_outer_dist) / (z_size * 0.5));
            return min_mesh_size + max_dist_ratio * (max_mesh_size-min_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 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.4, 0.75);
    size_1.SetMeshSize(10, 30);

    Mesh_criteria criteria(facet_angle = 25, facet_size = 10, facet_distance = 5,
        cell_radius_edge_ratio = 3, cell_size = size_1);

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

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

    return 0;
}

Environment

citystrawman commented 1 year ago

I tested again and this time I can get the mesh