void Viewer::triangulate_facet()
{


    pos_facets.resize(0);
    flat_normals.resize(0);
    smooth_normals.resize(0);
    colors.resize(0);

    LCC &lcc = *scene->lcc;

    for (LCC::Attribute_range<3>::type::iterator
         it=lcc.attributes<3>().begin(),
         itend=lcc.attributes<3>().end(); it!=itend; ++it )
    {
        if ( it->info().is_visible() )
        {
            for(LCC::One_dart_per_incident_cell_range<2,3>::iterator
                dartIter=lcc.one_dart_per_incident_cell<2,3>
                (lcc.dart_of_attribute<3>(it)).begin(); dartIter.cont(); ++dartIter)
            {

                //Computes the normal of the facet
                Traits::Vector_3 normal = CGAL::compute_normal_of_cell_2(lcc,dartIter);
                normal = normal/(CGAL::sqrt(normal*normal));

                P_traits cdt_traits(normal);
                CDT cdt(cdt_traits);

               LCC::Dart_of_orbit_range<1>::const_iterator
                 he_circ = lcc.darts_of_orbit<1>(dartIter).begin(),
                 he_circ_end = lcc.darts_of_orbit<1>(dartIter).end();

                // Iterates on the vector of facet handles
                CDT::Vertex_handle previous, first;
                do {
                    CDT::Vertex_handle vh = cdt.insert(lcc.point(he_circ));
                    if(first == 0) {
                        first = vh;
                    }
                    //vh->info() = he_circ;
                    if(previous != 0 && previous != vh) {
                        cdt.insert_constraint(previous, vh);
                    }
                    previous = vh;
                } while( ++he_circ != he_circ_end );
                cdt.insert_constraint(previous, first);

                // sets mark is_external
                for(CDT::All_faces_iterator
                    fit = cdt.all_faces_begin(),
                    end = cdt.all_faces_end();
                    fit != end; ++fit)
                {
                    fit->info().is_external = false;
                }
                //check if the facet is external or internal
                std::queue<CDT::Face_handle> face_queue;
                face_queue.push(cdt.infinite_vertex()->face());
                while(! face_queue.empty() ) {
                    CDT::Face_handle fh = face_queue.front();
                    face_queue.pop();
                    if(fh->info().is_external) continue;
                    fh->info().is_external = true;
                    for(int i = 0; i <3; ++i) {
                        if(!cdt.is_constrained(std::make_pair(fh, i)))
                        {
                            face_queue.push(fh->neighbor(i));
                        }
                    }
                }

                //iterates on the internal faces to add the vertices to the positions
                //and the normals to the appropriate vectors
                for(CDT::Finite_faces_iterator
                    ffit = cdt.finite_faces_begin(),
                    end = cdt.finite_faces_end();
                    ffit != end; ++ffit)
                {
                    if(ffit->info().is_external)
                        continue;

                    //compute normals (no smooth for non-triangle facets objects
                    LCC::Vector normal = CGAL::compute_normal_of_cell_2(lcc,dartIter);
                    normal = normal/(CGAL::sqrt(normal*normal));

                     smooth_normals.push_back(normal.x());smooth_normals.push_back(normal.y());smooth_normals.push_back(normal.z());
                     smooth_normals.push_back(normal.x());smooth_normals.push_back(normal.y());smooth_normals.push_back(normal.z());
                     smooth_normals.push_back(normal.x());smooth_normals.push_back(normal.y());smooth_normals.push_back(normal.z());

                     flat_normals.push_back(normal.x());flat_normals.push_back(normal.y());flat_normals.push_back(normal.z());
                     flat_normals.push_back(normal.x());flat_normals.push_back(normal.y());flat_normals.push_back(normal.z());
                     flat_normals.push_back(normal.x());flat_normals.push_back(normal.y());flat_normals.push_back(normal.z());



                    pos_facets.push_back(ffit->vertex(0)->point().x()); pos_facets.push_back(ffit->vertex(0)->point().y()); pos_facets.push_back(ffit->vertex(0)->point().z());
                    pos_facets.push_back(ffit->vertex(1)->point().x()); pos_facets.push_back(ffit->vertex(1)->point().y()); pos_facets.push_back(ffit->vertex(1)->point().z());
                     pos_facets.push_back(ffit->vertex(2)->point().x()); pos_facets.push_back(ffit->vertex(2)->point().y()); pos_facets.push_back(ffit->vertex(2)->point().z());

//.........这里部分代码省略.........

 void FixedPlaneMesh::TrianglatePolygon(const double3& normal, std::vector<Point3D>& pts, ListOfvertices& results)
 {
    if (pts.size() < 3)
        return ; 
    if (pts.size() ==3)
    {
       VertexInfo vi1(pts[0], normal, mColor);
       VertexInfo vi2(pts[1], normal, mColor);
       VertexInfo vi3(pts[2], normal, mColor);
       results.push_back(vi1);
       results.push_back(vi2);
       results.push_back(vi3);
       return ;
    }
    typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

    typedef CGAL::Triangulation_vertex_base_2<K>                     Vb;
    typedef CGAL::Constrained_triangulation_face_base_2<K>           Fb;
    typedef CGAL::Triangulation_data_structure_2<Vb,Fb>              TDS;
     //typedef CGAL::Exact_predicates_tag                               Itag;
    typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, CGAL::No_intersection_tag> CDT;

    vec3<double> origin = pts[0];
    vec3<double>  N = normal;
    vec3<double>  U = normalize(pts[1] - origin);
    vec3<double>  V = cross(N, U);
    CDT cdt;
    CDT::Vertex_handle vh1, vh2, vh3;
    vec3<double> v0 = PosToLocal(U, V, N, origin, pts[0]);
    CDT::Point p0(v0.x, v0.y);
    vh1 = vh3 = cdt.insert(p0);
    for ( int i = 1; i< pts.size() ; i++)
    {
        vec3<double> v1 = PosToLocal(U, V, N, origin, pts[i]);
        CDT::Point p1(v1.x, v1.y);
        vh2 = cdt.insert(p1);
        cdt.insert_constraint(vh1, vh2);
        vh1 = vh2;
    }
    cdt.insert_constraint(vh2, vh3);
    int count = cdt.number_of_faces() ; 
    results.reserve(count*3);
    for (CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();
       fit != cdt.finite_faces_end(); ++fit)
   {
	   vec2<double> v0(fit->vertex(2)->point().x(),fit->vertex(2)->point().y() );
	   vec2<double> v1(fit->vertex(1)->point().x(),fit->vertex(1)->point().y() );
	   vec2<double> v2(fit->vertex(0)->point().x(),fit->vertex(0)->point().y() );   
	   if (IsEqual(cross(v0- v2, v1-v2), (double)0.,  (double)EPSF ))
		   continue; //
       vec3<double > p0(v0, 0.0);
       vec3<double > p1(v1, 0.0);
       vec3<double > p2(v2, 0.0);
       p0 = PosToGlobal(U, V, N, origin, p0);
       p1 = PosToGlobal(U, V, N, origin, p1);
       p2 = PosToGlobal(U, V, N, origin, p2);
       VertexInfo vi1(p0, N, mColor);
       VertexInfo vi2(p1, N, mColor);
       VertexInfo vi3(p2, N, mColor);
       results.push_back(vi1);
       results.push_back(vi2);
       results.push_back(vi3);
   }
  
 }

void Scene_nef_polyhedron_item::compute_normals_and_vertices(void)
{
     int count = 0;
    positions_facets.resize(0);
    positions_points.resize(0);
    color_lines.resize(0);
    color_facets.resize(0);
    color_points.resize(0);
    normals.resize(0);
    positions_lines.resize(0);
    //The Facets
    {
        for(Nef_polyhedron::Halffacet_const_iterator
            f = nef_poly->halffacets_begin (),
            end = nef_poly->halffacets_end();
            f != end; ++f)
        {
            if(f->is_twin()) continue;
            count++;
            Nef_polyhedron::Vector_3 v = f->plane().orthogonal_vector();
            P_traits cdt_traits(v);
            CDT cdt(cdt_traits);

            for(Nef_polyhedron::Halffacet_cycle_const_iterator
                fc = f->facet_cycles_begin(),
                end = f->facet_cycles_end();
                fc != end; ++fc)
            {
                if ( fc.is_shalfedge() )
                {

                    Nef_polyhedron::SHalfedge_const_handle h = fc;
                    Nef_polyhedron::SHalfedge_around_facet_const_circulator hc(h), he(hc);

                    CDT::Vertex_handle previous, first;

                    do {
                        Nef_polyhedron::SVertex_const_handle v = hc->source();
                        const Nef_polyhedron::Point_3& point = v->source()->point();
                        CDT::Vertex_handle vh = cdt.insert(point);
                        if(first == 0) {
                            first = vh;
                        }
                        vh->info() = hc->source();
                        if(previous != 0 && previous != vh) {
                            cdt.insert_constraint(previous, vh);
                        }
                        previous = vh;
                    } while( ++hc != he );

                    cdt.insert_constraint(previous, first);

                    // sets mark is_external
                    for(CDT::All_faces_iterator
                        fit = cdt.all_faces_begin(),
                        end = cdt.all_faces_end();
                        fit != end; ++fit)
                    {
                        fit->info().is_external = false;

                    }
                    //check if the facet is external or internal
                    std::queue<CDT::Face_handle> face_queue;
                    face_queue.push(cdt.infinite_vertex()->face());

                    while(! face_queue.empty() ) {
                        CDT::Face_handle fh = face_queue.front();
                        face_queue.pop();
                        if(fh->info().is_external) continue;
                        fh->info().is_external = true;
                        for(int i = 0; i <3; ++i) {
                            if(!cdt.is_constrained(std::make_pair(fh, i)))
                            {
                                face_queue.push(fh->neighbor(i));
                            }
                        }

                    }
                    //iterates on the internal faces to add the vertices to the positions
                    //and the normals to the appropriate vectors

                    for(CDT::Finite_faces_iterator
                        ffit = cdt.finite_faces_begin(),
                        end = cdt.finite_faces_end();
                        ffit != end; ++ffit)
                    {


                        if(ffit->info().is_external){ continue;}
                        for(int i = 0; i<3; i++)
                        {
                            positions_facets.push_back(CGAL::to_double(ffit->vertex(i)->point().x()));
                            positions_facets.push_back(CGAL::to_double(ffit->vertex(i)->point().y()));
                            positions_facets.push_back(CGAL::to_double(ffit->vertex(i)->point().z()));

                        }



                        Nef_polyhedron::Vector_3 v = f->plane().orthogonal_vector();
//.........这里部分代码省略.........

void
Scene_polyhedron_item::triangulate_facet(Facet_iterator fit) const
{
    //Computes the normal of the facet
    Traits::Vector_3 normal =
            CGAL::Polygon_mesh_processing::compute_face_normal(fit,*poly);
    //check if normal contains NaN values
    if (normal.x() != normal.x() || normal.y() != normal.y() || normal.z() != normal.z())
    {
        qDebug()<<"Warning : normal is not valid. Facet not displayed";
        return;
    }
    P_traits cdt_traits(normal);
    CDT cdt(cdt_traits);

    Facet::Halfedge_around_facet_circulator
            he_circ = fit->facet_begin(),
            he_circ_end(he_circ);

    // Iterates on the vector of facet handles
    CDT::Vertex_handle previous, first;
    do {
        CDT::Vertex_handle vh = cdt.insert(he_circ->vertex()->point());
        if(first == 0) {
            first = vh;
        }
        vh->info() = he_circ;
        if(previous != 0 && previous != vh) {
            cdt.insert_constraint(previous, vh);
        }
        previous = vh;
    } while( ++he_circ != he_circ_end );
    cdt.insert_constraint(previous, first);
    // sets mark is_external
    for(CDT::All_faces_iterator
        fit2 = cdt.all_faces_begin(),
        end = cdt.all_faces_end();
        fit2 != end; ++fit2)
    {
        fit2->info().is_external = false;
    }
    //check if the facet is external or internal
    std::queue<CDT::Face_handle> face_queue;
    face_queue.push(cdt.infinite_vertex()->face());
    while(! face_queue.empty() ) {
        CDT::Face_handle fh = face_queue.front();
        face_queue.pop();
        if(fh->info().is_external) continue;
        fh->info().is_external = true;
        for(int i = 0; i <3; ++i) {
            if(!cdt.is_constrained(std::make_pair(fh, i)))
            {
                face_queue.push(fh->neighbor(i));
            }
        }
    }
    //iterates on the internal faces to add the vertices to the positions
    //and the normals to the appropriate vectors
    for(CDT::Finite_faces_iterator
        ffit = cdt.finite_faces_begin(),
        end = cdt.finite_faces_end();
        ffit != end; ++ffit)
    {
        if(ffit->info().is_external)
            continue;

        double vertices[3][3];
        vertices[0][0] = ffit->vertex(0)->point().x();
        vertices[0][1] = ffit->vertex(0)->point().y();
        vertices[0][2] = ffit->vertex(0)->point().z();

        vertices[1][0] = ffit->vertex(1)->point().x();
        vertices[1][1] = ffit->vertex(1)->point().y();
        vertices[1][2] = ffit->vertex(1)->point().z();

        vertices[2][0] = ffit->vertex(2)->point().x();
        vertices[2][1] = ffit->vertex(2)->point().y();
        vertices[2][2] = ffit->vertex(2)->point().z();

        positions_facets.push_back( vertices[0][0]);
        positions_facets.push_back( vertices[0][1]);
        positions_facets.push_back( vertices[0][2]);
        positions_facets.push_back(1.0);

        positions_facets.push_back( vertices[1][0]);
        positions_facets.push_back( vertices[1][1]);
        positions_facets.push_back( vertices[1][2]);
        positions_facets.push_back(1.0);

        positions_facets.push_back( vertices[2][0]);
        positions_facets.push_back( vertices[2][1]);
        positions_facets.push_back( vertices[2][2]);
        positions_facets.push_back(1.0);

        typedef Kernel::Vector_3	    Vector;
        Vector n = CGAL::Polygon_mesh_processing::compute_face_normal(fit, *poly);
        normals_flat.push_back(n.x());
        normals_flat.push_back(n.y());
        normals_flat.push_back(n.z());

//.........这里部分代码省略.........

void triangulate(const Polygon_2& polygon, 
		 Cut_iter cuts_begin, Cut_iter cuts_end,
		 const boost::unordered_map<Point_3, boost::unordered_set<Segment_3_undirected> >& point2edges,
		 Out_iter triangles)
{
  typedef CGAL::Triangulation_vertex_base_2<Kernel>                     Vb;
  typedef CGAL::Triangulation_vertex_base_with_info_2<Point_3, Kernel, Vb>     Info;
  typedef CGAL::Constrained_triangulation_face_base_2<Kernel>           Fb;
  typedef CGAL::Triangulation_data_structure_2<Info,Fb>              TDS;
  typedef CGAL::Exact_predicates_tag                               Itag;
  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel, TDS, Itag> CDT;
  typedef CDT::Vertex_handle Vertex_handle;

  static log4cplus::Logger logger = log4cplus::Logger::getInstance("polygon_utils");

  Polygon_2 p(polygon);
  LOG4CPLUS_TRACE(logger, "Triangulating " << pp(p));
  if (p.size() < 3) return;

  bool vertical = is_vertical(p);
  if (vertical)
  {
    LOG4CPLUS_TRACE(logger, "Polygon is vertical.  Rotating.");
    p = yz_swap_neg(p);
  }

  bool reverse = !p.is_counterclockwise_oriented();
  if (reverse)
    p.reverse_orientation();

  CDT cdt;

  boost::unordered_map<Point_3, Vertex_handle> point2handle;
  for (Polygon_2::Vertex_iterator it = p.vertices_begin(); it != p.vertices_end(); ++it)
  {
    Vertex_handle h = cdt.insert(*it);
    point2handle[*it] = h;
    h->info() = *it;//it->z();
  }

  Polygon_2::Vertex_circulator start = p.vertices_circulator();
  Polygon_2::Vertex_circulator c = start;
  Polygon_2::Vertex_circulator n = c;
  ++n;
  do
  {
    Vertex_handle ch = point2handle[*c];//cdt.insert(*c);
    Vertex_handle nh = point2handle[*n];//cdt.insert(*n);
//     ch->info() = c->z();
//     nh->info() = n->z();
//     cdt.insert_constraint(*c, *n);
    cdt.insert_constraint(ch, nh);
    ++c;
    ++n;
  } while (c != start);

  for (Cut_iter c_it = cuts_begin; c_it != cuts_end; ++c_it)
  {
    Polyline_2 cut = *c_it;
    LOG4CPLUS_TRACE(logger, "Adding cut: " << pp(cut));
    if (vertical)
      cut = yz_swap_neg(cut);
    for (Polyline_2::const_iterator c = cut.begin(); c != cut.end(); ++c)
    {
      Polyline_2::const_iterator n = c;
      ++n;
      if (n != cut.end())
      {
	const Point_3& cp = *c;
	const Point_3& np = *n;
	if (point2handle.find(cp) == point2handle.end())
	{
	  Vertex_handle h = cdt.insert(cp);
	  point2handle[cp] = h;
	  h->info() = cp;//cp.z();
	}
	if (point2handle.find(np) == point2handle.end())
	{
	  Vertex_handle h = cdt.insert(np);
	  point2handle[np] = h;
	  h->info() = np;//np.z();
	}

	Vertex_handle ch = point2handle[*c];//cdt.insert(*c);
	Vertex_handle nh = point2handle[*n];//cdt.insert(*n);
// 	ch->info() = c->z();
// 	nh->info() = n->z();
// 	cdt.insert_constraint(*c, *n);
	cdt.insert_constraint(ch, nh);
	LOG4CPLUS_TRACE(logger, "  " << pp(Segment_2(*c, *n)));
      }
    }
  }

  // Loop through the triangulation and store the vertices of each triangle
  for (CDT::Finite_faces_iterator ffi = cdt.finite_faces_begin();
       ffi != cdt.finite_faces_end();
       ++ffi)
  {
    Triangle t;
//.........这里部分代码省略.........

void Viewer::compute_face(Dart_handle dh, LCC::size_type markface)
{
  LCC &lcc = *scene->lcc;

  CGAL::mark_cell<LCC, 2>(lcc, dh, markface);

  double r = (double)lcc.info<3>(dh).color().r()/255.0;
  double g = (double)lcc.info<3>(dh).color().g()/255.0;
  double b = (double)lcc.info<3>(dh).color().b()/255.0;
  if ( !lcc.is_free(dh, 3) )
  {
    r += (double)lcc.info<3>(lcc.beta(dh,3)).color().r()/255.0;
    g += (double)lcc.info<3>(lcc.beta(dh,3)).color().g()/255.0;
    b += (double)lcc.info<3>(lcc.beta(dh,3)).color().b()/255.0;
    r /= 2; g /= 2; b /= 2;
  }

  //compute flat normals
  LCC::Vector normal = CGAL::compute_normal_of_cell_2(lcc,dh);
  normal = normal/(CGAL::sqrt(normal*normal));

  if (lcc.beta<1,1,1>(dh)!=dh)
  {
    P_traits cdt_traits(normal);
    CDT cdt(cdt_traits);

    // Iterates on the vector of facet handles
    CDT::Vertex_handle previous = NULL, first = NULL;
    for (LCC::Dart_of_orbit_range<1>::const_iterator
           he_circ = lcc.darts_of_orbit<1>(dh).begin(),
           he_circ_end = lcc.darts_of_orbit<1>(dh).end();
         he_circ!=he_circ_end; ++he_circ)
    {
      CDT::Vertex_handle vh = cdt.insert(lcc.point(he_circ));
      if(first == NULL)
      { first = vh; }
      vh->info().v = CGAL::compute_normal_of_cell_0<LCC>(lcc, he_circ);
      if(previous!=NULL && previous != vh)
      { cdt.insert_constraint(previous, vh); }
      previous = vh;
    }
    if (previous!=NULL)
      cdt.insert_constraint(previous, first);

    // sets mark is_external
    for(CDT::All_faces_iterator fit = cdt.all_faces_begin(),
          fitend = cdt.all_faces_end(); fit!=fitend; ++fit)
    {
      fit->info().is_external = true;
      fit->info().is_process = false;
    }
    //check if the facet is external or internal
    std::queue<CDT::Face_handle> face_queue;
    CDT::Face_handle face_internal = NULL;
    face_queue.push(cdt.infinite_vertex()->face());
    while(! face_queue.empty() )
    {
      CDT::Face_handle fh = face_queue.front();
      face_queue.pop();
      if(!fh->info().is_process)
      {
        fh->info().is_process = true;
        for(int i = 0; i <3; ++i)
        {
          if(!cdt.is_constrained(std::make_pair(fh, i)))
          {
            face_queue.push(fh->neighbor(i));
          }
          else if (face_internal==NULL)
          {
            face_internal = fh->neighbor(i);
          }
        }
      }
    }

    if ( face_internal!=NULL )
      face_queue.push(face_internal);

    while(! face_queue.empty() )
    {
      CDT::Face_handle fh = face_queue.front();
      face_queue.pop();
      if(!fh->info().is_process)
      {
        fh->info().is_process = true;
        fh->info().is_external = false;
        for(int i = 0; i <3; ++i)
        {
          if(!cdt.is_constrained(std::make_pair(fh, i)))
          {
            face_queue.push(fh->neighbor(i));
          }
        }
      }
    }

    //iterates on the internal faces to add the vertices to the positions
    //and the normals to the appropriate vectors
    for(CDT::Finite_faces_iterator ffit = cdt.finite_faces_begin(),
//.........这里部分代码省略.........

int main() {
	//construct two non-intersecting nested polygons
	Polygon_2 polygon1;
	polygon1.push_back(Point_2(0.0, 0.0));
	polygon1.push_back(Point_2(2.0, 0.0));
	polygon1.push_back(Point_2(1.7, 1.0));
	polygon1.push_back(Point_2(2.0, 2.0));
	polygon1.push_back(Point_2(0.0, 2.0));
	Polygon_2 polygon2;
	polygon2.push_back(Point_2(0.5, 0.5));
	polygon2.push_back(Point_2(1.5, 0.5));
	polygon2.push_back(Point_2(1.5, 1.5));
	polygon2.push_back(Point_2(0.5, 1.5));

	//Insert the polyons into a constrained triangulation
	CDT cdt;
	insert_polygon(cdt, polygon1);
	insert_polygon(cdt, polygon2);

	//Extract point and provide the an index
	std::vector< triangulation_point > points ;
	for ( CDT::Vertex_iterator it = cdt.vertices_begin(); it != cdt.vertices_end(); ++it ){
		it->info() = points.size() ;
		points.push_back( it->point() );
	}


	//Mark facets that are inside the domain bounded by the polygon
	mark_domains(cdt);

	//
	int count = 0;
	for (CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(); fit != cdt.finite_faces_end(); ++fit) {
		if (fit->info().in_domain()){
			++count;
		}
	}


	/*
	 * export
	 */

	std::ofstream ofs("polygon_triangulation2.obj");
	if ( ! ofs.good() ){
		std::cout << "can't open file" << std::endl;
		return 1 ;
	}

	//-- print vertices
	ofs << "# " << points.size() << " vertices"<< std::endl ;
	for ( size_t i = 0; i < points.size(); i++ ){
		ofs << "v " << points[i] << " 0.0" << std::endl;
	}

	//-- print faces
	ofs << "# " << cdt.number_of_faces() << " faces"<< std::endl ;
	// warning : Delaunay_triangulation_2::All_faces_iterator iterator over infinite faces
	for ( CDT::Finite_faces_iterator it = cdt.finite_faces_begin(); it != cdt.finite_faces_end(); ++it )
	{
		//ignore holes
		if ( ! it->info().in_domain() ){
			continue ;
		}
		size_t ia = it->vertex(0)->info();
		size_t ib = it->vertex(1)->info();
		size_t ic = it->vertex(2)->info();

		assert( it->is_valid() );
		//assert ( ia < cdt.number_of_vertices() || ib < tri.number_of_vertices() || ic < tri.number_of_vertices() ) ;

		ofs << "f " << ( ia + 1 ) << " " << ( ib + 1 ) << " " << ( ic + 1 ) << std::endl;
	}

	return 0;
}

//.........这里部分代码省略.........
    
    vertex_range_type vertices = viennagrid::elements<viennagrid::vertex_tag>(plc);
    for (vertex_range_iterator it = vertices.begin(); it != vertices.end(); ++it)
    {
        vertex_handle_type const & vtx_handle = it.handle();
        vertex_type const & vtx = *it;
        point_type const & vgrid_point = viennagrid::point( mesh, vtx );
        
        Vertex_handle handle = cdt.insert( Point(vgrid_point[0], vgrid_point[1]) );
        
        vertex_handle_map[vtx_handle] = handle;
    }

    
    line_range_type lines = viennagrid::elements<viennagrid::line_tag>(plc);
    for (line_range_iterator it = lines.begin(); it != lines.end(); ++it)
    {
        line_type & line = *it;
        
        vertex_handle_type vgrid_v0 = viennagrid::elements<viennagrid::vertex_tag>(line).handle_at(0);
        vertex_handle_type vgrid_v1 = viennagrid::elements<viennagrid::vertex_tag>(line).handle_at(1);
        
        Vertex_handle cgal_v0 = vertex_handle_map[vgrid_v0];
        Vertex_handle cgal_v1 = vertex_handle_map[vgrid_v1];
        
        cdt.insert_constraint(cgal_v0, cgal_v1);
    }
    
    std::vector<point_type> & vgrid_list_of_holes = viennagrid::hole_points(plc);
    std::list<Point> cgal_list_of_holes;
    
    for (std::vector<point_type>::iterator it = vgrid_list_of_holes.begin(); it != vgrid_list_of_holes.end(); ++it)
        cgal_list_of_holes.push_back( Point( (*it)[0], (*it)[1] ) );
    
    CGAL::refine_Delaunay_mesh_2(cdt, cgal_list_of_holes.begin(), cgal_list_of_holes.end(), Criteria());
    
    std::cout << "Number of vertices: " << cdt.number_of_vertices() << std::endl;
    std::cout << "Number of finite faces: " << cdt.number_of_faces() << std::endl;
    
    
    
    
    
    typedef viennagrid::triangular_2d_mesh triangle_mesh_type;
    triangle_mesh_type triangle_mesh;
    
    typedef viennagrid::result_of::point<triangle_mesh_type>::type triangle_point_type;
     
    typedef viennagrid::result_of::element<triangle_mesh_type, viennagrid::vertex_tag>::type triangle_vertex_type;
    typedef viennagrid::result_of::handle<triangle_mesh_type, viennagrid::vertex_tag>::type triangle_vertex_handle_type;
    
    typedef viennagrid::result_of::element<triangle_mesh_type, viennagrid::line_tag>::type triangle_line_type;
    typedef viennagrid::result_of::handle<triangle_mesh_type, viennagrid::line_tag>::type triangle_line_handle_type;
    
    typedef viennagrid::result_of::element<triangle_mesh_type, viennagrid::triangle_tag>::type triangle_triangle_type;
    typedef viennagrid::result_of::handle<triangle_mesh_type, viennagrid::triangle_tag>::type triangle_triangle__handle_type;
    
    
    std::map<Point, triangle_vertex_handle_type> points;
    
    int mesh_faces_counter = 0;
    for(CDT::Finite_faces_iterator fit = cdt.finite_faces_begin(); fit != cdt.finite_faces_end(); ++fit) 
    {
        if(fit->is_in_domain())
        {
            typedef CDT::Triangle Triangle;
            Triangle tri = cdt.triangle(fit);
            
            triangle_vertex_handle_type vgrid_vtx[3];
            
            for (int i = 0; i < 3; ++i)
            {
                std::map<Point, triangle_vertex_handle_type>::iterator pit = points.find( tri[i] );
                if (pit == points.end())
                {
                    vgrid_vtx[i] = viennagrid::make_vertex( triangle_mesh, triangle_point_type(tri[i].x(), tri[i].y()) );
                    points[ tri[i] ] = vgrid_vtx[i];
                }
                else
                    vgrid_vtx[i] = pit->second;
            }
            
            viennagrid::make_element<triangle_triangle_type>( triangle_mesh, vgrid_vtx, vgrid_vtx+3 );
                
            
            std::cout << tri << std::endl;
            ++mesh_faces_counter;
        }
    }
    std::cout << "Number of faces in the mesh mesh: " << mesh_faces_counter << std::endl;
    

    
    std::copy( viennagrid::elements<triangle_triangle_type>(triangle_mesh).begin(), viennagrid::elements<triangle_triangle_type>(triangle_mesh).end(), std::ostream_iterator<triangle_triangle_type>(std::cout, "\n") );
    

    viennagrid::io::vtk_writer<triangle_mesh_type> vtk_writer;
    vtk_writer(triangle_mesh, "test");
    
}