CGAL users discussion list

[<>]

Hello, 
trying to get the connectivity of a 2d triangulation which is created out of a
set of given data points, I wrote the code which is posted below. It is
actually partly a reimplementation of the header file
'CGAL/examples/Polyhedron_IO/triangulation_print_OFF.h'.  The program reads 
the
coordinates of the nodes, written in a file 'data.cin' just like in the
triangulation examples, and gives back the indexed nodes and the element
connectivities. In my own implementation the mesh is generated by using the 
library DOLFIN, and the specific code lines are commented out.  


I would be happy to have any feedback regarding the code, its accuracy and of
course any improvement, or different implementation (maybe using the function
id() just like in the example 'dijkstra_with_internal_properties.cpp' ??)

#include <fstream>
#include <iostream>
#include <map>

//#include <dolfin.h>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_2.h>

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_id_2.h>

#include <CGAL/IO/File_writer_OFF.h>
#include "triangulation_print_OFF.h"

struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};

typedef CGAL::Triangulation_vertex_base_with_id_2<K> Tvb;
typedef CGAL::Triangulation_face_base_2<K> Tfb;
typedef CGAL::Triangulation_data_structure_2<Tvb,Tfb> Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds> Triangulation;

typedef K::Point_2 Point;
typedef Triangulation::Vertex           Vertex;
typedef Triangulation::Vertex_iterator  Vertex_iterator;
typedef Triangulation::Face_iterator    Face_iterator;

int main()
{
  std::ifstream in("./data2d.cin");
  std::istream_iterator<Point> begin(in);
  std::istream_iterator<Point> end;

  Triangulation t;
  t.insert(begin, end);
  uint number_of_vertices = t.number_of_vertices();
  uint number_of_elements = t.number_of_faces();
  std::cout <<"nr of nodes   =  "<< number_of_vertices <<std::endl;
  std::cout <<"nr of elements=  "<< number_of_elements <<std::endl;
  
  //to test the mesh, output a .off file
  std::ofstream outfile("mesh.off");
  triangulation_print_OFF(outfile, t);

  // Create Mesh
  //dolfin::Mesh mesh;

  // Create a MeshEditor (using dolfin)
  //dolfin::MeshEditor mesh_editor;

  // Initialise mesh editor (using dolfin)
  //mesh_editor.open(mesh,"triangle", 2, 2);
 
  // Specify number of vertices (using dolfin)
  //mesh_editor.init_vertices(number_of_vertices);
  
  // Specify number of elements (using dolfin)
  //mesh_editor.init_cells(number_of_elements);
  
  // Build a map from vertex pointers to vertex indices,
  // and create mesh vertices. (using dolfin)
  std::map<const Vertex*,std::size_t, std::less<const Vertex*> > mapping;
  std::size_t vn = 0;
       for (Vertex_iterator vertex = t.vertices_begin(); vertex !=
t.vertices_end(); ++vertex)
  {
    CGAL_assertion(!t.is_infinite(vertex));
    mapping[&*vertex] = vn;

    const double x = vertex->point().x();
    const double y = vertex->point().y();
    // add vertex to the mesh (using dolfin)
    //mesh_editor.add_vertex(vn,x, y);

    std::cout <<"x and y coordiantes of the node " << vn <<std::endl;
    std::cout << x << std::endl;
    std::cout << y << std::endl;

    vn++;
  }  

   // calculate the connectivity of the faces
   // and add data to the mesh
  std::size_t fn = 0;
  Face_iterator fi = t.faces_begin();
  for ( ; fi != t.faces_end(); ++fi)
  {
    CGAL_assertion( ! t.is_infinite( fi));

    //std::cout << *(fi->vertex(1)) << std::endl;
    CGAL_assertion( mapping.find(&*(fi->vertex(0))) != mapping.end());
    CGAL_assertion( mapping.find(&*(fi->vertex(1))) != mapping.end());
    CGAL_assertion( mapping.find(&*(fi->vertex(2))) != mapping.end());

    uint cell_index = fn;
    //unsigned int v0 = static_cast<unsigned int>(*(fi->vertex(0)));
    //uint v0 = mapping[&*(fi->vertex(0))];
    //uint v1 = mapping[&*(fi->vertex(1))];
    //uint v2 = mapping[&*(fi->vertex(2))];
    // add triangular element (using dolfin)
    //mesh_editor.add_cell(cell_index, v0, v1, v2);
    
    std::cout << "connectivity of the element "<< fn << std::endl;
    std::cout << mapping[&*(fi->vertex(0))] << std::endl;
    std::cout << mapping[&*(fi->vertex(1))] << std::endl;
    std::cout << mapping[&*(fi->vertex(2))] << std::endl;
     
        
   fn++;
  }

  // Close mesh editor
  //mesh_editor.close();
  //plot(mesh);

    return 0;
}