Check edge cases of CGAL::Polygon_mesh_processing::corefine_and_compute_union

/*
  https://github.com/openscad/openscad/pull/3641

  g++ \
    -stdlib=libc++ -std=c++1y \
    -I../CGAL-5.2/include \
    -lgmp -lmpfr \
    test_pmp.cc \
    -o test_pmp && ./test_pmp
*/

#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/draw_polyhedron.h>
#include <fstream>
// typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
typedef CGAL::Epeck Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::Halfedge_handle Halfedge_handle;

// Adapted from https://doc.cgal.org/latest/Polyhedron/Polyhedron_2polyhedron_prog_cube_8cpp-example.html
// (added ugly dx, dy, dz)
template <class Poly>
typename Poly::Halfedge_handle make_cube_3( Poly& P, double sx, double sy, double sz, double dx, double dy, double dz) {
    // appends a cube of size [0,1]^3 to the polyhedron P.
    CGAL_precondition( P.is_valid());
    typedef typename Poly::Point_3         Point;
    typedef typename Poly::Halfedge_handle Halfedge_handle;
    Halfedge_handle h = P.make_tetrahedron( Point( sx + dx, 0 + dy, 0 + dz),
                                            Point( 0 + dx, 0 + dy, sz + dz),
                                            Point( 0 + dx, 0 + dy, 0 + dz),
                                            Point( 0 + dx, sy + dy, 0 + dz));
    Halfedge_handle g = h->next()->opposite()->next();             // Fig. (a)
    P.split_edge( h->next());
    P.split_edge( g->next());
    P.split_edge( g);                                              // Fig. (b)
    h->next()->vertex()->point()     = Point( sx + dx, 0 + dy, sz + dz);
    g->next()->vertex()->point()     = Point( 0 + dx, sy + dy, sz + dz);
    g->opposite()->vertex()->point() = Point( sx + dx, sy + dy, 0 + dz);            // Fig. (c)
    Halfedge_handle f = P.split_facet( g->next(),
                                       g->next()->next()->next()); // Fig. (d)
    Halfedge_handle e = P.split_edge( f);
    e->vertex()->point() = Point( sx + dx, sy + dy, sz + dz);                        // Fig. (e)
    P.split_facet( e, f->next()->next());                          // Fig. (f)
    CGAL_postcondition( P.is_valid());
    return h;
}


int main(int argc, char* argv[])
{
  Polyhedron P1; make_cube_3(P1, 1, 1, 1, 0, 0, 0); // cube(1);

  Polyhedron P2; make_cube_3(P2, 1, 1, 1, 1, 1, 0); // translate([1, 1, 0]) cube(1); // BREAKS
  // Polyhedron P2; make_cube_3(P2, 0.5, 0.5, 0.5, 1, 1, 1); // translate([1, 1, 1]) cube(0.5); // OK
  // Polyhedron P2; make_cube_3(P2, 0.5, 0.5, 0.5, 1, 1, 0.25); // translate([1, 1, 0.25]) cube(0.5); // BREAKS
  // Polyhedron P2; make_cube_3(P2, 0.5, 0.5, 0.5, 1, 0, 0.25); // translate([1, 0, 0.25]) cube(0.5); // BREAKS

  CGAL::Polygon_mesh_processing::triangulate_faces(P1);
  CGAL::Polygon_mesh_processing::triangulate_faces(P2);

  CGAL::Polygon_mesh_processing::corefine_and_compute_union(P1, P2, P1);

  std::ofstream out("out.off");
  out << P1;
  return EXIT_SUCCESS;
}