CGAL users discussion list

[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/