As far as I understand how a facet works in CGAL (using v5.1), it is composed of a cell handle and an index to the vertex opposing the facet in that cell handle. One can thus obtains all the vertex of the considered cell using (I am working using alpha-shapes in 3D):
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT FT;
typedef CGAL::Triangulation_vertex_base_with_info_3<std::size_t, Kernel> Vb3;
typedef CGAL::Fixed_alpha_shape_vertex_base_3<Kernel, Vb3> asVb3;
typedef CGAL::Fixed_alpha_shape_cell_base_3<Kernel> asCb3;
typedef CGAL::Triangulation_data_structure_3<asVb3, asCb3> asTds3;
typedef CGAL::Delaunay_triangulation_3<Kernel, asTds3, CGAL::Fast_location> asTriangulation_3;
typedef CGAL::Fixed_alpha_shape_3<asTriangulation_3> Alpha_shape_3;
typedef Kernel::Point_3 Point_3;
...
Alpha_shape_3::Facet facet;
...
Alpha_shape_3::Vertex_handle facetV1 = facet.first->vertex((facet.second+1)%4);
Alpha_shape_3::Vertex_handle facetV2 = facet.first->vertex((facet.second+2)%4);
Alpha_shape_3::Vertex_handle facetV3 = facet.first->vertex((facet.second+3)%4);
Alpha_shape_3::Vertex_handle facetOppositeVertex = facet.first->vertex((facet.second)%4);
I assumed the opposite vertex is inside a cell, and the facet is also part of this cell. Using this, I can compute the exterior-pointing normals to the surface of my alpha-shape (if there is a way to do this automatically using CGAL, I would like to know it). To compute the normal of a point on the free surface, I do an average of the the different facet normals sharing that point, weighted by the angle of the facet at that point. To orient my normals exteriorly to the mesh, I simply check the dot product between the normal and a vector from one point of the facet to the opposite vertex.
The problem is that, sometimes the opposite vertex seems to be in a complete different cell, leading in my FEM computations to bad results. Here are some figures of what happen:
initially Normal to free surface (good)
when the mesh slightly deforms: Normal to free surface (bad)
Actually what happens is that the opposite vertex of some facet does not lie in the only cell sharing that free surface facet, but inside another cell far away, such that the orientation of the facet normal is not correctly computed, and that some facet normals cancel each other, leading thus to bad normals.
Is there:
- a way to compute those normals automatically in cgal ?
- a problem in my understanding of how a facet works ?
NB: I already check when the opposite vertex is infinite and in that case I mirror the facet using the mirror_facet function
EDIT:
As mentionned in a comment, I tried to use different functions to access the vertex:
std::cout << "--------------------------------------------------" << std::endl;
std::cout << facet.first->vertex(asTriangulation_3::vertex_triple_index(facet.second, 0))->info() << std::endl; //facetV1
std::cout << facet.first->vertex(asTriangulation_3::vertex_triple_index(facet.second, 1))->info() << std::endl; //facetV2
std::cout << facet.first->vertex(asTriangulation_3::vertex_triple_index(facet.second, 2))->info() << std::endl; //facetV3
std::cout << facet.first->vertex(asTriangulation_3::vertex_triple_index(facet.second, 3))->info() << std::endl; //opposite vertex
std::cout << "--------------------------------------------------" << std::endl;
But it sometimes gives me:
--------------------------------------------------
117
25
126
117
--------------------------------------------------
So an opposite vertex with the same index as one of the facet's vertex.
