CGAL users discussion list

[Patrick <>]

Hi,

I've been using CGAL to mesh (approximations to) fattened gyroid surfaces, confined to a box from an implicit function. Attached is a simple example code. Unfortunately, this results in a surface mesh containing a non-manifold edge (with 4 faces attached to it). this problem seems to persist in the tet mesh that was generated (I wrote a little program to recalculate its surface mesh). Am I making some silly mistake here? Is there any way to fix it?

Thanks,

Patrick

#include <iostream>
#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/Implicit_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
//#include "utils.hpp"

using std::abs;
using std::min;
using std::max;
using std::string;

const double pi = 3.14159265358979323846;

template <class type>
inline type
min( type a, type b, type c ) {
  const type m = std::min( a, b );
  return std::min( m, c );
}

template <class type>
inline type
max( type a, type b, type c ) {
  const type m = std::max( a, b );
  return std::max( m, c );
}

template <class type>
inline type
sqr( type a ) {
   return a*a;
}


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef FT (Function)(const Point&);
typedef CGAL::Implicit_mesh_domain_3<Function,K> Mesh_domain;
//#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
//#else
//typedef CGAL::Sequential_tag Concurrency_tag;
//#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;

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

// ---- Start here to define volume ---- 

const FT box_size_x = 4.0; // box goes from x=0 to x=box_size_x!
const FT box_size_y = 4.0; // box goes from y=0 to y=box_size_y!
const FT box_size_z = 4.0; // box goes from z=0 to z=box_size_z!

const FT factor = 1.0;

const FT num_unit_x = 1.0;
const FT num_unit_y = 1.0;
const FT num_unit_z = 1.0;

const FT thickness = 0.8;

FT gyroid (const Point &p) {
  
  const FT x = p.x(), y = p.y(), z = p.z();
  
  const FT scale_x = num_unit_x*2.0*pi / box_size_x;
  const FT scale_y = num_unit_y*2.0*pi / box_size_y;
  const FT scale_z = num_unit_z*2.0*pi / box_size_z;
  
  return abs(cos(scale_x*x) * sin(scale_y*y) + cos(scale_y*y)*sin(scale_z*z) + cos(scale_z*z)*sin(scale_x*x)) - thickness;
}

FT box (const Point &p) {
	const FT x = p.x(), y = p.y(), z = p.z();
	const FT x_center = box_size_x/2.0, y_center = box_size_y/2.0, z_center = box_size_z/2.0;
	return max(2.0/(factor*box_size_x) * abs(x-x_center), 2.0/(factor*box_size_y) * abs(y-y_center), 2.0/(factor*box_size_z) * abs(z-z_center))-1.0;
}

FT impl_function (const Point&p) {
	return max(gyroid(p),box(p));
    //return min(box(p), schwp(p));
	//return min(box(p), ellipse(p));
}

// parameters
const FT fa = 30.0;  // This parameter controls the shape of surface facets. Actually, 
                        // it is a lower bound for the angle (in degree) of 
	                    // surface facets. When boundary surfaces are smooth, the 
                        // termination of the meshing process is guaranteed if the angular bound is at most 30 degrees.
const FT fs = 0.1;    // This parameter controls the size of surface facets. 
                        // Actually, each surface facet has a surface 
                        // Delaunay ball which is a ball circumscribing the surface 
                        // facet and centered on the surface patch. The parameter facet_size 
                        // is either a constant or a spatially variable scalar field, providing an 
                        // upper bound for the radii of surface Delaunay balls.
const FT fd = 0.025; // This parameter controls the approximation error of boundary and 
	                      // subdivision surfaces. Actually, it is either a constant or a spatially
	                      // variable scalar field. It provides an upper bound for the distance 
	                      // between the circumcenter of a surface facet and the center of a 
	                      // surface Delaunay ball of this facet.
// facet_topology = 1.0   // This parameters controls the set of topological constraints which have to 
	                      // be verified by each surface facet. By default, each vertex of a surface facet has 
	                      // to be located on a surface patch, on a curve segment, or on a corner. It can also 
	                      // be set to check whether the three vertices of a surface facet belongs to the same 
	                      // surface patch. This has to be done cautiously, as such a criteria needs that 
	                      // each surface patches intersection is an input 1-dimensional feature.
const FT crer = 1.8; // This parameter controls the shape of mesh cells (but can't filter slivers, as 
	                            // we discussed earlier). Actually, it is an upper bound for the ratio between the 
	                            // circumradius of a mesh tetrahedron and its shortest edge. There is a theoretical 
	                            // bound for this parameter: the Delaunay refinement process is guaranteed to 
	                            // terminate for values of cell_radius_edge_ratio bigger than 2.
const FT cs = 0.1;             // This parameter controls the size of mesh tetrahedra. It is either a scalar or a 
                                // spatially variable scalar field. It provides an upper bound on the circumradii of
                                // the mesh tetrahedra.

int main()
{
  // Domain (Warning: Sphere_3 constructor uses squared radius !)
  const Point center(box_size_x/2.0, box_size_y/2.0, box_size_z/2.0);
  const FT radsq = sqr(factor*box_size_x/2.0) + sqr(factor*box_size_y/2.0) + sqr(factor*box_size_z/2.0)+0.1;
  //std::cerr << impl_function(center);
  Mesh_domain domain(impl_function,
                     K::Sphere_3(center, radsq));
  // Mesh criteria
  Mesh_criteria criteria(facet_angle=fa, facet_size=fs, facet_distance=fd,
                         cell_radius_edge_ratio=crer, cell_size=cs);
  
  // Mesh generation
  //C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, odt(), perturb(), exude());
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
  // Output
  std::ofstream medit_file("out.mesh");
  c3t3.output_to_medit(medit_file);
  medit_file.close();
  std::ofstream off_file("out.off");
  c3t3.output_boundary_to_off(off_file);
  off_file.close();
  return 0;
}