]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added Implementation for Triangulation<2,3> 13743/head
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 17 May 2022 16:29:09 +0000 (19:29 +0300)
committerMarco Feder <marco.feder@sissa.it>
Sat, 4 Jun 2022 14:32:10 +0000 (16:32 +0200)
doc/doxygen/images/grid_generator_implicit_function_2d.png [new file with mode: 0644]
include/deal.II/cgal/utilities.h
include/deal.II/grid/grid_generator.h
tests/cgal/cgal_surface_triangulation_03.cc [new file with mode: 0644]
tests/cgal/cgal_surface_triangulation_03.output [new file with mode: 0644]
tests/cgal/implicit_triangulation_01.cc
tests/cgal/implicit_triangulation_02.cc [new file with mode: 0644]
tests/cgal/implicit_triangulation_02.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/grid_generator_implicit_function_2d.png b/doc/doxygen/images/grid_generator_implicit_function_2d.png
new file mode 100644 (file)
index 0000000..a61a23e
Binary files /dev/null and b/doc/doxygen/images/grid_generator_implicit_function_2d.png differ
index 5758429575a3189980648cc8b019a5ee90feef7d..ae8b7178f9797d8798e22f21fd36486a069b02b1 100644 (file)
@@ -113,6 +113,102 @@ namespace CGALWrappers
     compute_union = 1 << 3,
   };
 
+
+
+  /**
+   * Struct that must 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.)
+   *
+   * The arguments allow for fine control on the size, quality, and distribution
+   * of the cells of the final triangulation. CGAL uses Boost named parameters
+   * for these arguments in dimension three, i.e., they must be specified with
+   * the syntax `CGAL::parameters::parameter_name=parameter_value`, irrespective
+   * of their order. Accepted parameters are:
+   *
+   * - `CGAL::parameters::edge_size`: a constant
+   * providing a uniform upper bound for the lengths of
+   * curve edges. This parameter has to be set to a positive value when
+   * 1-dimensional features protection is used.
+   * - `CGAL::parameters::facet_angle`: a lower bound for the angles (in
+   * degrees) of the surface mesh facets.
+   * - `CGAL::parameters::facet_size`: a constant
+   * describing a uniform upper bound for the radii
+   * of the surface Delaunay balls.
+   * - `CGAL::parameters::facet_distance`: a constant
+   * describing a uniform upper bound for the distance
+   * between the facet circumcenter and the center of its surface Delaunay ball.
+   * - `CGAL::parameters::facet_topology`: the set of topological constraints
+   * which have to be verified by each surface facet. The default value is
+   * `CGAL::FACET_VERTICES_ON_SURFACE`. See CGAL::Mesh_facet_topology manual
+   * page to get all possible values.
+   * - `CGAL::parameters::cell_radius_edge_ratio`: an upper bound for the
+   * radius-edge ratio of the mesh tetrahedra.
+   * - `CGAL::parameters::cell_size`: a constant
+   * describing a uniform upper bound for the
+   * circumradii of the mesh tetrahedra.
+   *
+   * @note This struct must be instantiated with `dim=3`.
+   */
+  template <int dim>
+  struct AdditionalData
+  {
+    double facet_angle;
+    double facet_size;
+    double facet_distance;
+    double cell_radius_edge_ratio;
+    double cell_size;
+
+    AdditionalData(
+      double facet_a            = CGAL::parameters::is_default_parameter(true),
+      double facet_s            = CGAL::parameters::is_default_parameter(true),
+      double facet_d            = CGAL::parameters::is_default_parameter(true),
+      double cell_radius_edge_r = CGAL::parameters::is_default_parameter(true),
+      double cell_s             = CGAL::parameters::is_default_parameter(true))
+    {
+      AssertThrow(
+        dim == 3,
+        ExcMessage(
+          "These struct can be instantiated with 3D Triangulations only."));
+      facet_angle            = facet_a;
+      facet_size             = facet_s;
+      facet_distance         = facet_d;
+      cell_radius_edge_ratio = cell_radius_edge_r;
+      cell_size              = cell_s;
+    }
+  };
+
+
+
+  /**
+   * Specialization of the above struct when the object to be constructed is a
+   * 2D triangulation embedded in the 3D space, i.e. a Triangulation<2,3>.
+   * Only three parameters are accepted:
+   * - `angular_bound` is a lower bound in degrees for the angles of mesh
+   * facets.
+   * - `radius_bound` is an upper bound on the radii of surface Delaunay balls.
+   * A surface Delaunay ball is a ball circumscribing a mesh facet and centered
+   * on the surface.
+   * - `distance_bound` is an upper bound for the distance between the
+   * circumcenter of a mesh facet and the center of a surface Delaunay ball of
+   * this facet.
+   */
+  template <>
+  struct AdditionalData<2>
+  {
+    double angular_bound, radius_bound, distance_bound;
+    AdditionalData(double angular_b  = 0.,
+                   double radius_b   = 0.,
+                   double distance_b = 0.)
+    {
+      angular_bound  = angular_b;
+      radius_bound   = radius_b;
+      distance_bound = distance_b;
+    }
+  };
+
+
+
   /**
    * Convert from a deal.II Point to any compatible CGAL point.
    *
index 258e97b2860061b261ed76bca86c3301f70a3aca..1035e1f0921e034f78b3c5a3136fda0e6054a39e 100644 (file)
 
 #ifdef DEAL_II_WITH_CGAL
 // All functions needed by the CGAL mesh generation utilities
+#  include <CGAL/Complex_2_in_triangulation_3.h>
+#  include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
+#  include <CGAL/Implicit_surface_3.h>
 #  include <CGAL/Labeled_mesh_domain_3.h>
 #  include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 #  include <CGAL/Mesh_criteria_3.h>
 #  include <CGAL/Mesh_triangulation_3.h>
+#  include <CGAL/Surface_mesh.h>
+#  include <CGAL/Surface_mesh_default_triangulation_3.h>
 #  include <CGAL/make_mesh_3.h>
+#  include <CGAL/make_surface_mesh.h>
 #  include <deal.II/cgal/triangulation.h>
 #  include <deal.II/cgal/utilities.h>
 #endif
@@ -1765,70 +1771,58 @@ namespace GridGenerator
    * 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, and positive outside, and the
-   * triangulation that is generated covers the volume bounded by the zero level
-   * set of the implicit function and the boundary of a ball of radius,
-   * @p outer_ball_radius where the @p implicit_function is negative.
+   * 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, with or without boundary, depending
-   * on whether the zero level set intersects the boundary of the outer ball.
+   * 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 orientation is such that the surface triangulation has normals pointing
-   * towards the region where @p implicit_function is positive.
-   *
-   * The optional variable arguments @p cgal_args can be used to pass additional
+   * 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.
-   *
-   * The arguments allow fine control on the size, quality, and distribution of
-   * the cells of the final triangulation. CGAL uses named parameters for these
-   * arguments, i.e., they must be specified with the syntax
-   * `CGAL::parameters::parameter_name=parameter_value`. Accepted parameters
-   * are:
-   *
-   * - CGAL::parameters::facet_angl,
-   * - CGAL::parameters::facet_size,
-   * - CGAL::parameters::facet_distance,
-   * - CGAL::parameters::cell_radius_edge_ratio,
-   * - CGAL::parameters::cell_size;
+   * https://doc.cgal.org/latest/Mesh_3/index.html for more information.)
    *
    * An example usage of this function is given by
    *
    * @code
-   * Triangulation<3>  tria;
-   * FunctionParser<3> implicit_function("(1-sqrt(x^2+y^2))^2+z^2-.25");
-   * GridGenerator::implicit_function(
-   *     tria,
-   *     implicit_function,
-   *     Point<3>(1, 0, 0),
-   *     10.0,
-   *     CGAL::parameters::cell_size = 0.2);
+   * 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 image:
+   * 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
-   * @param[in] cgal_args Additional parameters to pass to the CGAL::make_mesh_3
-   * function and to the CGAL::make_surface_mesh functions
    */
-  template <int dim, typename... Args>
+  template <int dim>
   void
-  implicit_function(Triangulation<dim, 3> &tria,
-                    const Function<3> &    implicit_function,
-                    const Point<3> &       interior_point    = Point<3>(),
-                    const double &         outer_ball_radius = 1.0,
-                    Args... cgal_args);
+  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);
 #endif
   ///@}
 
@@ -2740,18 +2734,18 @@ namespace GridGenerator
                         const bool);
 
 #  ifdef DEAL_II_WITH_CGAL
-  template <int dim, typename... Args>
+  template <int dim>
   void
   implicit_function(Triangulation<dim, 3> &tria,
-                    const Function<3> &    implicit_function,
-                    const Point<3> &       interior_point,
-                    const double &         outer_ball_radius,
-                    Args... cgal_args)
+                    const Function<3> &    dealii_implicit_function,
+                    const CGALWrappers::AdditionalData<dim> &data,
+                    const Point<3> &                         interior_point,
+                    const double &                           outer_ball_radius)
   {
-    Assert(implicit_function.n_components == 1,
+    Assert(dealii_implicit_function.n_components == 1,
            ExcMessage(
              "The implicit function must have exactly one component."));
-    Assert(implicit_function.value(interior_point) < 0,
+    Assert(dealii_implicit_function.value(interior_point) < 0,
            ExcMessage(
              "The implicit function must be negative at the interior point."));
     Assert(outer_ball_radius > 0,
@@ -2759,14 +2753,14 @@ namespace GridGenerator
     Assert(tria.n_active_cells() == 0,
            ExcMessage("The triangulation must be empty."));
 
-    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>;
-
     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,
@@ -2777,7 +2771,7 @@ namespace GridGenerator
 
         auto cgal_implicit_function = [&](const Point_3 &p) {
           return NumberType(
-            implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
+            dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
         };
 
         Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
@@ -2786,13 +2780,57 @@ namespace GridGenerator
             Point_3(interior_point[0], interior_point[1], interior_point[2]),
             outer_ball_radius * outer_ball_radius));
 
-        Mesh_criteria criteria(cgal_args...);
+        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));
diff --git a/tests/cgal/cgal_surface_triangulation_03.cc b/tests/cgal/cgal_surface_triangulation_03.cc
new file mode 100644 (file)
index 0000000..5edd6e4
--- /dev/null
@@ -0,0 +1,63 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Convert a closed CGAL::Surface_mesh or a closed CGAL::Polyhedron_3 to a
+// deal.II triangulation. The input mesh is a tripod.
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <CGAL/IO/io.h>
+#include <deal.II/cgal/triangulation.h>
+
+#include <fstream>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+using Kernel     = CGAL::Simple_cartesian<double>;
+using CGALPoint  = CGAL::Point_3<Kernel>;
+using Mesh       = CGAL::Surface_mesh<CGALPoint>;
+using Polyhedron = CGAL::Polyhedron_3<Kernel>;
+
+int
+main()
+{
+  initlog();
+  Triangulation<2, 3> tria;
+  GridOut             go;
+  {
+    deallog << "Surface mesh" << std::endl;
+    Mesh          sm;
+    std::ifstream input_sm(SOURCE_DIR "/input_grids/tripod.off");
+    input_sm >> sm;
+    cgal_surface_mesh_to_dealii_triangulation(sm, tria);
+    go.write_msh(tria, deallog.get_file_stream());
+  }
+  tria.clear();
+
+  // Now do the same with a Polyhedron
+  deallog << "Polyhedron test" << std::endl;
+  {
+    Polyhedron    P;
+    std::ifstream input_p(SOURCE_DIR "/input_grids/tripod.off");
+    input_p >> P;
+    cgal_surface_mesh_to_dealii_triangulation(P, tria);
+    go.write_msh(tria, deallog.get_file_stream());
+  }
+}
diff --git a/tests/cgal/cgal_surface_triangulation_03.output b/tests/cgal/cgal_surface_triangulation_03.output
new file mode 100644 (file)
index 0000000..e69de29
index beeb93604f2449eb1e12ac9892d6ceeeecf6782f..43c8201681c7ddab1a411efeafa3830ef2c4b569 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-// Convert a cgal Delaunay_triangulation_3 to a
-// dealii::Triangulation<dim, spacedim>
-// Non trivial case of a sphere. We need Exact predicates, inexact construction
-// kernel here, since Simple_cartesian will throw an exception when trying to
-// add some almost coplanar points.
+// Create a Triangulation<3,3> from an implicit function
 
 #include <deal.II/base/config.h>
 
@@ -38,14 +34,16 @@ main()
   // Build a deal.II triangulation
   Triangulation<3>  tria;
   FunctionParser<3> implicit_function("(1-sqrt(x^2+y^2))^2+z^2-.25");
+  CGALWrappers::AdditionalData<3> data;
+  data.cell_size = 0.4;
   GridGenerator::implicit_function(
-    tria, implicit_function, Point<3>(1, 0, 0), 10.0, cell_size = 0.4);
-  GridOut go;
+    tria, implicit_function, data, Point<3>(1, 0, 0), 10.0);
   {
+    GridOut       go;
     std::ofstream of("tria.vtk");
     go.write_vtk(tria, of);
   }
-  remove("tria.vtk");
-  // If we got here, everything was ok, including writing the grid.
+  // remove(/"tria.vtk");
+  //  If we got here, everything was ok, including writing the grid.
   deallog << "OK" << std::endl;
 }
diff --git a/tests/cgal/implicit_triangulation_02.cc b/tests/cgal/implicit_triangulation_02.cc
new file mode 100644 (file)
index 0000000..a0a76b0
--- /dev/null
@@ -0,0 +1,51 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Create a Triangulation<2,3> from an implicit function
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/function_parser.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+using namespace CGAL::parameters;
+
+int
+main()
+{
+  initlog();
+  // Build a deal.II triangulation
+  Triangulation<2, 3> tria;
+  FunctionParser<3>   implicit_function("(1-sqrt(x^2+y^2))^2+z^2-.25");
+  CGALWrappers::AdditionalData<2> data;
+  data.angular_bound  = 30.;
+  data.radius_bound   = .1;
+  data.distance_bound = .1;
+  GridGenerator::implicit_function(
+    tria, implicit_function, data, Point<3>(1, 0, 0), 10.);
+  {
+    GridOut       go;
+    std::ofstream of("tria_ancora.vtk");
+    go.write_vtk(tria, of);
+  }
+  remove("tria.vtk");
+  //  If we got here, everything was ok, including writing the grid.
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/cgal/implicit_triangulation_02.output b/tests/cgal/implicit_triangulation_02.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.