Problems with Sharp Features in 3d Mesh Generation

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

Problems with Sharp Features in 3d Mesh Generation

MariusZeinhofer
Dear all,

I am trying to mesh a 3d mesh using the create_implicit_mesh_domain method
and I am running into problems with the sharp representations of 1d
features.

My 3d object is a set union of a cylinder with a segement of another
cylinder of larger radius (see the image attached), and I want these edges
to be sharp. I am working my way along the example 3.6.2 in the CGAL manual
(Implicit Domain with 1D Features, the intersecting balls). The problem I'm
getting is that my edges are "too sharp" represented which makes the mesh
too big for my purposes. The mesh I get is extremly fine at the edges, see
the following picture:

<http://cgal-discuss.949826.n4.nabble.com/file/t376248/CylinderSmall.png>

What confuses me is, that not all of the sharp edges are problematic and
also if I am just meshing the smaller cylinder I'm not getting this issue at
the top and bottom of the cylinder, even though I am using the same code.

My mesh criteria are:
Mesh_criteria criteria(edge_size = 0.29, facet_angle = 25, facet_size = 0.3,
cell_radius_edge_ratio = 3.5, cell_size = .5, facet_distance = 0.4);

The way I'm introducing the sharp features is:

typedef std::vector<Point>        Polyline_3;
typedef std::list<Polyline_3>       Polylines;

for (int i = 0; i < 360; ++i)
    {
        Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i *
CGAL_PI / 360));
        polyline.push_back(p);
    }
    polyline.push_back(polyline.front()); // close the line
    polylines.push_back(polyline);

and in a similar fashion for the other sharp features.

Thank you very much for your help in advance!
Cheers
Marius



--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

Sebastien Loriot (GeometryFactory)
Could you provide a minimal example that we can run so that we can
reproduce the issue?

I don't see any issue in the snippet provided.

Thanks,

Sebastien.

On 9/7/20 8:33 AM, MariusZeinhofer wrote:

> Dear all,
>
> I am trying to mesh a 3d mesh using the create_implicit_mesh_domain method
> and I am running into problems with the sharp representations of 1d
> features.
>
> My 3d object is a set union of a cylinder with a segement of another
> cylinder of larger radius (see the image attached), and I want these edges
> to be sharp. I am working my way along the example 3.6.2 in the CGAL manual
> (Implicit Domain with 1D Features, the intersecting balls). The problem I'm
> getting is that my edges are "too sharp" represented which makes the mesh
> too big for my purposes. The mesh I get is extremly fine at the edges, see
> the following picture:
>
> <http://cgal-discuss.949826.n4.nabble.com/file/t376248/CylinderSmall.png>
>
> What confuses me is, that not all of the sharp edges are problematic and
> also if I am just meshing the smaller cylinder I'm not getting this issue at
> the top and bottom of the cylinder, even though I am using the same code.
>
> My mesh criteria are:
> Mesh_criteria criteria(edge_size = 0.29, facet_angle = 25, facet_size = 0.3,
> cell_radius_edge_ratio = 3.5, cell_size = .5, facet_distance = 0.4);
>
> The way I'm introducing the sharp features is:
>
> typedef std::vector<Point>        Polyline_3;
> typedef std::list<Polyline_3>       Polylines;
>
> for (int i = 0; i < 360; ++i)
>      {
>          Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i *
> CGAL_PI / 360));
>          polyline.push_back(p);
>      }
>      polyline.push_back(polyline.front()); // close the line
>      polylines.push_back(polyline);
>
> and in a similar fashion for the other sharp features.
>
> Thank you very much for your help in advance!
> Cheers
> Marius
>
>
>
> --
> Sent from: http://cgal-discuss.949826.n4.nabble.com/
>

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

andreas.fabri
Qu'est-ce qu'on peut faire avec Surface_mesher qu'on arrive pas a faire avec Mesh_3?

Meme si on peut tout faire avec Mesh_3 ca vaut peut-etre d'avoir un chapitre à part?
Et peut-etre deriver Surface_mesher de Mesh_3 ?


On 9/7/2020 9:24 AM, "Sebastien Loriot (GeometryFactory)" ([hidden email] via cgal-discuss Mailing List) wrote:
Could you provide a minimal example that we can run so that we can
reproduce the issue?

I don't see any issue in the snippet provided.

Thanks,

Sebastien.

On 9/7/20 8:33 AM, MariusZeinhofer wrote:
Dear all,

I am trying to mesh a 3d mesh using the create_implicit_mesh_domain method
and I am running into problems with the sharp representations of 1d
features.

My 3d object is a set union of a cylinder with a segement of another
cylinder of larger radius (see the image attached), and I want these edges
to be sharp. I am working my way along the example 3.6.2 in the CGAL manual
(Implicit Domain with 1D Features, the intersecting balls). The problem I'm
getting is that my edges are "too sharp" represented which makes the mesh
too big for my purposes. The mesh I get is extremly fine at the edges, see
the following picture:

<http://cgal-discuss.949826.n4.nabble.com/file/t376248/CylinderSmall.png>

What confuses me is, that not all of the sharp edges are problematic and
also if I am just meshing the smaller cylinder I'm not getting this issue at
the top and bottom of the cylinder, even though I am using the same code.

My mesh criteria are:
Mesh_criteria criteria(edge_size = 0.29, facet_angle = 25, facet_size = 0.3,
cell_radius_edge_ratio = 3.5, cell_size = .5, facet_distance = 0.4);

The way I'm introducing the sharp features is:

typedef std::vector<Point>        Polyline_3;
typedef std::list<Polyline_3>       Polylines;

for (int i = 0; i < 360; ++i)
     {
         Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i *
CGAL_PI / 360));
         polyline.push_back(p);
     }
     polyline.push_back(polyline.front()); // close the line
     polylines.push_back(polyline);

and in a similar fashion for the other sharp features.

Thank you very much for your help in advance!
Cheers
Marius



--
Sent from: http://cgal-discuss.949826.n4.nabble.com/


-- 
Andreas Fabri, PhD
Chief Officer, GeometryFactory
Editor, The CGAL Project

phone: +33.492.954.912    skype: andreas.fabri

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

MariusZeinhofer
In reply to this post by Sebastien Loriot (GeometryFactory)
Of course.
And thank you very much for your help!

Cheers Marius


#include <vector>

#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/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <cmath>

// Kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

// Domain
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT(Function)(const Point&);
typedef
CGAL::Mesh_domain_with_polyline_features_3<CGAL::Labeled_mesh_domain_3&lt;K>
> Mesh_domain;

// Polyline
typedef std::vector<Point>        Polyline_3;
typedef std::list<Polyline_3>       Polylines;

#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,
Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3
 Mesh_criteria;

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

// Functions
// describes cylinder: \{ (x,y,z) \in \mathbb{R}^3 \mid |(y,z)|_euclid <=
10, -15 <= x <= 15  \}
FT cylinder_function(const Point& p)
{
        const FT x2 = p.x() * p.x();
        const FT y2 = p.y() * p.y();
        const FT z2 = p.z() * p.z();
       
        FT aux = y2 + z2 - 100;
       
        if (p.x() > 15 && aux < 0)
                aux = -aux;

        if (p.x() < -15 && aux < 0)
                aux = -aux;

        return aux;
}

// section of cylinder of larger radius: \{ (x,y,z) \in \mathbb{R}^3 \mid
|(y,z)|_euclid <= 12, -15 <= x <= 15, angle in (y,z) plane between 0 and
\pi/6  \}
FT cylinder_section_function(const Point& p)
{
        const FT x2 = p.x() * p.x();
        const FT y2 = p.y() * p.y();
        const FT z2 = p.z() * p.z();
        FT aux = y2 + z2 - 144;
        if (p.x() > 15 && aux < 0)
                aux = -aux;

        if (p.x() < -15 && aux < 0)
                aux = -aux;

        if ( (aux < 0 && !(0 < std::atan(p.z() / p.y()) && std::atan(p.z() / p.y())
< CGAL_PI / 6.) ) || p.y() < 0 || p.z() < 0 )
                aux = 1;

        return aux;
}

FT set_union(const Point& p)
{
        if (cylinder_section_function(p) < 0 || cylinder_function(p) < 0)
                return -1;
        else
                return 1;
}

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

//---------------------------------------------------------------------generate-.mesh-file-----------------------------------------------------------------------------------------------------------//

        // Domain (Warning: Sphere_3 constructor uses squared radius !)
        Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(set_union,
K::Sphere_3(CGAL::ORIGIN, 400));

        // Mesh criteria
        Mesh_criteria criteria(edge_size = 0.29, facet_angle = 25, facet_size =
0.3, cell_radius_edge_ratio = 3.5, cell_size = .5, facet_distance = 0.4);

        // Create edge that we want to preserve
        Polylines polylines;
        //Polyline_3& polyline = polylines.front();
        //Polyline_3& polyline2 = polylines.back();

        Polyline_3 polyline;
        Polyline_3 polyline2;
        Polyline_3 polyline3;
        Polyline_3 polyline4;
        Polyline_3 polyline5;
        Polyline_3 polyline6;
        Polyline_3 polyline7;
        Polyline_3 polyline8;

        //top of cylinder
        for (int i = 0; i < 360; ++i)
        {
                Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
/ 360));
                polyline.push_back(p);
        }
        polyline.push_back(polyline.front()); // close the line
        polylines.push_back(polyline);

        //bottom of cylinder
        for (int i = 0; i < 360; ++i)
        {
                Point p(-15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
/ 360));
                polyline2.push_back(p);
        }
        polyline2.push_back(polyline2.front()); // close the line
        polylines.push_back(polyline2);


        //fixateur line: (-15 < x < 15, y = 12, z = 0)
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * (30.) / 99., 12., 0.);
                polyline3.push_back(p);
        }
        polylines.push_back(polyline3);

        //fixateur line: (-15 < x < 15, y = 10, z = 0)
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * ((30.) / 99.), 10., 0.);
                polyline4.push_back(p);
        }
        polylines.push_back(polyline4);


        //fixateur line: (-15 < x < 15, y = 12 * cos(pi/6), z = 12 * sin(pi/6) )
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * ((30.) / 99.), 12 * std::cos(CGAL_PI / 6), 12 *
std::sin(CGAL_PI / 6));
                polyline5.push_back(p);
        }
        polylines.push_back(polyline5);

        //fixateur line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
std::sin(CGAL_PI / 6));
                polyline6.push_back(p);
        }
        polylines.push_back(polyline6);


        //top of cylinder segment
        for (int i = 0; i < 61; ++i)
        {
                Point p(15.0, 12. * std::cos(i * CGAL_PI / 360), 12. * std::sin(i *
CGAL_PI / 360));
                polyline7.push_back(p);
        }
        polylines.push_back(polyline7);

        //bottom of cylinder segment
        for (int i = 0; i < 61; ++i)
        {
                Point p(-15.0, 12. * std::cos(i * CGAL_PI / 360), 12. * std::sin(i *
CGAL_PI / 360));
                polyline8.push_back(p);
        }
        polylines.push_back(polyline8);

        // Insert edge in domain
        domain.add_features(polylines.begin(), polylines.end());

        // Mesh generation without feature preservation
        C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
CGAL::parameters::no_features());

        // Output
        std::ofstream medit_file("out-no-protection.mesh");
        c3t3.output_to_medit(medit_file);
        medit_file.close();
        c3t3.clear();

        // Mesh generation with feature preservation
        c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);

        // Output
        medit_file.open("out.mesh");
        c3t3.output_to_medit(medit_file);
        medit_file.close();

        return 0;
}



--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

Sebastien Loriot (GeometryFactory)
Except if I made a mistake dumping your polylines, they are not what you
expect them to be as you can see in the attached screenshot. If you have
2 very close (or identical polylines) then the protection algorithm will
try to seperate them which is not possible.

Could you please check them?

Thanks,

Sebastien.

PS: Code used to dump them in a format readable by CGAL the polyhedron demo

int ggg=0;
for (auto p : polylines)
{
   std::ofstream outpol("poly_"+std::to_string(ggg)+".polylines.txt");
   outpol << p.size();
   for (auto oo : p)
     outpol << " " << oo;
   outpol << "\n";
   ++ggg;
}


On 9/7/20 10:04 AM, MariusZeinhofer wrote:

> Of course.
> And thank you very much for your help!
>
> Cheers Marius
>
>
> #include <vector>
>
> #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/Mesh_domain_with_polyline_features_3.h>
> #include <CGAL/make_mesh_3.h>
> #include <cmath>
>
> // Kernel
> typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
>
> // Domain
> typedef K::FT FT;
> typedef K::Point_3 Point;
> typedef FT(Function)(const Point&);
> typedef
> CGAL::Mesh_domain_with_polyline_features_3<CGAL::Labeled_mesh_domain_3&lt;K>
>> Mesh_domain;
>
> // Polyline
> typedef std::vector<Point>        Polyline_3;
> typedef std::list<Polyline_3>       Polylines;
>
> #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,
> Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
>
> // Criteria
> typedef CGAL::Mesh_criteria_3
>   Mesh_criteria;
>
> // To avoid verbose function and named parameters call
> using namespace CGAL::parameters;
>
> // Functions
> // describes cylinder: \{ (x,y,z) \in \mathbb{R}^3 \mid |(y,z)|_euclid <=
> 10, -15 <= x <= 15  \}
> FT cylinder_function(const Point& p)
> {
> const FT x2 = p.x() * p.x();
> const FT y2 = p.y() * p.y();
> const FT z2 = p.z() * p.z();
>
> FT aux = y2 + z2 - 100;
>
> if (p.x() > 15 && aux < 0)
> aux = -aux;
>
> if (p.x() < -15 && aux < 0)
> aux = -aux;
>
> return aux;
> }
>
> // section of cylinder of larger radius: \{ (x,y,z) \in \mathbb{R}^3 \mid
> |(y,z)|_euclid <= 12, -15 <= x <= 15, angle in (y,z) plane between 0 and
> \pi/6  \}
> FT cylinder_section_function(const Point& p)
> {
> const FT x2 = p.x() * p.x();
> const FT y2 = p.y() * p.y();
> const FT z2 = p.z() * p.z();
> FT aux = y2 + z2 - 144;
> if (p.x() > 15 && aux < 0)
> aux = -aux;
>
> if (p.x() < -15 && aux < 0)
> aux = -aux;
>
> if ( (aux < 0 && !(0 < std::atan(p.z() / p.y()) && std::atan(p.z() / p.y())
> < CGAL_PI / 6.) ) || p.y() < 0 || p.z() < 0 )
> aux = 1;
>
> return aux;
> }
>
> FT set_union(const Point& p)
> {
> if (cylinder_section_function(p) < 0 || cylinder_function(p) < 0)
> return -1;
> else
> return 1;
> }
>
> int main(int argc, char* argv[])
> {
>
> //---------------------------------------------------------------------generate-.mesh-file-----------------------------------------------------------------------------------------------------------//
>
> // Domain (Warning: Sphere_3 constructor uses squared radius !)
> Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(set_union,
> K::Sphere_3(CGAL::ORIGIN, 400));
>
> // Mesh criteria
> Mesh_criteria criteria(edge_size = 0.29, facet_angle = 25, facet_size =
> 0.3, cell_radius_edge_ratio = 3.5, cell_size = .5, facet_distance = 0.4);
>
> // Create edge that we want to preserve
> Polylines polylines;
> //Polyline_3& polyline = polylines.front();
> //Polyline_3& polyline2 = polylines.back();
>
> Polyline_3 polyline;
> Polyline_3 polyline2;
> Polyline_3 polyline3;
> Polyline_3 polyline4;
> Polyline_3 polyline5;
> Polyline_3 polyline6;
> Polyline_3 polyline7;
> Polyline_3 polyline8;
>
> //top of cylinder
> for (int i = 0; i < 360; ++i)
> {
> Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
> / 360));
> polyline.push_back(p);
> }
> polyline.push_back(polyline.front()); // close the line
> polylines.push_back(polyline);
>
> //bottom of cylinder
> for (int i = 0; i < 360; ++i)
> {
> Point p(-15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
> / 360));
> polyline2.push_back(p);
> }
> polyline2.push_back(polyline2.front()); // close the line
> polylines.push_back(polyline2);
>
>
> //fixateur line: (-15 < x < 15, y = 12, z = 0)
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * (30.) / 99., 12., 0.);
> polyline3.push_back(p);
> }
> polylines.push_back(polyline3);
>
> //fixateur line: (-15 < x < 15, y = 10, z = 0)
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * ((30.) / 99.), 10., 0.);
> polyline4.push_back(p);
> }
> polylines.push_back(polyline4);
>
>
> //fixateur line: (-15 < x < 15, y = 12 * cos(pi/6), z = 12 * sin(pi/6) )
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * ((30.) / 99.), 12 * std::cos(CGAL_PI / 6), 12 *
> std::sin(CGAL_PI / 6));
> polyline5.push_back(p);
> }
> polylines.push_back(polyline5);
>
> //fixateur line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
> std::sin(CGAL_PI / 6));
> polyline6.push_back(p);
> }
> polylines.push_back(polyline6);
>
>
> //top of cylinder segment
> for (int i = 0; i < 61; ++i)
> {
> Point p(15.0, 12. * std::cos(i * CGAL_PI / 360), 12. * std::sin(i *
> CGAL_PI / 360));
> polyline7.push_back(p);
> }
> polylines.push_back(polyline7);
>
> //bottom of cylinder segment
> for (int i = 0; i < 61; ++i)
> {
> Point p(-15.0, 12. * std::cos(i * CGAL_PI / 360), 12. * std::sin(i *
> CGAL_PI / 360));
> polyline8.push_back(p);
> }
> polylines.push_back(polyline8);
>
> // Insert edge in domain
> domain.add_features(polylines.begin(), polylines.end());
>
> // Mesh generation without feature preservation
> C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
> CGAL::parameters::no_features());
>
> // Output
> std::ofstream medit_file("out-no-protection.mesh");
> c3t3.output_to_medit(medit_file);
> medit_file.close();
> c3t3.clear();
>
> // Mesh generation with feature preservation
> c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
>
> // Output
> medit_file.open("out.mesh");
> c3t3.output_to_medit(medit_file);
> medit_file.close();
>
> return 0;
> }
>
>
>
> --
> Sent from: http://cgal-discuss.949826.n4.nabble.com/
>
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss



Screenshot from 2020-09-08 09-14-05.png (326K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

MariusZeinhofer
Thanks for your reply. I modified my code in the wrong way prior to posting
here therefore these polylines are not what I wanted them to be. I have
tried again to provide a minimal example which illustrates the issue at the
end of this message, which produces the following image:

<http://cgal-discuss.949826.n4.nabble.com/file/t376248/CylinderProblem.jpg>

But from your answer I may have understood what the problem is:

Say I parametrize two polylines:

        //top of cylinder
        for (int i = 31; i < 360; ++i)
        {
                Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
/ 180));
                polyline.push_back(p);
        }
        polylines.push_back(polyline);

        //line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
std::sin(CGAL_PI / 6));
                polyline6.push_back(p);
        }
        polylines.push_back(polyline6);

those two lines (almost) meet at one point, namely: (15.0, cos(pi/6),
sin(pi/6) ) which causes problems with the protection algorithm? So if I
have several sharp 1D features that I want to preserve I should be careful
with my polylines not sharing points?

Thank you again and sorry for the faulty minimal example.
Cheers Marius

Minimal Example:

#include <vector>

#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/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <cmath>

// Kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

// Domain
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT(Function)(const Point&);
typedef
CGAL::Mesh_domain_with_polyline_features_3<CGAL::Labeled_mesh_domain_3&lt;K>
> Mesh_domain;

// Polyline
typedef std::vector<Point>        Polyline_3;
typedef std::list<Polyline_3>       Polylines;

#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,
Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3
 Mesh_criteria;

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

// Functions
// describes cylinder: \{ (x,y,z) \in \mathbb{R}^3 \mid |(y,z)|_euclid <=
10, -15 <= x <= 15  \}
FT cylinder_function(const Point& p)
{
        const FT x2 = p.x() * p.x();
        const FT y2 = p.y() * p.y();
        const FT z2 = p.z() * p.z();
       
        FT aux = y2 + z2 - 100;
       
        if (p.x() > 15 && aux < 0)
                aux = -aux;

        if (p.x() < -15 && aux < 0)
                aux = -aux;

        return aux;
}

// section of cylinder of larger radius: \{ (x,y,z) \in \mathbb{R}^3 \mid
|(y,z)|_euclid <= 12, -15 <= x <= 15, angle in (y,z) plane between 0 and
\pi/6  \}
FT cylinder_section_function(const Point& p)
{
        const FT x2 = p.x() * p.x();
        const FT y2 = p.y() * p.y();
        const FT z2 = p.z() * p.z();
        FT aux = y2 + z2 - 144;
        if (p.x() > 15 && aux < 0)
                aux = -aux;

        if (p.x() < -15 && aux < 0)
                aux = -aux;

        if ( (aux < 0 && !(0 < std::atan(p.z() / p.y()) && std::atan(p.z() / p.y())
< CGAL_PI / 6.) ) || p.y() < 0 || p.z() < 0 )
                aux = 1;

        return aux;
}

FT set_union(const Point& p)
{
        if (cylinder_section_function(p) < 0 || cylinder_function(p) < 0)
                return -1;
        else
                return 1;
}

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

//---------------------------------------------------------------------generate-.mesh-file-----------------------------------------------------------------------------------------------------------//

        // Domain (Warning: Sphere_3 constructor uses squared radius !)
        Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(set_union,
K::Sphere_3(CGAL::ORIGIN, 400));

        // Mesh criteria
        Mesh_criteria criteria(edge_size = 0.49, facet_angle = 25, facet_size =
0.5, cell_radius_edge_ratio = 3.5, cell_size = .8, facet_distance = 0.6);

        // Create edge that we want to preserve
        Polylines polylines;
        //Polyline_3& polyline = polylines.front();
        //Polyline_3& polyline2 = polylines.back();

        Polyline_3 polyline;
        Polyline_3 polyline6;

        //top of cylinder
        for (int i = 31; i < 360; ++i)
        {
                Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
/ 180));
                polyline.push_back(p);
        }
        polylines.push_back(polyline);

        //fixateur line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
        for (int i = 0; i < 100; ++i)
        {
                Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
std::sin(CGAL_PI / 6));
                polyline6.push_back(p);
        }
        polylines.push_back(polyline6);

        // Insert edge in domain
        domain.add_features(polylines.begin(), polylines.end());

        // Mesh generation without feature preservation
        C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
CGAL::parameters::no_features());

        // Output
        std::ofstream medit_file("out-no-protection.mesh");
        c3t3.output_to_medit(medit_file);
        medit_file.close();
        c3t3.clear();

        // Mesh generation with feature preservation
        c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);

        // Output
        medit_file.open("out.mesh");
        c3t3.output_to_medit(medit_file);
        medit_file.close();

        return 0;
}



--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Problems with Sharp Features in 3d Mesh Generation

Sebastien Loriot (GeometryFactory)
You have to make the end points of those two polylines identical,
otherwise small protecting balls will be inserted so as to separate
then. Then there is a interpolation of the size of the protecting balls
along polylines.

Best,

Sebastien.

On 9/8/20 1:53 PM, MariusZeinhofer wrote:

> Thanks for your reply. I modified my code in the wrong way prior to posting
> here therefore these polylines are not what I wanted them to be. I have
> tried again to provide a minimal example which illustrates the issue at the
> end of this message, which produces the following image:
>
> <http://cgal-discuss.949826.n4.nabble.com/file/t376248/CylinderProblem.jpg>
>
> But from your answer I may have understood what the problem is:
>
> Say I parametrize two polylines:
>
>          //top of cylinder
> for (int i = 31; i < 360; ++i)
> {
> Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
> / 180));
> polyline.push_back(p);
> }
> polylines.push_back(polyline);
>
> //line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
> std::sin(CGAL_PI / 6));
> polyline6.push_back(p);
> }
> polylines.push_back(polyline6);
>
> those two lines (almost) meet at one point, namely: (15.0, cos(pi/6),
> sin(pi/6) ) which causes problems with the protection algorithm? So if I
> have several sharp 1D features that I want to preserve I should be careful
> with my polylines not sharing points?
>
> Thank you again and sorry for the faulty minimal example.
> Cheers Marius
>
> Minimal Example:
>
> #include <vector>
>
> #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/Mesh_domain_with_polyline_features_3.h>
> #include <CGAL/make_mesh_3.h>
> #include <cmath>
>
> // Kernel
> typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
>
> // Domain
> typedef K::FT FT;
> typedef K::Point_3 Point;
> typedef FT(Function)(const Point&);
> typedef
> CGAL::Mesh_domain_with_polyline_features_3<CGAL::Labeled_mesh_domain_3&lt;K>
>> Mesh_domain;
>
> // Polyline
> typedef std::vector<Point>        Polyline_3;
> typedef std::list<Polyline_3>       Polylines;
>
> #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,
> Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3;
>
> // Criteria
> typedef CGAL::Mesh_criteria_3
>   Mesh_criteria;
>
> // To avoid verbose function and named parameters call
> using namespace CGAL::parameters;
>
> // Functions
> // describes cylinder: \{ (x,y,z) \in \mathbb{R}^3 \mid |(y,z)|_euclid <=
> 10, -15 <= x <= 15  \}
> FT cylinder_function(const Point& p)
> {
> const FT x2 = p.x() * p.x();
> const FT y2 = p.y() * p.y();
> const FT z2 = p.z() * p.z();
>
> FT aux = y2 + z2 - 100;
>
> if (p.x() > 15 && aux < 0)
> aux = -aux;
>
> if (p.x() < -15 && aux < 0)
> aux = -aux;
>
> return aux;
> }
>
> // section of cylinder of larger radius: \{ (x,y,z) \in \mathbb{R}^3 \mid
> |(y,z)|_euclid <= 12, -15 <= x <= 15, angle in (y,z) plane between 0 and
> \pi/6  \}
> FT cylinder_section_function(const Point& p)
> {
> const FT x2 = p.x() * p.x();
> const FT y2 = p.y() * p.y();
> const FT z2 = p.z() * p.z();
> FT aux = y2 + z2 - 144;
> if (p.x() > 15 && aux < 0)
> aux = -aux;
>
> if (p.x() < -15 && aux < 0)
> aux = -aux;
>
> if ( (aux < 0 && !(0 < std::atan(p.z() / p.y()) && std::atan(p.z() / p.y())
> < CGAL_PI / 6.) ) || p.y() < 0 || p.z() < 0 )
> aux = 1;
>
> return aux;
> }
>
> FT set_union(const Point& p)
> {
> if (cylinder_section_function(p) < 0 || cylinder_function(p) < 0)
> return -1;
> else
> return 1;
> }
>
> int main(int argc, char* argv[])
> {
>
> //---------------------------------------------------------------------generate-.mesh-file-----------------------------------------------------------------------------------------------------------//
>
> // Domain (Warning: Sphere_3 constructor uses squared radius !)
> Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(set_union,
> K::Sphere_3(CGAL::ORIGIN, 400));
>
> // Mesh criteria
> Mesh_criteria criteria(edge_size = 0.49, facet_angle = 25, facet_size =
> 0.5, cell_radius_edge_ratio = 3.5, cell_size = .8, facet_distance = 0.6);
>
> // Create edge that we want to preserve
> Polylines polylines;
> //Polyline_3& polyline = polylines.front();
> //Polyline_3& polyline2 = polylines.back();
>
> Polyline_3 polyline;
> Polyline_3 polyline6;
>
> //top of cylinder
> for (int i = 31; i < 360; ++i)
> {
> Point p(15.0, 10 * std::cos(i * CGAL_PI / 180), 10 * std::sin(i * CGAL_PI
> / 180));
> polyline.push_back(p);
> }
> polylines.push_back(polyline);
>
> //fixateur line: (-15 < x < 15, y = 10 * cos(pi/6), z = 10 * sin(pi/6) )
> for (int i = 0; i < 100; ++i)
> {
> Point p(-15.0 + i * ((30.) / 99.), 10 * std::cos(CGAL_PI / 6), 10 *
> std::sin(CGAL_PI / 6));
> polyline6.push_back(p);
> }
> polylines.push_back(polyline6);
>
> // Insert edge in domain
> domain.add_features(polylines.begin(), polylines.end());
>
> // Mesh generation without feature preservation
> C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria,
> CGAL::parameters::no_features());
>
> // Output
> std::ofstream medit_file("out-no-protection.mesh");
> c3t3.output_to_medit(medit_file);
> medit_file.close();
> c3t3.clear();
>
> // Mesh generation with feature preservation
> c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
>
> // Output
> medit_file.open("out.mesh");
> c3t3.output_to_medit(medit_file);
> medit_file.close();
>
> return 0;
> }
>
>
>
> --
> Sent from: http://cgal-discuss.949826.n4.nabble.com/
>
--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss



Screenshot from 2020-09-08 14-49-18.png (12K) Download Attachment