CGAL users discussion list

[Maarten Moesen <>]

Dear all,

I'm currently trying out the sharp-feature
preservation of CGAL's 3D mesh generator on a 32-bit machine with MSVC
2008 (version 3.8, with Laurent Rineau's Triangulation_data_structure_3.h
patch applied). My test geometry is a box having a big sphere cut out from
its center. I add  as 1D sharp features  using line segments
the circles where the sphere cuts the sides of the box, as well as the
ribs of the box. Unfortunately, the mesher does not preserve all the features
near the circle... most of the times a few segments are missing. The sample
code is below. An example mesh is also attached.

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

#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>

namespace POMO
{
const double PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798;
}

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel
K;
typedef K::FT FT;
typedef K::Point_3 Point;

class Cubic_function : public std::unary_function<Point,FT>

{
public:
  Cubic_function(double const&
radius = 0.6): m_x0(0), m_y0(0), m_z0(0), m_x1(1), m_y1(1), m_z1(1), m_xc(0.5),
m_yc(0.5), m_zc(0.5), m_r(radius) {}

  FT const& x_min() const {
return m_x0; }
  FT const& x_max() const {
return m_x1; }
  FT const& y_min() const {
return m_y0; }
  FT const& y_max() const {
return m_y1; }
  FT const& z_min() const {
return m_z0; }
  FT const& z_max() const {
return m_z1; }
  FT const& x_cen() const {
return m_xc; }
  FT const& y_cen() const {
return m_yc; }
  FT const& z_cen() const {
return m_zc; }

  FT const& radius() const
{ return m_r; }
  
  FT square(FT const& x) const
  {
    return x*x;
  }

  FT operator()(const Point&
p) const
  {
    return std::max(std::max(std::max(std::max(x_min()-p.x(),
p.x()-x_max()), std::max(y_min()-p.y(), p.y()-y_max())), std::max(z_min()-p.z(),
p.z()-z_max())),
           
        square(radius()) - square(x_cen()-p.x()) -
square(y_cen()-p.y()) - square(z_cen()-p.z()) ); 
  }

  K::Sphere_3 bounding_sphere()
const
  {
    return K::Sphere_3(Point(0.04,
0.04, 0.04), 3.0);
    
  }

private:
  FT m_x0, m_x1, m_y0, m_y1, m_z0,
m_z1, m_xc, m_yc, m_zc, m_r;
};

//typedef FT (Function)(const Point&);
typedef CGAL::Implicit_mesh_domain_3<Cubic_function,K>
Implicit_mesh_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Implicit_mesh_domain>
Mesh_domain;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type
Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index>
C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr>
Mesh_criteria;

typedef std::vector<Point>  
     Polyline;
typedef std::vector<Polyline>
    Polylines;

int main()
{
  CGAL::Timer timer;

  Cubic_function fun;
  Mesh_domain domain(fun, fun.bounding_sphere(),
1e-6);

  double edge_length = 0.05;
  // Mesh criteria
  Mesh_criteria criteria(CGAL::parameters::facet_angle=30,
    CGAL::parameters::edge_size=edge_length,
    CGAL::parameters::facet_size=edge_length*2,

    CGAL::parameters::facet_distance=0.001,
    CGAL::parameters::cell_radius_edge_ratio=2,
    CGAL::parameters::cell_size=edge_length*2);
  

  // The ribs of the cube
  Point p000(fun.x_min(), fun.y_min(),
fun.z_min());
  Point p001(fun.x_min(), fun.y_min(),
fun.z_max());  
  Point p010(fun.x_min(), fun.y_max(),
fun.z_min());
  Point p011(fun.x_min(), fun.y_max(),
fun.z_max());
  Point p100(fun.x_max(), fun.y_min(),
fun.z_min());
  Point p101(fun.x_max(), fun.y_min(),
fun.z_max());  
  Point p110(fun.x_max(), fun.y_max(),
fun.z_min());
  Point p111(fun.x_max(), fun.y_max(),
fun.z_max());
  
  Polylines polylines(18);
  polylines[0].push_back(p000);
polylines[0].push_back(p100);
  polylines[1].push_back(p001);
polylines[1].push_back(p101);
  polylines[2].push_back(p010);
polylines[2].push_back(p110);
  polylines[3].push_back(p011);
polylines[3].push_back(p111);
  polylines[4].push_back(p000);
polylines[4].push_back(p010);
  polylines[5].push_back(p001);
polylines[5].push_back(p011);
  polylines[6].push_back(p100);
polylines[6].push_back(p110);
  polylines[7].push_back(p101);
polylines[7].push_back(p111);
  polylines[8].push_back(p000);
polylines[8].push_back(p001);
  polylines[9].push_back(p010);
polylines[9].push_back(p011);
  polylines[10].push_back(p100);
polylines[10].push_back(p101);
  polylines[11].push_back(p110);
polylines[11].push_back(p111);

  double rad = CGAL::to_double(fun.radius());
  rad = std::sqrt(rad*rad - 0.25);

  int n = std::max(8, static_cast<int>(std::ceil(2*POMO::PI*rad/edge_length)));

  double* cs = new double[n+1];
  double* ss = new double[n+1];
  for (int ci = 0; ci < n; ++ci)
  {
    double angle = 2*POMO::PI*ci/n;
    cs[ci] = 0.5 + rad*std::cos(angle);
    ss[ci] = 0.5 + rad*std::sin(angle);
  }
  cs[n] = cs[0];
  ss[n] = ss[0];

  std::cout << "Radius:
" << rad << "\nNumber of points: " <<
n << std::endl;

  for (int ci = 0; ci <= n;
++ci)
  {
    //std::cout << fun(Point(0.0,
cs[ci], ss[ci])) << "   " << fun(Point(1.0,
cs[ci], ss[ci])) << std::endl;
    polylines[12].push_back(Point(0.0,
cs[ci], ss[ci]));
    polylines[13].push_back(Point(1.0,
cs[ci], ss[ci]));
    polylines[14].push_back(Point(cs[ci],
0.0, ss[ci]));
    polylines[15].push_back(Point(cs[ci],
1.0, ss[ci]));
    polylines[16].push_back(Point(cs[ci],
ss[ci], 0.0));
    polylines[17].push_back(Point(cs[ci],
ss[ci], 1.0));
  }
  
  delete [] cs;
  delete [] ss;

  domain.add_features(polylines.begin(),
polylines.end());
  
  // Mesh generation
  std::cout << "Generating
the mesh: ";
  timer.reset();timer.start();
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
criteria, CGAL::parameters::no_exude(), CGAL::parameters::no_perturb());
  timer.stop(); 
  std::cout << timer.time()
<< " s." << std::endl;

  /*
  std::cout << "ODT
Smoothing: ";
  timer.reset();timer.start();
  CGAL::Mesh_optimization_return_code
code = CGAL::odt_optimize_mesh_3(c3t3, domain);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;

  std::cout << "Perturbation:
";
  timer.reset();timer.start();
  code = CGAL::perturb_mesh_3(c3t3,
domain);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;

  std::cout << "Sliver
exudation: ";
  timer.reset();timer.start();
  code = CGAL::exude_mesh_3(c3t3);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;
  */
  // Output
  std::ofstream medit_file("sphere_in_box.mesh");
  c3t3.output_to_medit(medit_file);
  medit_file.close();
  std::cout << "Mesh
saved to sphere_in_box.mesh" << std::endl;

  std::cout << "Number
of cells in complex: " << c3t3.number_of_cells() << std::endl;
  std::cout << "Number
of cells in triangulation: " << c3t3.triangulation().number_of_cells()
<< std::endl;
  std::cout << "Number
of finite cells in triangulation: " << c3t3.triangulation().number_of_finite_cells()
<< std::endl;
  std::cout << "Number
of vertices in triangulation: " << c3t3.triangulation().number_of_vertices()
<< std::endl;
  
  return 0;
}

--
Maarten Moesen, PHD
Polymer Physicist Computational Modeling

HUNTSMAN (Europe) BVBA
Everslaan 45
B-3078 Everberg

Office Phone: 0032 (0) 2 758 9962