Add remesh_surface utility

include/deal.II/cgal/triangulation.h

@@ -27,6 +27,7 @@
 
 #ifdef DEAL_II_WITH_CGAL
 
+
 #  include <boost/hana.hpp>
 
 #  include <CGAL/Complex_2_in_triangulation_3.h>
@@ -181,6 +182,84 @@ namespace CGALWrappers
   cgal_surface_mesh_to_dealii_triangulation(const CGAL_MeshType &cgal_mesh,
                                             Triangulation<2, 3> &triangulation);
 
+  /**
+   * Generate a Triangulation from the zero level set of an implicit function,
+   * using the CGAL library.
+   *
+   * This function is only implemented for `dim` equal to two or three, and
+   * requires that deal.II is configured using `DEAL_II_WITH_CGAL`. When `dim`
+   * is equal to three, the @p implicit_function is supposed to be negative in
+   * the interior of the domain, positive outside, and to be entirely enclosed
+   * in a ball of radius @p outer_ball_radius centered at the point
+   * @p interior_point. The triangulation that is generated covers the volume
+   * bounded by the zero level set of the implicit function  where the
+   * @p implicit_function is negative.
+   *
+   * When `dim` is equal to two, the generated surface triangulation is the zero
+   * level set of the @p implicit_function, oriented such that the surface
+   * triangulation has normals pointing towards the region where
+   * @p implicit_function is positive.
+   *
+   * The struct @p data can be used to pass additional
+   * arguments to the CGAL::Mesh_criteria_3 class (see
+   * https://doc.cgal.org/latest/Mesh_3/index.html for more information.)
+   *
+   * An example usage of this function is given by
+   *
+   * @code
+   * Triangulation<dim, 3>  tria;
+   * FunctionParser<3> my_function("(1-sqrt(x^2+y^2))^2+z^2-.25");
+   * GridGenerator::implicit_function( tria, my_function,
+   *      Point<3>(.5, 0, 0), 1.0, cell_size = 0.2);
+   * @endcode
+   *
+   * The above snippet of code generates the following grid for `dim` equal to
+   * two and three respectively
+   *
+   * @image html grid_generator_implicit_function_2d.png
+   *
+   * @image html grid_generator_implicit_function_3d.png
+   *
+   * @ingroup simplex
+   *
+   * @param[out] tria The output triangulation
+   * @param[in] implicit_function The implicit function
+   * @param[in] data Additional parameters to pass to the CGAL::make_mesh_3
+   * function and to the CGAL::make_surface_mesh functions
+   * @param[in] interior_point A point in the interior of the domain, for which
+   * @p implicit_function is negative
+   * @param[in] outer_ball_radius The radius of the ball that will contain the
+   * generated Triangulation object
+   */
+  template <int dim>
+  void
+  implicit_function(Triangulation<dim, 3> &                  tria,
+                    const Function<3> &                      implicit_function,
+                    const CGALWrappers::AdditionalData<dim> &data =
+                      CGALWrappers::AdditionalData<dim>{},
+                    const Point<3> &interior_point    = Point<3>(),
+                    const double &  outer_ball_radius = 1.0);
+
+  /**
+   * Create a deal.II Triangulation<3> out of a deal.II Triangulation<2,3>
+   * by filling it with tetrahedra.
+   *
+   * The last optional argument @p data can be used to pass additional
+   * arguments to the CGAL::Mesh_criteria_3 class (see
+   * https://doc.cgal.org/latest/Mesh_3/index.html for more information).
+   *
+   *
+   * @param [in] surface_tria The input deal.II Triangulation<2,3>.
+   * @param [out] vol_tria The output deal.II Triangulation<3>.
+   * @param[in] data Additional parameters to pass to the CGAL::make_mesh_3
+   * function and to the CGAL::make_surface_mesh functions.
+   */
+  void
+  surface_mesh_to_volumetric_mesh(const Triangulation<2, 3> &surface_tria,
+                                  Triangulation<3> &         vol_tria,
+                                  const CGALWrappers::AdditionalData<3> &data =
+                                    CGALWrappers::AdditionalData<3>{});
+
 
 
 #  ifndef DOXYGEN

include/deal.II/cgal/utilities.h

@@ -39,16 +39,21 @@
 #  include <CGAL/Mesh_triangulation_3.h>
 #  include <CGAL/Polygon_mesh_processing/corefinement.h>
 #  include <CGAL/Polygon_mesh_processing/measure.h>
+#  include <CGAL/Polygon_mesh_processing/remesh.h>
 #  include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 #  include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
 #  include <CGAL/Simple_cartesian.h>
 #  include <CGAL/Surface_mesh.h>
 #  include <CGAL/Triangulation_3.h>
+#  include <CGAL/boost/graph/copy_face_graph.h>
+#  include <CGAL/boost/graph/selection.h>
 #  include <CGAL/convex_hull_3.h>
 #  include <CGAL/make_mesh_3.h>
 #  include <CGAL/make_surface_mesh.h>
 #  include <deal.II/cgal/surface_mesh.h>
+//#  include <CGAL/tetrahedral_remeshing.h> REQUIRES CGAL_VERSION>=5.1.5
 
+#  include <fstream>
 #  include <type_traits>
 
 
@@ -255,6 +260,39 @@ namespace CGALWrappers
     const unsigned int                                                   degree,
     const Mapping<dim0, spacedim> &mapping0,
     const Mapping<dim1, spacedim> &mapping1);
+
+  /**
+   * Remesh a CGAL::Surface_mesh.
+   *
+   * If the domain has 1-dimensional exposed features, the criteria includes a
+   * sizing field to guide the sampling of 1-dimensional features with
+   * protecting balls centers.
+   * - CGAL::parameters::edge_size.
+   *
+   * This utility should be exploited to improve the quality of a mesh coming
+   * from boolean operations. As an example, consider two CGAL::Surface_mesh
+   * objects describing two hyper_spheres that intersect non-trivially. After
+   * computing their corefinement and union with compute_boolean_operation(),
+   * the resulting mesh is the following:
+   *
+   * @image html boolean_union_hyper_spheres.png
+   *
+   * Clearly, elements (triangles) like the ones arising along the intersection
+   * of the two spheres should be avoided as they're quite degenerate and would
+   * decrease the accuracy of the results. A call to the present function will
+   * try to remesh the surface, according to the given named parameters, to
+   * improve its quality. In the present case, the result is the following:
+   *
+   * @image html boolean_union_hyper_spheres_remeshed.png
+   *
+   * @param surface_mesh The input CGAL::Surface_mesh.
+   * @param [in] data AdditionalData object to pass to the CGAL::make_mesh_3 function. See the documentation
+   * of that struct for a description of those parameters.
+   */
+  template <typename CGALPointType>
+  void
+  remesh_surface(CGAL::Surface_mesh<CGALPointType> &surface_mesh,
+                 const AdditionalData<3> &          data = AdditionalData<3>{});
 } // namespace CGALWrappers
 
 #  ifndef DOXYGEN
@@ -528,6 +566,80 @@ namespace CGALWrappers
     tr.insert(dummy.points().begin(), dummy.points().end());
     return compute_quadrature(tr, degree);
   }
+
+
+
+  template <typename CGALPointType>
+  void
+  remesh_surface(CGAL::Surface_mesh<CGALPointType> &cgal_mesh,
+                 const AdditionalData<3> &          data)
+  {
+    using K           = CGAL::Exact_predicates_inexact_constructions_kernel;
+    using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
+    // Polyhedron type
+    using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
+    // Triangulation
+    using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
+    using C3t3 =
+      CGAL::Mesh_complex_3_in_triangulation_3<Tr,
+                                              Mesh_domain::Corner_index,
+                                              Mesh_domain::Curve_index>;
+    // Criteria
+    using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
+
+    Polyhedron poly;
+    CGAL::copy_face_graph(cgal_mesh, poly);
+    // Create a vector with only one element: the pointer to the
+    // polyhedron.
+    std::vector<Polyhedron *> poly_ptrs_vector(1, &poly);
+    // Create a polyhedral domain, with only one polyhedron,
+    // and no "bounding polyhedron", so the volumetric part of the
+    // domain will be empty.
+    Mesh_domain domain(poly_ptrs_vector.begin(), poly_ptrs_vector.end());
+    // Get sharp features
+    domain.detect_features(); // includes detection of borders
+
+    // Mesh criteria
+    if (data.edge_size > 0)
+      {
+        Mesh_criteria criteria(CGAL::parameters::edge_size   = data.edge_size,
+                               CGAL::parameters::facet_size  = data.facet_size,
+                               CGAL::parameters::facet_angle = data.facet_angle,
+                               CGAL::parameters::facet_distance =
+                                 data.facet_distance,
+                               CGAL::parameters::cell_radius_edge_ratio =
+                                 data.cell_radius_edge_ratio,
+                               CGAL::parameters::cell_size = data.cell_size);
+        // Mesh generation
+        C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
+                                            criteria,
+                                            CGAL::parameters::no_perturb(),
+                                            CGAL::parameters::no_exude());
+        cgal_mesh.clear();
+        CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
+      }
+    else if (data.edge_size == CGAL::parameters::is_default_parameter(true))
+      {
+        Mesh_criteria criteria(CGAL::parameters::facet_size  = data.facet_size,
+                               CGAL::parameters::facet_angle = data.facet_angle,
+                               CGAL::parameters::facet_distance =
+                                 data.facet_distance,
+                               CGAL::parameters::cell_radius_edge_ratio =
+                                 data.cell_radius_edge_ratio,
+                               CGAL::parameters::cell_size = data.cell_size);
+        // Mesh generation
+        C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
+                                            criteria,
+                                            CGAL::parameters::no_perturb(),
+                                            CGAL::parameters::no_exude());
+        cgal_mesh.clear();
+        CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
+      }
+    else
+      {
+        Assert(false, ExcInternalError());
+      }
+  }
 } // namespace CGALWrappers
 #  endif
 

include/deal.II/grid/grid_generator.h

@@ -1820,11 +1820,7 @@ namespace GridGenerator
    * @param [in] surface_tria The input deal.II Triangulation<2,3>.
    * @param [out] vol_tria The output deal.II Triangulation<3>.
    * @param[in] data Additional parameters to pass to the CGAL::make_mesh_3
-<<<<<<< HEAD
    * function.
-=======
-   * function and to the CGAL::make_surface_mesh functions.
->>>>>>> 2337fe411f (Surface deal.II tria to volumetric deal.II tria)
    */
   void
   surface_mesh_to_volumetric_mesh(const Triangulation<2, 3> &surface_tria,

tests/cgal/cgal_remesh_surface.cc

@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Remesh the union of two deal.II hyper_spheres in order to avoid bad
+// triangles.
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/surface_mesh.h>
+#include <deal.II/cgal/triangulation.h>
+#include <string.h>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+using K         = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_3<K>;
+using CGALMesh  = CGAL::Surface_mesh<CGALPoint>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl;
+
+  Triangulation<spacedim> tria0, tria1;
+  Triangulation<2, 3>     tria_out;
+  GridOut                 go;
+  CGALMesh                surface_mesh0, surface_mesh1, out_mesh;
+
+  GridGenerator::hyper_ball(tria0, {0., 0., 0.}, 0.4);
+  GridGenerator::hyper_ball(tria1, {0.3, 0.3, 0.3}, 0.4);
+  tria0.refine_global(3);
+  tria1.refine_global(3);
+
+  // Move to CGAL surfaces
+  dealii_tria_to_cgal_surface_mesh(tria0, surface_mesh0);
+  dealii_tria_to_cgal_surface_mesh(tria1, surface_mesh1);
+
+  // close the surfaces
+  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh0);
+  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh1);
+
+  CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0);
+  CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1);
+
+  compute_boolean_operation(surface_mesh0,
+                            surface_mesh1,
+                            BooleanOperation::compute_union,
+                            out_mesh);
+  AdditionalData<3> data;
+  data.edge_size      = .02;
+  data.facet_angle    = 25;
+  data.facet_size     = 0.05;
+  data.facet_distance = 0.05;
+
+  remesh_surface<CGALPoint>(out_mesh, data);
+
+  //  Now back to deal.II
+  cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out);
+  std::ofstream out_name_spheres("remeshed_union_spheres.vtk");
+  go.write_vtk(tria_out, out_name_spheres);
+  deallog << "OK" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+  test<3, 3>();
+}

tests/cgal/cgal_remesh_surface.output

@@ -0,0 +1,3 @@
+
+DEAL::dim= 3,	 spacedim= 3
+DEAL::OK

tests/cgal/cgal_surface_mesh_05.cc

@@ -70,10 +70,10 @@ test()
                             out_mesh);
   // Now back to deal.II
   cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out);
-  std::ofstream out_name_spheres("boolean_triangulation_spheres.vtk");
+  std::ofstream out_name_spheres("boolean_intersection_hyper_spheres.vtk");
   go.write_vtk(tria_out, out_name_spheres);
   deallog << "OK" << std::endl;
-  remove("tria_out.vtk");
+  remove("boolean_intersection_hyper_spheres.vtk");
 
   // Clear everything
   tria0.clear();
@@ -106,10 +106,10 @@ test()
                             out_mesh);
   // Now back to deal.II
   cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out);
-  std::ofstream out_name_cubes("boolean_triangulation_cubes.vtk");
+  std::ofstream out_name_cubes("boolean_intersection_cubes.vtk");
   go.write_vtk(tria_out, out_name_cubes);
-  remove("tria_out.vtk");
   deallog << "OK" << std::endl;
+  remove("boolean_intersection_cubes.vtk");
 }
 
 int