inter_tets.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Surface_mesh.h>
#include <math.h>
#include <fstream>
#include <iostream>

typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef K::Point_3 Point_3;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;

namespace PMP = CGAL::Polygon_mesh_processing;

int main()
{
  double phi = (1.0 + sqrt(5.0)) / 2.0;
  double a = 1.0 / sqrt(3.0);
  double b = a / phi;
  double c = a * phi;

  Mesh mesh1, mesh2, mesh3, mesh4, mesh5;

  
  std::vector<Point_3> g_pts =
  {
        Point_3(a, a, a),
        Point_3(a, a, -a),
        Point_3(a, -a, a),
        Point_3(-a, -a, a),
        Point_3(-a, a, -a),
        Point_3(-a, a, a),
        Point_3(0, b, -c),
        Point_3(0, -b, -c),
        Point_3(0, -b, c),
        Point_3(c, 0, -b),
        Point_3(-c, 0, -b),
        Point_3(-c, 0, b),
        Point_3(b, c, 0),
        Point_3(b, -c, 0),
        Point_3(-b, -c, 0),
        Point_3(-b, c, 0),
        Point_3(0, b, c),
        Point_3(a, -a, -a),
        Point_3(c, 0, b),
        Point_3(-a, -a, -a),
  };
  
  std::vector< std::array<int,3> >poly=
  {
    {0, 1, 2}, 
    {2, 1, 3}, 
    {3, 1, 0}, 
    {0, 2, 3},
  };
  
  std::vector< std::array<int,4> > pts_ids = {
    {16, 13, 1, 10},
    {17, 0, 3, 4},
    {18, 5, 14, 6},
    {2, 12, 11, 7},
    {19, 15, 9, 8}
  };

  std::vector<Mesh> meshes(5);
  for (int i=0; i<5; ++i)
  {
    std::vector<Point_3> pts;
    for (int k : pts_ids[i])
      pts.push_back(g_pts[k]);
    CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(pts, poly,meshes[i]);
    std::ofstream("mesh_"+std::to_string(i)+".off") << std::setprecision(17) << meshes[i];
  }

  Mesh res=meshes[0];
  for (int i=1;i<5;++i)
    PMP::corefine_and_compute_intersection(meshes[i], res, res);


  std::ofstream("res.off") << std::setprecision(17) << res;
}