CGAL users discussion list

[Benjamin Kehlet <>]

Hello CGAL people!

I got a report of some unexpected results in mshr (a mesh generation
component in the FEniCS project which uses CGAL heavily). I have
simplified it a bit and reproduced the issue as a standalone CGAL
program.

When meshing a box with a crude triangulation of a sphere subtracted,
the sphere is sometimes "missing" in the resulting mesh.

When meshing the polyhedron in smaller.off (which is basically an axis
aligned box with corners at (-1, -1, -.4) and (1,1,1) with a sphere
with center at origin and radius .3 subtracted) it works fine. See the
surface plot of the resulting mesh in smaller.png.

However, with bigger.off (where the first corner is at (-1,-1,-.7) and
the rest is equal) the inner sphere can not be seen in the plot. See
bigger.png.

I confirmed that it is not a plotting issue by calculating the volume
of the mesh. For bigger.off the result is 6.8 which is equal to the
volume of the box (2*2*1.7). For smaller.off it calculates the volume
to 5.5223. I assume is the correct solution since it is smaller than
5.6 which is the volume of the box (2*2*1.4=5.6). The program source
code is attached. The issue seems to be insensitive to the meshing
parameters.

Best regards

Benjamin Kehlet
#define CGAL_MESH_3_VERBOSE
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>

#include <cstdlib>

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
typedef CGAL::Mesh_polyhedron_3<K>::type MeshPolyhedron_3;
typedef K::Vector_3 Vector_3;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index> C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

typedef K::Tetrahedron_3 Tetrahedron_3;

// To avoid verbose function and named parameters call
using namespace CGAL::parameters;

int main(int argc, char** argv)
{
  MeshPolyhedron_3 P;
  {
    std::cout << "Reading form file: " << argv[1] << std::endl;
    std::ifstream infile(argv[1]);
    infile >> P;
  }

  // Create domain
  Mesh_domain domain(P);

  // Get sharp features
  domain.detect_features();

  const double cs = 0.22;

  // Mesh criteria
  Mesh_criteria criteria(edge_size = cs,
                         facet_angle = 30, 
                         facet_size = cs, 
                         facet_distance = cs/10.0,
                         cell_radius_edge_ratio = 3.0,
                         cell_size = cs);

  // Mesh generation
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, 
                                      criteria,
                                      CGAL::parameters::no_perturb(),
                                      CGAL::parameters::no_exude());


  // Compute volume
  double mesh_volume = .0;
  std::size_t mesh_cell_count = 0;
  for(C3t3::Cells_in_complex_iterator cit = c3t3.cells_in_complex_begin();
      cit != c3t3.cells_in_complex_end();
      ++cit)
  {

    Tetrahedron_3 t(cit->vertex(0)->point(),
                    cit->vertex(1)->point(),
                    cit->vertex(2)->point(),
                    cit->vertex(3)->point());
    mesh_volume += std::abs(CGAL::to_double(t.volume()));
    mesh_cell_count++;
  }

  std::cout << "Volume of mesh: "       << mesh_volume << std::endl;
  std::cout << "Number of cells:     "  << mesh_cell_count << std::endl;
}

Attachment: bigger.off
Description: Binary data

Attachment: smaller.off
Description: Binary data

Attachment: bigger.png
Description: PNG image

Attachment: smaller.png
Description: PNG image