CGAL users discussion list

[Stephane Tayeb <>]

Ianic wrote:
> I am still in the process of getting the neighbors list for all cells and
> facets of a 3D mesh. Here is the code that I've written (I simply added it
> to an example with a liver image):
>
>
> ...
> typedef Tr::Cell_handle Cell_handle;
> typedef C3t3::Cell_iterator Cell_iterator;
> typedef Tr::Facet Facet;
> typedef C3t3::Facet_iterator Facet_iterator;
> typedef Tr::Facet_circulator Facet_circulator;
> typedef Tr::Edge Edge;
> typedef Tr::Vertex_handle Vertex_handle;
> ...
>
>   //---------------------------
>   // A mesh connectivity file
>   //---------------------------
>   std::ofstream neigh_file;
>   neigh_file.open("liver.neigh");
>
>   //---------------------------
>   // 1.0 Cell connectivity
>   //---------------------------
>
>   // Enumerating cells
>   std::map<Cell_handle, int> C;
>   int icell = 1;
>   for( Cell_iterator cit = c3t3.cells_begin() ;
>        cit != c3t3.cells_end() ;
>        ++cit )
>     C[cit] = icell++;
>
>   // Writing cell connectivity data
>   neigh_file << "Cell neighbours" << std::endl
>      << c3t3.number_of_cells() << std::endl;
>
>   for( Cell_iterator cit = c3t3.cells_begin() ;
>        cit != c3t3.cells_end() ;
>        ++cit ) {
>     const Cell_handle& c = cit;
>     for ( int i=0 ; i<4 ; ++i ) {
>       const Cell_handle& neighbor = c->neighbor(i);
>       if ( c3t3.is_in_complex(neighbor) )
>         neigh_file << C[neighbor] << " ";
>       else
>         neigh_file << -1 << " ";
>     }
>     neigh_file << std::endl;
>   }
>
>   //---------------------------
>   // 2.0 Facet connectivity
>   //---------------------------
>
>   // Enumerating facets
>   std::map<Facet, int> F;
>   int ifacet = 1;
>   for( Facet_iterator fit = c3t3.facets_begin() ;
>        fit != c3t3.facets_end() ;
>        ++fit )
>     F[*fit] = ifacet++;
>
>   // Writing facet connectivity data
>   neigh_file << "Facet neighbours" << std::endl
>      << c3t3.number_of_facets() << std::endl;
>
>   std::vector<int> vert;
>   vert.resize(3);
>   std::vector<Edge> edge;
>   edge.resize(3);
>
>   for( Facet_iterator fit = c3t3.facets_begin();
>        fit != c3t3.facets_end();
>        ++fit ) {
>   
>     const Cell_handle& cell = fit->first;
>     int indexInCell = fit->second;
>     vert[0] = (indexInCell+1)%4;
>     vert[1] = (indexInCell+2)%4;
>     vert[2] = (indexInCell+3)%4;
>
>     edge[0] = Edge(cell, vert[1], vert[2]);
>     edge[1] = Edge(cell, vert[2], vert[0]);
>     edge[2] = Edge(cell, vert[0], vert[1]);
>
>     for (unsigned int i = 0; i < edge.size(); ++i) {
>       Facet_circulator fc = c3t3.triangulation().incident_facets(edge[i],
> *fit);
>       Facet_circulator done(fc);
>       do {
>         if ( c3t3.is_in_complex(*fc) )
>           neigh_file << F[*fc] << " ";
>       } while (++fc != done);
>     }
>     neigh_file << std::endl;
>   }
>   neigh_file << "End" << std::endl;
>   neigh_file.close();
>
> The cell neighbors list seems to be correct, though I am still wondering
> what the neighbor(i) holds for a boundary cell, a NULL pointer?

Hi,

No. It is a Cell_handle pointing to:
 + an infinite cell, i.e. a cell which is incident to the infinite 
vertex (you could have a look in 3D Triangulation package documentation, 
and also here: 
http://www.cgal.org/Manual/last/doc_html/cgal_manual/TriangulationDS_3/Chapter_main.html).

 + a finite cell which is not part of the domain discretization

> The facet code doesn't work. The idea essentially is that I take all 3 edges
> of each facet, circulate over their incident facets that belong to a c3t3
> complex (should be 2) and for each edge I write out the facet number which
> is not the current one (this check is omitted in this code just now).
> However, when I run this code, I get lots of zeros from this line:
>
> neigh_file << F[*fc] << " ";
>
> Why does it happen? Can it be because I assign IDs to facets using typedef
> C3t3::Facet_iterator, while when I go over the edge facets, I use typedef
> Tr::Facet_circulator?

I didn't analyze deeply your code, but maybe this could help you.

In CGAL 3D Triangulation, Facet type is std::pair<Cell_handle,int>. 
Thus, each facet of the triangulation can be seen from two side (the 
side of each cell containing the facet). You could have a look at the 
method mirror_facet(Facet) of the triangulation API.
In your code, it seems that you don't take this into account. For 
example, you could fill your map with the facet seen from the cell which 
has the smallest Cell_handle, and do the same when you retrieve the data 
in your loop.

Best,
Stéphane.

--
Stephane Tayeb
Software engineer - INRIA Sophia Antipolis
Geometrica Project-Team