search

CGAL / cgal

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

Hybrid domain meshing leads to erroneus elements close to the edges (Unexpected behavior) #8212

Closed f-p-b closed 1 month ago

f-p-b commented 1 month ago

Issue Details

When creating a volume mesh using a hybrid domain, elements are regularly missing or assigned to the wrong side of the implicit surface as can be seen in the figure.

The source code provided bellow should reproduce the issue if given the following polylines to protect: polylines.txt and the input .off file.

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.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/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Simplicial_mesh_cell_base_3.h>
#include <CGAL/Simplicial_mesh_vertex_base_3.h>
#include <CGAL/Surface_mesh/Surface_mesh.h>
#include <CGAL/IO/File_medit.h>
#include <CGAL/SMDS_3/Dump_c3t3.h>

#include "read_polylines.h"

typedef CGAL::Parallel_if_available_tag Concurrency_tag;

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

typedef CGAL::Surface_mesh<Point_3> SurfaceMesh;

// Implicit Domain
typedef CGAL::Labeled_mesh_domain_3<Kernel> Implicit_domain;

// Polyhedral Domain
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, Kernel> Polyhedron_domain;

using Subdomain_index = int;
using Surface_patch_index = unsigned char;
using Curve_index = char;
using Corner_index = short;
using Cb = CGAL::Simplicial_mesh_cell_base_3<Subdomain_index, Surface_patch_index>;
using Vb = CGAL::Simplicial_mesh_vertex_base_3<Kernel, Subdomain_index, Surface_patch_index,
                                               Curve_index, Corner_index>;
using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb, Concurrency_tag>;
using Triangulation = CGAL::Triangulation_3<Kernel, Tds>;

class Hybrid_domain {
    const Implicit_domain& implicit_domain;
    const Polyhedron_domain& polyhedron_domain;

    public:
    Hybrid_domain(const Implicit_domain& implicit_domain, const Polyhedron_domain& polyhedron_domain)
                  : implicit_domain(implicit_domain) , polyhedron_domain(polyhedron_domain) {}

    // types required by the 'MeshDomain_3' concept
    typedef int Surface_patch_index;
    typedef int Subdomain_index;
    typedef int Index;
    typedef Kernel R;
    typedef Kernel::Point_3 Point_3;
    typedef std::tuple<Point_3, Index, int> Intersection;

    CGAL::Bbox_3 bbox() const { return implicit_domain.bbox() + polyhedron_domain.bbox(); }

    struct Construct_initial_points {
        Construct_initial_points(const Hybrid_domain& domain) : r_domain_(domain) {}
        template<class OutputIterator> OutputIterator operator()(OutputIterator pts, const int n = 20) const {
            //construct initial points on implicit domain
            typedef Implicit_domain::Index Implicit_Index;
            std::vector<std::pair<Point_3, Implicit_Index> > implicit_points_vector;
            Implicit_domain::Construct_initial_points cstr_implicit_initial_points =
                r_domain_.implicit_domain.construct_initial_points_object();
            cstr_implicit_initial_points(std::back_inserter(implicit_points_vector), n/2);
            for(std::size_t i = 0, end = implicit_points_vector.size(); i<end; ++i) {
                *pts++ = std::make_pair(implicit_points_vector[i].first, 2);
            }
            //construct initial points on polyhedral domain
            typedef Polyhedron_domain::Index Polyhedron_Index;
            std::vector<std::pair<Point_3, Polyhedron_Index> > polyhedron_points_vector;
            Polyhedron_domain::Construct_initial_points cstr_polyhedron_initial_points =
                r_domain_.polyhedron_domain.construct_initial_points_object();
            cstr_polyhedron_initial_points(std::back_inserter(polyhedron_points_vector), n/2);
            for(std::size_t i = 0, end = polyhedron_points_vector.size(); i<end; ++i) {
                *pts++ = std::make_pair(polyhedron_points_vector[i].first, 1);
            }
            return pts;
        }

    private:
    const Hybrid_domain& r_domain_;
    }; // end Construct_initial_points_object

    Construct_initial_points construct_initial_points_object() const {
        return Construct_initial_points(*this);
    }

    struct Is_in_domain {
        Is_in_domain(const Hybrid_domain& domain) : r_domain_(domain) {}
        boost::optional<Subdomain_index> operator()(const Kernel::Point_3& p) const {
            boost::optional<Subdomain_index> implicit_subdomain_index = 
                r_domain_.implicit_domain.is_in_domain_object()(p);
            boost::optional<Subdomain_index> polyhedron_subdomain_index = 
                r_domain_.polyhedron_domain.is_in_domain_object()(p);

            if(!implicit_subdomain_index && polyhedron_subdomain_index)
                return 0;
            else
                return r_domain_.polyhedron_domain.is_in_domain_object()(p);

        }

        private:
        const Hybrid_domain& r_domain_;
    }; // end Is_in_domain_object

    Is_in_domain is_in_domain_object() const { return Is_in_domain(*this); }

    struct Construct_intersection {
        Construct_intersection(const Hybrid_domain& domain) : r_domain_(domain) {}
        template <typename Query> Intersection operator()(const Query& query) const {

            using boost::get;

            // intersection with implicit domain (within the polyhedral domain)
            Implicit_domain::Intersection implicit_inter =
                r_domain_.implicit_domain.construct_intersection_object()(query);
            // if found, return it
            if ( get<2>(implicit_inter) != 0 ) {
                const Point_3 inter_point = get<0>(implicit_inter);
                if ( r_domain_.polyhedron_domain.is_in_domain_object()(inter_point) ) // Inner implicit surface
                    return Intersection(inter_point, 4, get<2>(implicit_inter));
            }

            // intersection with polyhedral domain
            Polyhedron_domain::Intersection polyhedron_inter =
                r_domain_.polyhedron_domain.construct_intersection_object()(query);
            // if found, return it
            if ( get<2>(polyhedron_inter) != 0 ) {
                const Point_3 inter_point = get<0>(polyhedron_inter);
                if ( r_domain_.implicit_domain.is_in_domain_object()(inter_point) ) // Scaffold
                    return Intersection(inter_point, 4, get<2>(polyhedron_inter));
                else // Void
                    return Intersection(inter_point, 5, get<2>(polyhedron_inter));
            }

            //no intersection found
            return Intersection();
        }

        private:
        const Hybrid_domain& r_domain_;
    }; // end Construct_intersection_object

    Construct_intersection construct_intersection_object() const {
        return Construct_intersection(*this);
    }

    //Index types converters
    Index index_from_surface_patch_index(const Surface_patch_index& index) const
        { return index; }
    Index index_from_subdomain_index(const Subdomain_index& index) const
        { return index; }
    Surface_patch_index surface_patch_index(const Index& index) const
        { return index; }
    Subdomain_index subdomain_index(const Index& index) const
        { return index; }
}; // end Hybrid_domain class

typedef CGAL::Mesh_domain_with_polyline_features_3<Hybrid_domain> H_domain;
typedef CGAL::Mesh_triangulation_3<H_domain, CGAL::Default, Concurrency_tag>::type H_Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<H_Tr> C3t3;
typedef CGAL::Mesh_criteria_3<H_Tr>     H_Mesh_criteria;
typedef H_Mesh_criteria::Edge_criteria  H_Edge_criteria;
typedef H_Mesh_criteria::Facet_criteria H_Facet_criteria;
typedef H_Mesh_criteria::Cell_criteria  H_Cell_criteria;

FT my_implicit_surface (const Point_3& p) {
    const double x2 = p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
    double a = (x2+y2+z2-0.8); // Outer shape

    double scaling = (2*3.14159265359)/1.5;
    double b = sin(scaling*p.x())*cos(scaling*p.y()) +
        sin(scaling*p.y())*cos(scaling*p.z()) +
        sin(scaling*p.z())*cos(scaling*p.x()) -
         0;

    if (a > b) {return a;} else {return b;}
}

int main() {
  const std::string fname = CGAL::data_file_path("cube.off");

  // Create input polyhedron
  Polyhedron polyhedron;
  std::ifstream input(fname);
  input >> polyhedron;
  if(input.bad()){
    std::cerr << "Error: Cannot read file " <<  fname << std::endl;
    return EXIT_FAILURE;
  }
  input.close();

  // Domains
  Polyhedron_domain polyhedron_domain(polyhedron);
  Implicit_domain sphere_domain =
    Implicit_domain::create_implicit_mesh_domain(my_implicit_surface,
                                                 Kernel::Sphere_3(Point_3(0, 0, 0), FT(9)));
  H_domain domain(sphere_domain, polyhedron_domain);

  // Polyline features
  const char* lines_fname = "polylines.txt";
  std::vector<std::vector<Point_3> > featured_curves;
  if (!read_polylines(lines_fname, featured_curves)) {
    std::cerr << "Error: Cannot read file " << lines_fname << std::endl;
    return EXIT_FAILURE;
  }

  // Add features for protection
  domain.add_features(featured_curves.begin(), featured_curves.end());

  // Criteria
  H_Edge_criteria edge_criteria(0.05);
  H_Facet_criteria facet_criteria(30, 0.05, 0.025); // angle, size, approximation
  H_Cell_criteria cell_criteria(2, 0.05); // radius-edge ratio, size
  H_Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria);

  // Mesh generation (without optimization)
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
                                      CGAL::parameters::no_perturb().no_exude());
  // Output
  dump_c3t3(c3t3, "out");

  return EXIT_SUCCESS;
}

Environment

sloriot commented 1 month ago

With cgal 5.6.x branch I have the following: image

using https://gist.github.com/sloriot/f7ea17c69c722acc54bef824e1017b57

If I add .manifold() I have:

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
File: /home/sloriot/CGAL/git/integration/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
Line: 859
Explanation: (-0.3633066446675759 -0.5 -0.29269817778953555 0) is already inserted on surface.

Aborted
f-p-b commented 1 month ago

Although the spheres representing the constrained edges obscure most of the affected areas in the image, based on the few spots with issues visible (marked by the red arrows) it seems to reproduce the issue I am encountering. The two tetra sticking out at the top should not be there as they don't belong to the domain being meshed and at the bottom corner there is an element missing from the domain:

When running with .manifold() I have: 

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
File: /mnt/c/Workdir/Develop/TestBox/cgal_test/libs/CGAL-5.6/include/CGAL/Mesh_3/Refine_facets_3.h
Line: 1500
Explanation: Mesh_3 ERROR: A facet is not in conflict with its refinement point!
Debugging information:
  Facet: (0x55f6743245a0, 3) = (0.5 0.24802094429145205 0.31521894221984847 0, 0.5 0.24802094735800306 0.31521895230952635 0, 0.5 0.24802094229966587 0.315218953368535 0)
  Dual: Segment(0.50239510072181526 0.24802094379119005 0.31521894788273952 , 0.49527770566129697 0.2480344687443283 0.31526551456807811)
  Refinement point: 0.5 0.2480254951227866 0.31523461820885157 0
  Cells adjacent to facet:
    ( 0.5 0.248021 0.315219 0 , 0.5 0.248021 0.315219 0 , 0.5 0.248021 0.315219 0 , 0.500025 0.254315 0.355228 0.00164025 )
    ( 0.499965 0.287755 0.323072 0.00164025 , 0.5 0.248021 0.315219 0 , 0.5 0.248021 0.315219 0 , 0.5 0.248021 0.315219 0 )

Aborted

Am I making a mistake when defining my hybrid domain or polylines, or is it more an issue on CGAL's end?

lrineau commented 1 month ago

Anonymous (@Another-Anonymous-User, really?! why that name?),

I have reproduced your issue and I propose a fix. I published my version as a branch in my fork: fix-issue-8212-hydrib-domain.

You can see the four commits I have added, and the diff compared to your version: https://github.com/lrineau/cgal/compare/e088fc4dfbc...lrineau:cgal:fix-issue-8212-hydrib-domain

https://github.com/lrineau/cgal/compare/e088fc4dfbc is your version of the code.

That issue is not really a bug, as far as we know. That is the expected behavior of Mesh_3, given the data, the domain, and the mesh criteria. I will close this issue, but we can still discuss in it.

f-p-b commented 1 month ago

Thanks you for helping out!

Regarding your edit to line 102 (return 2 instead of 0). I was returning a zero because I don’t need the mesh of both domains and assumed for large problems keeping only the required domain will reduce the required resources noticeably since only the mesh for the domain of interest is kept in memory. Given you changed it, are there risks in removing the domain by returning a 0 instead of a 2?

Regarding the adjustments to the relative_error_bound and facet_distance. I have tested some other cases and adjusting them (especially the relative_error_bound) has helped fix similar issues. I am still a bit puzzled by it though. At which stage does relative_error_bound come into play? Because sure, I get its the intersection between the implicit surface and query segments. However, how closely the surface of the mesh approximates the implicit surface is controlled by facet_distance if I understand correctly. Does this means that it’s the query made to determine at which side of the surface an element lies? Or is it some other query at some other stage of the process?

why that name?

@lrineau, It started as a joke and then I was too lazy to change it if I am honest ^_^'. Will probably get to changing it one of these days…

lrineau commented 1 month ago

Thanks you for helping out!

Regarding your edit to line 102 (return 2 instead of 0). I was returning a zero because I don’t need the mesh of both domains and assumed for large problems keeping only the required domain will reduce the required resources noticeably since only the mesh for the domain of interest is kept in memory. Given you changed it, are there risks in removing the domain by returning a 0 instead of a 2?

There is no risk, but the right value is not return 0. I have pushed a new commit, see: https://github.com/lrineau/cgal/compare/e088fc4dfbc...lrineau:cgal:fix-issue-8212-hydrib-domain

With return 0, the result is a non-null boost::optional<int> containing the value 0. And that is forbidden by the CGAL Mesh_3 documentation (the subdomain index cannot be null).

With return boost::none, the result is a null boost::optional, that is the way to tell the Mesh_3 engine that the query point is outside the mesh domain.

Regarding the adjustments to the relative_error_bound and facet_distance. I have tested some other cases and adjusting them (especially the relative_error_bound) has helped fix similar issues. I am still a bit puzzled by it though. At which stage does relative_error_bound come into play? Because sure, I get its the intersection between the implicit surface and query segments. However, how closely the surface of the mesh approximates the implicit surface is controlled by facet_distance if I understand correctly. Does this means that it’s the query made to determine at which side of the surface an element lies? Or is it some other query at some other stage of the process?

The code you are working on has two distinct parts that are following the abstraction used in CGAL Mesh_3:

The parameter facet_distance is in the mesh criteria, and only used by the mesh engine to decide if a given 3D triangle is "close enough" to the actual surface of the mesh domain.

The parameter relative_error_bound is a parameter for the mesh domain class, and rules the bisection algorithm used to compute intersections between a query segment and a given implicit surface. If the relative_error_bound is too big, the mesh domain implementation will compute points on the implicit surface, with not enough precision. That means that the surface the mesh engine is constructing will be noisy, or blurry. And then, with a small facet_distance, the mesh engine has close to zero chance to converge to a state where the facet_distance criterion is satisfied for all 3D triangles of the surface mesh.

why that name?

@lrineau, It started as a joke and then I was too lazy to change it if I am honest ^_^'. Will probably get to changing it one of these days…

By the way, what are you working on? Unless you really want to stay anonymous, maybe tell us more about you and the actual goal of the CGAL code you are working on. Statistically, you are probably a student, probably a PhD student, who needs that CGAL code for a research purpose.

f-p-b commented 1 month ago

Thanks for the additional explanations.

By the way, what are you working on?

This is what I am working on: https://github.com/tpms-lattice/ASLI/tree/develop

lrineau commented 4 weeks ago

why that name?

@lrineau, It started as a joke and then I was too lazy to change it if I am honest ^_^'. Will probably get to changing it one of these days…

You are probably https://www.mech.kuleuven.be/en/bme/bme-people/00116401

I wonder if your lab has had collaboration with the CGAL project in the past. Maybe we can help you more about your code. Have you already tried CGAL 3d periodic meshes? Maybe you could discuss with @MaelRL and @MoniqueTeillaud who are still active contributors to the CGAL project.

There was actually a collaboration between CGAL (actually @MoniqueTeillaud and KU Leuven in 2008-2009:

f-p-b commented 3 weeks ago

Bingo, I see you were really curious.

The group I am in (Biomech Research Unit) has not had any collaborations with the CGAL project as far as I know.

I did look into using the 3D periodic meshes of CGAL a while back for cases were the structure is uniform but was having some issues with elements missing from the triangulation. I was in contact about it with MaelRL trough stackoverflow but even after he made some improvements it was not completely eliminated. So I did not really pursue it further at the time since the code main use mode (at least for us) is non-uniform structures. Still want to try it out again at some point now that CGAL has had some updates as it would make generating uniform infills much faster. But I have to make some time for that.