GridGenerator: put CGAL implementation into own compilation unit

diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt
index 88a236d06eae6efd8788b5569d07b4310b49b78f..0d13509ef3cf3932b36263eb98b3886d02475749 100644 (file)
--- a/source/grid/CMakeLists.txt
+++ b/source/grid/CMakeLists.txt
@@ -21,12 +21,12 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 # 100k characters produced by one of the CGAL headers. Unfortunately,
 # guarding the include by DEAL_II_DISABLE_EXTRA_DIAGNOSTICS is not enough
 # due to a longstanding bug in gcc. Thus, simply wet
-# -Wno-misleading-indentation on the command line for the grid_generator.cc
-# compilation unit.
+# -Wno-misleading-indentation on the command line for the
+# grid_generator_cgal.cc compilation unit.
 #
 IF(DEAL_II_WITH_CGAL)
   ENABLE_IF_SUPPORTED(_flag "-Wno-misleading-indentation")
-  SET_PROPERTY(SOURCE "grid_generator.cc"
+  SET_PROPERTY(SOURCE "grid_generator_cgal.cc"
     APPEND_STRING PROPERTY COMPILE_FLAGS " ${_flag}"
     )
 ENDIF()
@@ -49,15 +49,16 @@ SET(_unity_include_src
   )
 
 SET(_separate_src
+  grid_generator.cc
+  grid_generator_cgal.cc
   grid_generator_from_name.cc
+  grid_generator_pipe_junction.cc
   grid_in.cc
   grid_out.cc
-  grid_generator.cc
-  grid_generator_pipe_junction.cc
-  grid_tools.cc
   grid_tools_cache.cc
-  grid_tools_nontemplates.cc
+  grid_tools.cc
   grid_tools_dof_handlers.cc
+  grid_tools_nontemplates.cc
   tria.cc
   )
 
@@ -72,8 +73,9 @@ SETUP_SOURCE_LIST("${_unity_include_src}"
 
 SET(_inst
   cell_id.inst.in
-  grid_generator.inst.in
+  grid_generator_cgal.inst.in
   grid_generator_from_name.inst.in
+  grid_generator.inst.in
   grid_generator_pipe_junction.inst.in
   grid_in.inst.in
   grid_out.inst.in

diff --git a/source/grid/grid_generator_cgal.cc b/source/grid/grid_generator_cgal.cc
new file mode 100644 (file)
index 0000000..ed9e748
--- /dev/null
+++ b/source/grid/grid_generator_cgal.cc
@@ -0,0 +1,155 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1999 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/ndarray.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#ifdef DEAL_II_WITH_CGAL
+// Functions needed by the CGAL mesh generation utilities are inside
+#  include <deal.II/cgal/triangulation.h>
+#endif
+
+
+DEAL_II_NAMESPACE_OPEN
+
+// work around the problem that doxygen for some reason lists all template
+// specializations in this file
+#ifndef DOXYGEN
+
+namespace GridGenerator
+{
+  template <int dim>
+  void
+  implicit_function(Triangulation<dim, 3> &tria,
+                    const Function<3> &    dealii_implicit_function,
+                    const CGALWrappers::AdditionalData<dim> &data,
+                    const Point<3> &                         interior_point,
+                    const double &                           outer_ball_radius)
+  {
+#  ifdef DEAL_II_WITH_CGAL
+    Assert(dealii_implicit_function.n_components == 1,
+           ExcMessage(
+             "The implicit function must have exactly one component."));
+    Assert(dealii_implicit_function.value(interior_point) < 0,
+           ExcMessage(
+             "The implicit function must be negative at the interior point."));
+    Assert(outer_ball_radius > 0,
+           ExcMessage("The outer ball radius must be positive."));
+    Assert(tria.n_active_cells() == 0,
+           ExcMessage("The triangulation must be empty."));
+
+    if constexpr (dim == 3)
+      {
+        using K          = CGAL::Exact_predicates_inexact_constructions_kernel;
+        using NumberType = K::FT;
+        using Point_3    = K::Point_3;
+        using Sphere_3   = K::Sphere_3;
+
+        using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
+        using Tr =
+          CGAL::Mesh_triangulation_3<Mesh_domain,
+                                     CGAL::Default,
+                                     CGALWrappers::ConcurrencyTag>::type;
+        using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
+        using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
+
+
+        auto cgal_implicit_function = [&](const Point_3 &p) {
+          return NumberType(
+            dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
+        };
+
+        Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
+          cgal_implicit_function,
+          K::Sphere_3(
+            Point_3(interior_point[0], interior_point[1], interior_point[2]),
+            outer_ball_radius * outer_ball_radius));
+
+        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);
+
+        auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
+        CGALWrappers::cgal_triangulation_to_dealii_triangulation(
+          cgal_triangulation, tria);
+      }
+    else if constexpr (dim == 2)
+      {
+        // default triangulation for Surface_mesher
+        using Tr       = CGAL::Surface_mesh_default_triangulation_3;
+        using C2t3     = CGAL::Complex_2_in_triangulation_3<Tr>;
+        using GT       = Tr::Geom_traits;
+        using Sphere_3 = GT::Sphere_3;
+        using Point_3  = GT::Point_3;
+        using FT       = GT::FT;
+        typedef FT (*Function)(Point_3);
+        using Surface_3    = CGAL::Implicit_surface_3<GT, Function>;
+        using Surface_mesh = CGAL::Surface_mesh<Point_3>;
+
+
+        auto cgal_implicit_function = [&](const Point_3 &p) {
+          return FT(
+            dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
+        };
+
+        Surface_3 surface(cgal_implicit_function,
+                          Sphere_3(Point_3(interior_point[0],
+                                           interior_point[1],
+                                           interior_point[2]),
+                                   outer_ball_radius * outer_ball_radius));
+
+        Tr           tr;
+        C2t3         c2t3(tr);
+        Surface_mesh mesh;
+
+        CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
+                                                           data.radius_bound,
+                                                           data.distance_bound);
+        CGAL::make_surface_mesh(c2t3,
+                                surface,
+                                criteria,
+                                CGAL::Non_manifold_tag());
+        CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
+        CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
+      }
+    else
+      {
+        Assert(false, ExcImpossibleInDim(dim));
+      }
+
+#  else
+
+    (void)tria;
+    (void)dealii_implicit_function;
+    (void)data;
+    (void)interior_point;
+    (void)outer_ball_radius;
+    AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));
+
+#  endif
+  }
+} // namespace GridGenerator
+
+// explicit instantiations
+#  include "grid_generator_cgal.inst"
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/source/grid/grid_generator_cgal.inst.in b/source/grid/grid_generator_cgal.inst.in
new file mode 100644 (file)
index 0000000..0db59c8
--- /dev/null
+++ b/source/grid/grid_generator_cgal.inst.in
@@ -0,0 +1,30 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1999 - 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.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    namespace GridGenerator
+    \{
+#if deal_II_dimension >= 2
+      template void
+      implicit_function<deal_II_dimension>(
+        Triangulation<deal_II_dimension, 3> &,
+        const Function<3> &,
+        const CGALWrappers::AdditionalData<deal_II_dimension> &,
+        const Point<3> &,
+        const double &);
+#endif
+    \}
+  }