From: Luca Heltai Date: Sat, 14 May 2022 10:32:02 +0000 (+0300) Subject: CGAL Triangulation from implicit function. X-Git-Tag: v9.4.0-rc1~70^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=377f52fb8e9f4be7ddeb349aedf1b11305fb6516;p=dealii.git CGAL Triangulation from implicit function. --- diff --git a/doc/doxygen/images/grid_generator_implicit_function_3d.png b/doc/doxygen/images/grid_generator_implicit_function_3d.png new file mode 100644 index 0000000000..7154a04c00 Binary files /dev/null and b/doc/doxygen/images/grid_generator_implicit_function_3d.png differ diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 97365420b3..258e97b286 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -28,6 +28,17 @@ #include #include +#ifdef DEAL_II_WITH_CGAL +// All functions needed by the CGAL mesh generation utilities +# include +# include +# include +# include +# include +# include +# include +#endif + #include #include @@ -1747,6 +1758,78 @@ namespace GridGenerator Triangulation &tria, const std::string & grid_generator_function_name, const std::string & grid_generator_function_arguments); + +#ifdef DEAL_II_WITH_CGAL + /** + * 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, 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. + * + * 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. + * + * 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 + * 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; + * + * 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); + * @endcode + * + * The above snippet of code generates the following image: + * + * @image html grid_generator_implicit_function_3d.png + * + * + * @param[out] tria The output triangulation + * @param[in] implicit_function The implicit function + * @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 + void + implicit_function(Triangulation &tria, + const Function<3> & implicit_function, + const Point<3> & interior_point = Point<3>(), + const double & outer_ball_radius = 1.0, + Args... cgal_args); +#endif ///@} /** @@ -2655,6 +2738,68 @@ namespace GridGenerator const unsigned int, const double, const bool); + +# ifdef DEAL_II_WITH_CGAL + template + void + implicit_function(Triangulation &tria, + const Function<3> & implicit_function, + const Point<3> & interior_point, + const double & outer_ball_radius, + Args... cgal_args) + { + Assert(implicit_function.n_components == 1, + ExcMessage( + "The implicit function must have exactly one component.")); + Assert(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.")); + + 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; + + if constexpr (dim == 3) + { + using Tr = + CGAL::Mesh_triangulation_3::type; + using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3; + using Mesh_criteria = CGAL::Mesh_criteria_3; + + + auto cgal_implicit_function = [&](const Point_3 &p) { + return NumberType( + 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_args...); + auto cgal_triangulation = CGAL::make_mesh_3(domain, criteria); + CGALWrappers::cgal_triangulation_to_dealii_triangulation( + cgal_triangulation, tria); + } + else if constexpr (dim == 2) + {} + else + { + Assert(false, ExcImpossibleInDim(dim)); + } + } +# endif + #endif } // namespace GridGenerator diff --git a/tests/cgal/implicit_triangulation_01.cc b/tests/cgal/implicit_triangulation_01.cc new file mode 100644 index 0000000000..beeb93604f --- /dev/null +++ b/tests/cgal/implicit_triangulation_01.cc @@ -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. +// +// --------------------------------------------------------------------- + +// Convert a cgal Delaunay_triangulation_3 to a +// dealii::Triangulation +// 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. + +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace CGAL::parameters; + +int +main() +{ + initlog(); + // Build a deal.II triangulation + 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, cell_size = 0.4); + 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. + deallog << "OK" << std::endl; +} diff --git a/tests/cgal/implicit_triangulation_01.output b/tests/cgal/implicit_triangulation_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/cgal/implicit_triangulation_01.output @@ -0,0 +1,2 @@ + +DEAL::OK