Issue with CGAL 3D periodic Delaunay triangulation with info

Question

I want to construct a periodic 3D Delaunay triangulation with info (in this case an integer) using CGAL. For 2D this works well if I construct a vector of pairs (point, info) and pass it to the triangulation function. The very analogous code for 3D won't compile however. Below is an example, where points from my own format "particle" are converted to the CGAL point format.

Is this an issue with CGAL or am I missing something?

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/periodic_3_triangulation_3_io.h>

#include <iostream>
#include <fstream>
#include <cassert>
#include <list>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel            K;
typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>            Gt;

typedef CGAL::Periodic_3_triangulation_ds_vertex_base_3<>              VbDS;
typedef CGAL::Triangulation_vertex_base_3<Gt, VbDS>                    Vb;
typedef CGAL::Triangulation_vertex_base_with_info_3<unsigned, Gt, Vb>  VbInfo;

typedef CGAL::Periodic_3_triangulation_ds_cell_base_3<>                CbDS;
typedef CGAL::Triangulation_cell_base_3<Gt, CbDS>                      Cb;

typedef CGAL::Triangulation_data_structure_3<VbInfo, Cb>               Tds;

typedef CGAL::Delaunay_triangulation_3<K, Tds>                         DT3;
typedef CGAL::Periodic_3_Delaunay_triangulation_3<Gt, Tds>             P3DT3;

typedef P3DT3::Point               Point;
typedef P3DT3::Iso_cuboid          Iso_cuboid;

void create_PDT_3D(const std::vector<particle>& grid )
{
     // The cube for the periodic domain
     Iso_cuboid domain(0,0,0,1,1,1); 

     // Convert "particle" format to "point with info" format
     std::vector<Point> points_bare;
     std::vector< std::pair<Point,unsigned> > points_info;

     points_bare.reserve(grid.size());
     points_info.reserve(grid.size());

     // the "info" is just given by the index of the point
     for (unsigned i = 0; i < grid.size(); ++i )
     {
          points_bare.emplace_back( Point(grid[i].x, grid[i].y, grid[i].z) );
          points_info.emplace_back( std::make_pair( Point(grid[i].x, grid[i].y, grid[i].z), i ) );
     }

     // Triangulate
     DT3 Ta    ( points_bare.begin(), points_bare.end() );         // working
     DT3 Tb    ( points_info.begin(), points_info.end() );         // working

     // Triangulate (periodic)
     P3DT3 TTa ( points_bare.begin(), points_bare.end(), domain ); // working
     P3DT3 TTb ( points_info.begin(), points_info.end(), domain ); // NOT working
}

where "particle" is a very simple class

class particle
{
     public:
          double x, y, z;
          particle(double xr=0, double yr=0, double zr=0) : x(xr), y(yr), z(zr) {};
};

Mael Mael · Accepted Answer · 2019-08-28T08:06:32

The 3D periodic triangulations in CGAL do not have the necessary code to insert ranges of points with info. It is not because of a strong obstacle making it impossible to implement, it simply just was not done.

It should however be almost immediate to translate the code from the 2D periodic triangulations to the 3D case: you can look at the function insert_with_info() and how it is used in the class Periodic_2_Delaunay_triangulation_2. The code will be almost identical for 3D.

Let me know how it goes, and if you need help.

Issue with CGAL 3D periodic Delaunay triangulation with info

1 Answers