CGAL users discussion list

[Stephane Tayeb <>]

Hi Ianic,

Ianic wrote:
> Dear Stephane,
>
> Thank you for your reply. I mesh a single-domain geometry that has one cell
> family and one facet family, i.e. boundary of my domain. I essentially
> follow the liver example in doing this, and I only change the meshing
> criteria but keep the optimization parameters at their default values. I
> believe it implies that there is a perturbation step at the end. 

You are right, there is a perturbation step (and an exudation one) in 
the default case.

> In my case I actually relied on vertices' tags, because I have to:
>
> - remove unreferenced vertices that are not part of any cell (our finite
> element code won't accept a mesh otherwise). I thought it would be easier to
> NOT output them in CGAL in the first place. For that I added a small if()
> check in the File_medit.h file to write out a vertex only if its tag is not
> -1.
> - I also need to identify what vertices lie on a boundary, i.e. belong to
> facet. My understanding was that if a node belongs to a boundary facet it
> means that, following your terminology, it has dim=2 and therefore it will
> have a facet tag (which is ok since I only have one facet family).
>
> As a reference, my output mesh had appr. 900 000 vertices, of which 430-ish
> had a tag of "-1" and of those only 2 belonged to a c3t3 cell.

I understand your needs. I did more experiments and we thought a bit 
about your problem. We found that, depending on the input data, our 
algorithm could generate isolated points.

I suggest that you use the function of the attached file to clean-up 
your c3t3, then output your mesh as you did. The vertices with -1 index 
in the generated file should be the vertices that are not incident to 
any tetrahedron of the mesh. The vertices with dim=2 (and thus your 
facet tag) should be the vertices which lie on the boundary of the mesh.
You could also directly use the various iterator of the c3t3 object (and 
its embedded triangulation) to write your output.

Best,
Stéphane.

PS: if your data is not confidential, I would be glad if you could send 
it to me so that I can analyze it.

--
Stephane Tayeb
Software engineer - INRIA Sophia Antipolis
Geometrica Project-Team
template <typename C3T3>
void clean_c3t3(C3T3& c3t3)
{
  typedef typename C3T3::Triangulation            Tr;
  typedef typename C3T3::Cell_handle              Cell_handle;
  typedef typename C3T3::Facet                    Facet;
  typedef typename Tr::Finite_vertices_iterator   Fvi;
  
  const Tr& tr = c3t3.triangulation();
  if ( tr.dimension() < 3 ) { return; }
  
  // Iterate on vertices and fix dimension
  for ( Fvi vit = tr.finite_vertices_begin(), end = tr.finite_vertices_end() ;
       vit != end ; ++vit )
  {
    typename std::vector<Facet> facets;
    tr.finite_incident_facets(vit, std::back_inserter(facets));
    CGAL_assertion(!facets.empty());
    
    // Look for the first surface facet
    typename std::vector<Facet>::iterator it_facet = facets.begin();
    while ( it_facet != facets.end() && !c3t3.is_in_complex(*it_facet++) ){}
    
    // vit is on the surface: set dimension and break
    if ( it_facet != facets.end() )
    {
      c3t3.set_dimension(vit,2);
      continue;
    }
    
    // Look for the first interior cell
    typename std::vector<Cell_handle> cells;
    tr.finite_incident_cells(vit,std::back_inserter(cells));
    CGAL_assertion(cells.begin() != cells.end());
    
    if ( c3t3.is_in_complex(*cells.begin()) )
    {
      c3t3.set_dimension(vit,3);
      continue;
    }
    
    // vit is not part of the complex
    c3t3.set_dimension(vit,-1);
  }
}