From: Luca Heltai Date: Tue, 17 May 2022 16:29:09 +0000 (+0300) Subject: Added Implementation for Triangulation<2,3> X-Git-Tag: v9.4.0-rc1~70^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a6141d74f7a134752523270bcc2905617bc9369f;p=dealii.git Added Implementation for Triangulation<2,3> --- 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 index 0000000000..a61a23e8bc Binary files /dev/null and b/doc/doxygen/images/grid_generator_implicit_function_2d.png differ diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 5758429575..ae8b7178f9 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -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 + 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. * diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 258e97b286..1035e1f092 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -30,11 +30,17 @@ #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 @@ -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 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 + 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); + implicit_function(Triangulation & tria, + const Function<3> & implicit_function, + const CGALWrappers::AdditionalData &data = + CGALWrappers::AdditionalData{}, + 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 + template void implicit_function(Triangulation &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 &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; - 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(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(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)); diff --git a/tests/cgal/cgal_surface_triangulation_03.cc b/tests/cgal/cgal_surface_triangulation_03.cc new file mode 100644 index 0000000000..5edd6e4ae7 --- /dev/null +++ b/tests/cgal/cgal_surface_triangulation_03.cc @@ -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 + +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace CGALWrappers; +using Kernel = CGAL::Simple_cartesian; +using CGALPoint = CGAL::Point_3; +using Mesh = CGAL::Surface_mesh; +using Polyhedron = CGAL::Polyhedron_3; + +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 index 0000000000..e69de29bb2 diff --git a/tests/cgal/implicit_triangulation_01.cc b/tests/cgal/implicit_triangulation_01.cc index beeb93604f..43c8201681 100644 --- a/tests/cgal/implicit_triangulation_01.cc +++ b/tests/cgal/implicit_triangulation_01.cc @@ -13,11 +13,7 @@ // // --------------------------------------------------------------------- -// 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. +// Create a Triangulation<3,3> from an implicit function #include @@ -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 index 0000000000..a0a76b0008 --- /dev/null +++ b/tests/cgal/implicit_triangulation_02.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. +// +// --------------------------------------------------------------------- + +// Create a Triangulation<2,3> from an implicit function + +#include + +#include + +#include +#include +#include + +#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 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/cgal/implicit_triangulation_02.output @@ -0,0 +1,2 @@ + +DEAL::OK