From bd2bf3b033f8704d5e6ae67ed35b82b9d6ab04fb Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Wed, 8 Jun 2022 14:28:23 +0200 Subject: [PATCH] Move implementation of implicit_function() into .cc --- include/deal.II/cgal/additional_data.h | 136 +++++++++++++++++++++++++ include/deal.II/cgal/triangulation.h | 27 +++-- include/deal.II/cgal/utilities.h | 96 +---------------- include/deal.II/grid/grid_generator.h | 124 +--------------------- source/grid/grid_generator.cc | 121 ++++++++++++++++++++++ source/grid/grid_generator.inst.in | 8 ++ 6 files changed, 290 insertions(+), 222 deletions(-) create mode 100644 include/deal.II/cgal/additional_data.h diff --git a/include/deal.II/cgal/additional_data.h b/include/deal.II/cgal/additional_data.h new file mode 100644 index 0000000000..2fcdadfede --- /dev/null +++ b/include/deal.II/cgal/additional_data.h @@ -0,0 +1,136 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_cgal_additional_data_h +#define dealii_cgal_additional_data_h + +#include + +#ifdef DEAL_II_WITH_CGAL + +# include + +DEAL_II_NAMESPACE_OPEN +namespace CGALWrappers +{ + /** + * 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 + 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; + } + }; + +} // namespace CGALWrappers +DEAL_II_NAMESPACE_CLOSE + +#else + +DEAL_II_NAMESPACE_OPEN +namespace CGALWrappers +{ + template + struct AdditionalData + {}; +} // namespace CGALWrappers + +DEAL_II_NAMESPACE_CLOSE +#endif + +#endif diff --git a/include/deal.II/cgal/triangulation.h b/include/deal.II/cgal/triangulation.h index b6aa1d2e9a..3bc434f5de 100644 --- a/include/deal.II/cgal/triangulation.h +++ b/include/deal.II/cgal/triangulation.h @@ -18,17 +18,32 @@ #include -#ifdef DEAL_II_WITH_CGAL +#include +#include + +#include -# include +#include + +#ifdef DEAL_II_WITH_CGAL # include +# include +# include +# include +# include +# include +# include +# include # include # include +# include # include # include -# include +# include +# include +# include DEAL_II_NAMESPACE_OPEN @@ -256,8 +271,8 @@ namespace CGALWrappers // A face is identified by a cell and by the index within the cell // of the opposite vertex. Loop over vertices, and retain only those // that belong to this face. - unsigned int j = 0; - for (unsigned int i = 0; i < 4; ++i) + int j = 0; + for (int i = 0; i < 4; ++i) if (i != cgal_vertex_face_index) dealii_face.vertices[j++] = cgal_to_dealii_vertex_map[cgal_cell->vertex(i)]; @@ -545,8 +560,8 @@ namespace CGALWrappers } triangulation.create_triangulation(vertices, cells, subcell_data); } -# endif } // namespace CGALWrappers +# endif // doxygen DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 719a3e7d7d..56d0e69953 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -18,6 +18,8 @@ #include +#include + #ifdef DEAL_II_WITH_CGAL # include @@ -115,100 +117,6 @@ namespace CGALWrappers - /** - * 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 - 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; - } - }; - - - /** * Given a closed CGAL::Surface_mesh, this function fills the * region bounded by the surface with tets. diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 1035e1f092..054e28f5b7 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -28,22 +28,7 @@ #include #include -#ifdef DEAL_II_WITH_CGAL -// All functions needed by the CGAL mesh generation utilities -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -#endif +#include #include #include @@ -1765,7 +1750,6 @@ namespace GridGenerator 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. @@ -1823,7 +1807,6 @@ namespace GridGenerator CGALWrappers::AdditionalData{}, const Point<3> &interior_point = Point<3>(), const double & outer_ball_radius = 1.0); -#endif ///@} /** @@ -2733,110 +2716,7 @@ namespace GridGenerator const double, const bool); -# ifdef DEAL_II_WITH_CGAL - template - void - implicit_function(Triangulation &tria, - const Function<3> & dealii_implicit_function, - const CGALWrappers::AdditionalData &data, - const Point<3> & interior_point, - const double & outer_ball_radius) - { - 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; - 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( - 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(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; - 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; - using Surface_mesh = CGAL::Surface_mesh; - - - 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 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)); - } - } -# endif + #endif } // namespace GridGenerator diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 3c3535dbfc..e4a2afbcdc 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -31,6 +31,11 @@ #include +#ifdef DEAL_II_WITH_CGAL +// Functions needed by the CGAL mesh generation utilities are inside +# include +#endif + #include #include #include @@ -8639,6 +8644,122 @@ namespace GridGenerator } } + + + template + void + implicit_function(Triangulation &tria, + const Function<3> & dealii_implicit_function, + const CGALWrappers::AdditionalData &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; + 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( + 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(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; + 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; + using Surface_mesh = CGAL::Surface_mesh; + + + 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 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 diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 05d1fff1c2..0ef6acbddf 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -228,6 +228,14 @@ for (deal_II_dimension : DIMENSIONS) hyper_ball_balanced(Triangulation &, const Point &, const double); + + template void + implicit_function( + Triangulation &, + const Function<3> &, + const CGALWrappers::AdditionalData &, + const Point<3> &, + const double &); #endif \} } -- 2.39.5