CGAL users discussion list

[Martin Baeker <>]

Dear all,

I am facing a problem which is probably due to the fact that i do not
sufficiently understand how the template-kernels really work.

The problem is the following:

I am operating on a Nef-Polyhedron and convert this to a Polyhedron.

The next step would be to create a Mesh on this polyhedron, but I
cannot get this to run because it seems that the Mesh_3-stuff and the
Nef-Polyhedron cannot work with the same kernel. (I found that my
Nef-calculation only works with exact_predicates_exact_constructions,
but Mesh_3 seems not to work with this?)

When compiling the program, I get (depending on my compiler version)
either a segmentation fault of gcc or a compiler error.  (Example error
»CGAL::Lazy_exact_nt<CGAL::Gmpq>« kann nicht nach »double« in return
umgewandelt werden -- meaning canno convert
CGAL::Lazy_exact_nt<CGAL::Gmpq> to double in return)

Here is an excerpt from the program (full program below). It is based
on the examples polygon_prog_cut_cube and mesh_polyhedral_domain:

//Typedefs  
typedef CGAL::Exact_predicates_exact_constructions_kernel  Kernel;
typedef CGAL::Mesh_3::Robust_intersection_traits_3<Kernel> Geom_traits;
typedef CGAL::Polyhedron_3<Geom_traits>                           Polyhedron;
typedef CGAL::SNC_indexed_items Items;
typedef CGAL::Nef_polyhedron_3<Kernel, Items>  Nef_polyhedron;
typedef Kernel::Vector_3  Vector_3;
typedef Kernel::Aff_transformation_3  Aff_transformation_3;

typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, Geom_traits> Mesh_domain;

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

// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;

...excerpt from main:

    Polyhedron P;
    Halfedge_handle h = make_cube_3( P);
    Polyhedron Q=P;
       Nef_polyhedron N1(P);
     Nef_polyhedron N2(Q);      
... some operations on N1 and N2
         N2.convert_to_Polyhedron(P);
   Mesh_domain domain(P);
  Facet_criteria facet_criteria(25, 0.15, 0.008); // angle, size, 
approximation
  Cell_criteria cell_criteria(4, 0.2); // radius-edge ratio, size
  Mesh_criteria criteria(facet_criteria, cell_criteria);

  // Mesh generation --- the next command seems to be the troublesome one
     C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);


So how can I create a mesh on a polyhedron that is made froma
Nef_Polyhedron? is it possible? 

Any help would be appreciated,

Martin.

PS: Here is the full program:

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/Nef_3/SNC_indexed_items.h>
#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
#include <CGAL/polyhedron_cut_plane_3.h>
#include <CGAL/IO/Polyhedron_geomview_ostream.h>
#include <iostream>
#include <CGAL/AABB_intersections.h>
#include <CGAL/Mesh_3/Robust_intersection_traits_3.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/refine_mesh_3.h>

  
typedef CGAL::Exact_predicates_exact_constructions_kernel  Kernel;
typedef Kernel::Point_3                                      Point;
typedef Kernel::Plane_3                                      Plane;
typedef CGAL::Mesh_3::Robust_intersection_traits_3<Kernel> Geom_traits;
typedef CGAL::Polyhedron_3<Geom_traits>                           Polyhedron;
typedef Polyhedron::Halfedge_handle                          Halfedge_handle;
typedef Polyhedron::Halfedge_iterator                        
Halfedge_iterator;

typedef CGAL::SNC_indexed_items Items;
typedef CGAL::Nef_polyhedron_3<Kernel, Items>  Nef_polyhedron;
typedef Kernel::Vector_3  Vector_3;
typedef Kernel::Aff_transformation_3  Aff_transformation_3;

typedef CGAL::Polyhedral_mesh_domain_3<Polyhedron, Geom_traits> Mesh_domain;

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

// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;



template <class Poly>
typename Poly::Halfedge_handle make_cube_3( Poly& P) {
    // appends a cube of size [0,1]^3 to the polyhedron P.
    CGAL_precondition( P.is_valid());
    typedef typename Poly::Point_3         Point;
    typedef typename Poly::Plane_3         Plane;
    typedef typename Poly::Halfedge_handle Halfedge_handle;
    Halfedge_handle h = P.make_tetrahedron( Point( 1, 0, 0),
                                            Point( 0, 0, 1),
                                            Point( 0, 0, 0),
                                            Point( 0, 1, 0));
    Halfedge_handle g = h->next()->opposite()->next();
    P.split_edge( h->next());
    P.split_edge( g->next());
    P.split_edge( g);
    h->next()->vertex()->point()     = Point( 1, 0, 1);
    g->next()->vertex()->point()     = Point( 0, 1, 1);
    g->opposite()->vertex()->point() = Point( 1, 1, 0);
    Halfedge_handle f = P.split_facet( g->next(), g->next()->next()->next());
    Halfedge_handle e = P.split_edge( f);
    e->vertex()->point() = Point( 1, 1, 1);
    P.split_facet( e, f->next()->next());
    CGAL_postcondition( P.is_valid());
    g = h;
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    g = h->opposite();
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    g = h->next()->opposite();
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    g = h->next()->next()->opposite();
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    g = h->next()->next()->next()->opposite();
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    g = g->next()->next()->opposite();
    g->facet()->plane() = Plane( g->vertex()->point(),
                                 g->next()->vertex()->point(),
                                 g->next()->next()->vertex()->point());
    return h;
}


int main() {
  CGAL::set_pretty_mode( std::cout);
    Polyhedron P;
    Halfedge_handle h = make_cube_3( P);
    std::cout << "Before" << P;
//     Plane pl = Plane( Point( 0.5, 0.0, 0.0),
//                       Point( 1, 0.0, 1),
//                       Point( 1, 1, 1.));
    // Simple cut in the middle
    
    Plane pl = Plane( Point( 0.5, 0.0, 0.0),
                      Point( 0.5, 0.0, 0.5),
                      Point( 0.5, 0.5, 0.));
    //!!! It seems that the first point of the plane must be on the
    //correct handle of the polygon!!??

    Polyhedron Q=P;
    // Find the corresponding handle to h on Q
    Halfedge_handle hq;
    for (Halfedge_iterator it=Q.halfedges_begin(); it != Q.halfedges_end(); 
++it)
      {
        if ((it->facet()->plane() == h->facet()->plane()) &&
             (it->vertex()->point() == h->vertex()->point()))
          {
            hq = it;
            std::cout << "Hit!\n";
            break;
          }
      }
    std::cout << pl << std::endl;
    CGAL::polyhedron_cut_plane_3( P, h, pl);
    CGAL::polyhedron_cut_plane_3( Q, hq->opposite(), pl.opposite());
    std::cout << P << Q;

       Nef_polyhedron N1(P);
     Nef_polyhedron N2(Q);
     
    Aff_transformation_3 aff(CGAL::TRANSLATION, Vector_3(0,0,0.5));
      N2.transform(aff);
      // Now we would need to join N1 and N2 somehow

      N2=(N1+N2) ;
     
      Polyhedron Res;
     if (N2.is_simple())
       {
         
         N2.convert_to_Polyhedron(Res);
         
         std::cout << "N2 is simple!!!" << Res;
    CGAL::Geomview_stream geo;
    geo << CGAL::GREEN << Res;

    // wait for a mouse click.
    Point click;
    geo >> click;
       }
     else
       std::cout << "N2 is not simple!!" << N2;


         N2.convert_to_Polyhedron(P);
   Mesh_domain domain(P);

//   // Set mesh criteria
  Facet_criteria facet_criteria(25, 0.15, 0.008); // angle, size, 
approximation
  Cell_criteria cell_criteria(4, 0.2); // radius-edge ratio, size
  Mesh_criteria criteria(facet_criteria, cell_criteria);

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

//   // Output
//   std::ofstream medit_file("out.mesh");
//   c3t3.output_to_medit(medit_file);
//   medit_file.close();
     
    return 0;
}



                   Priv.-Doz. Dr. Martin Bäker
                   Institut für Werkstoffe
                   Technische Universität Braunschweig
                   Langer Kamp 8
                   38106 Braunschweig
                   Germany
                   Tel.: 00-49-531-391-3073
                   Fax   00-49-531-391-3058
                   e-mail 
<>