meshing an ellipsoid

ellipsoid.cpp

#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 <CGAL/IO/output_surface_facets_to_polyhedron.h>
#include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
#include <CGAL/Surface_mesh.h>


// 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 CGAL::Surface_mesh<Point_3> Surface_mesh;


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

class Ellipsoid {
  double aa, bb, cc;

public:
  Ellipsoid(double a, double b, double c)
    : aa(a*a), bb(b*b), cc(c*c)
  {}
  
  double operator() (const Point_3& p) const
  {
    double FT x2=(p.x()*p.x())/(aa), y2=(p.y()*p.y())/(bb), z2=(p.z()*p.z())/(cc);
    return x2+y2+z2-1;
  }
};

int main() {

  Ellipsoid ellipsoid(4,2,1);
  
  Tr tr;            // 3D-Delaunay triangulation
  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation

  // defining the surface
  Surface_3 surface(ellipsoid,             // pointer to function
                    Sphere_3(CGAL::ORIGIN, 20.)); // bounding sphere
  // Note that "20." above is the *squared* radius of the 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());

  Surface_mesh sm;
  
  //  CGAL::output_surface_facets_to_smhedron(c2t3,sm);
  CGAL::facets_in_complex_2_to_triangle_mesh(c2t3,sm);

  typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
 
  BOOST_FOREACH(vertex_descriptor vd, vertices(sm)){
    Point_3 p = sm.point(vd);
    p = Point_3(p.x(), p.y(), p.z()+ 2);  // your transformation
    sm.point(vd) = p;
  }
  std::cout << sm << std::endl;

  return 0;
}