CGAL users discussion list

[weihy1984 <>]

Hi, I use CGAL generate a ellipsoid surface mesh, the ellipsoid is

x^2+y^2/4+z^2/9=1. But the program print the following error message.  

 

CGAL error: assertion violation!
Expr: false
File: C:\Program Files\CGAL-3.3.1\include\CGAL/Delaunay_triangulation_3.h
Line: 1068

And the line 1068 of Delaunay_triangulation_3.h is

 

 CGAL_triangulation_assertion(false);

 

it always is false if the code arrive this line. Why? Are there any limitations about the 3D surface mesh generation of CGAL?

 

The following is my code:

 

#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <fstream>

// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;

// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;

typedef Tr::Geom_traits GT;
typedef GT::Sphere_3 Sphere_3;
typedef GT::Point_3 Point_3;
typedef GT::FT FT;
typedef Tr::Finite_facets_iterator Fi;

//typedef C2t3::Facet_iterator Fi;
typedef Tr::Finite_vertices_iterator Vi;
typedef C2t3::Vertex_handle Vertex_handle;

typedef FT (*Function)(Point_3);

typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;

 

FT ellipsoid_function (Point_3 p){
  const FT x2=p.x()*p.x(), y2=p.y()*p.y()/4.0, z2=p.z()*p.z()/9.0;
  return x2+y2+z2-1;
}

int main() {

  Tr tr; // 3D-Delaunay triangulation

 C2t3 c2t3(tr); // 2D-complex in 3D-Delaunay triangulation

  // defining the surface
  Surface_3 surface(ellipsoid_function, // pointer to function
                  Sphere_3(CGAL::ORIGIN, 3.5)); // bounding sphere
  // defining meshing criteria
  CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30., // angular bound
                                                   0.1,   // radius bound
                                                   0.1);  // distance bound
  // meshing surface
  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());

   std::ofstream os("out.mesh");
  // outputs dimension and number of vertices
  int n = tr.number_of_vertices();
  int m =c2t3.number_of_facets();

  os << "OFF"<<std::endl<<n <<" "<< m<<" "<<0<<std::endl;

  int i = 0;
  std::map<Vertex_handle, int > V;
  // write the vertices
  for (Vi it=tr.finite_vertices_begin(); it!=tr.finite_vertices_end(); ++it)
  {
    os << it->point() << std::endl;
    V[it] = i++;
  }

  // write the facets
  for(Fi it =tr.finite_facets_begin(); it !=tr.finite_facets_end(); ++it)
  {
    if(c2t3.is_in_complex(*it))
 {
   os<<"3 ";
      for(int j = 0; j < 4; j++)
   {
        if ((*it).second!=j)
     {
          os << V[(*it).first->vertex(j)]<<" ";
        }
     }
     os << std::endl;
    }
  }

  os.close();

return 0;
}

---

穿越地震带 纪念汶川地震一周年