From: Marco Feder Date: Wed, 8 Jun 2022 08:12:25 +0000 (+0200) Subject: Add remesh_surface utility X-Git-Tag: v9.4.0-rc1~29^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F13771%2Fhead;p=dealii.git Add remesh_surface utility --- diff --git a/doc/doxygen/images/boolean_union_hyper_spheres.png b/doc/doxygen/images/boolean_union_hyper_spheres.png new file mode 100644 index 0000000000..c3c4979619 Binary files /dev/null and b/doc/doxygen/images/boolean_union_hyper_spheres.png differ diff --git a/doc/doxygen/images/boolean_union_hyper_spheres_remeshed.png b/doc/doxygen/images/boolean_union_hyper_spheres_remeshed.png new file mode 100644 index 0000000000..c802fa8b4c Binary files /dev/null and b/doc/doxygen/images/boolean_union_hyper_spheres_remeshed.png differ diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index f3627efbc3..44e46bc511 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -39,16 +39,21 @@ # include # include # include +# include # include # include # include # include # include +# include +# include # include # include # include # include +//# include REQUIRES CGAL_VERSION>=5.1.5 +# include # include @@ -255,6 +260,39 @@ namespace CGALWrappers const unsigned int degree, const Mapping &mapping0, const Mapping &mapping1); + + /** + * Remesh a CGAL::Surface_mesh. + * + * If the domain has 1-dimensional exposed features, the criteria includes a + * sizing field to guide the sampling of 1-dimensional features with + * protecting balls centers. + * - CGAL::parameters::edge_size. + * + * This utility should be exploited to improve the quality of a mesh coming + * from boolean operations. As an example, consider two CGAL::Surface_mesh + * objects describing two hyper_spheres that intersect non-trivially. After + * computing their corefinement and union with compute_boolean_operation(), + * the resulting mesh is the following: + * + * @image html boolean_union_hyper_spheres.png + * + * Clearly, elements (triangles) like the ones arising along the intersection + * of the two spheres should be avoided as they're quite degenerate and would + * decrease the accuracy of the results. A call to the present function will + * try to remesh the surface, according to the given named parameters, to + * improve its quality. In the present case, the result is the following: + * + * @image html boolean_union_hyper_spheres_remeshed.png + * + * @param surface_mesh The input CGAL::Surface_mesh. + * @param [in] data AdditionalData object to pass to the CGAL::make_mesh_3 function. See the documentation + * of that struct for a description of those parameters. + */ + template + void + remesh_surface(CGAL::Surface_mesh &surface_mesh, + const AdditionalData<3> & data = AdditionalData<3>{}); } // namespace CGALWrappers # ifndef DOXYGEN @@ -528,6 +566,82 @@ namespace CGALWrappers tr.insert(dummy.points().begin(), dummy.points().end()); return compute_quadrature(tr, degree); } + + + + template + void + remesh_surface(CGAL::Surface_mesh &cgal_mesh, + const AdditionalData<3> & data) + { + using K = CGAL::Exact_predicates_inexact_constructions_kernel; + using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3; + // Polyhedron type + using Polyhedron = CGAL::Mesh_polyhedron_3::type; + // Triangulation + using Tr = CGAL::Mesh_triangulation_3::type; + using C3t3 = + CGAL::Mesh_complex_3_in_triangulation_3; + // Criteria + using Mesh_criteria = CGAL::Mesh_criteria_3; + + Polyhedron poly; + CGAL::copy_face_graph(cgal_mesh, poly); + // Create a vector with only one element: the pointer to the + // polyhedron. + std::vector poly_ptrs_vector(1, &poly); + // Create a polyhedral domain, with only one polyhedron, + // and no "bounding polyhedron", so the volumetric part of the + // domain will be empty. + Mesh_domain domain(poly_ptrs_vector.begin(), poly_ptrs_vector.end()); + // Get sharp features + domain.detect_features(); // includes detection of borders + + // Mesh criteria + const double default_value_edge_size = std::numeric_limits::max(); + if (data.edge_size > 0 && + std::abs(data.edge_size - default_value_edge_size) > 1e-12) + { + Mesh_criteria criteria(CGAL::parameters::edge_size = data.edge_size, + 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); + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, + criteria, + CGAL::parameters::no_perturb(), + CGAL::parameters::no_exude()); + cgal_mesh.clear(); + CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh); + } + else if (std::abs(data.edge_size - default_value_edge_size) <= 1e-12) + { + 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); + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, + criteria, + CGAL::parameters::no_perturb(), + CGAL::parameters::no_exude()); + cgal_mesh.clear(); + CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh); + } + else + { + Assert(false, ExcInternalError()); + } + } } // namespace CGALWrappers # endif diff --git a/tests/cgal/cgal_remesh_surface.cc b/tests/cgal/cgal_remesh_surface.cc new file mode 100644 index 0000000000..6a1a8426dc --- /dev/null +++ b/tests/cgal/cgal_remesh_surface.cc @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Remesh the union of two deal.II hyper_spheres in order to avoid bad +// triangles. + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace CGALWrappers; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using CGALPoint = CGAL::Point_3; +using CGALMesh = CGAL::Surface_mesh; + +template +void +test() +{ + deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl; + + Triangulation tria0, tria1; + Triangulation<2, 3> tria_out; + GridOut go; + CGALMesh surface_mesh0, surface_mesh1, out_mesh; + + GridGenerator::hyper_ball(tria0, {0., 0., 0.}, 0.4); + GridGenerator::hyper_ball(tria1, {0.3, 0.3, 0.3}, 0.4); + tria0.refine_global(3); + tria1.refine_global(3); + + // Move to CGAL surfaces + dealii_tria_to_cgal_surface_mesh(tria0, surface_mesh0); + dealii_tria_to_cgal_surface_mesh(tria1, surface_mesh1); + + // close the surfaces + CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh0); + CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh1); + + CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0); + CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1); + + compute_boolean_operation(surface_mesh0, + surface_mesh1, + BooleanOperation::compute_union, + out_mesh); + AdditionalData<3> data; + data.edge_size = .02; + data.facet_angle = 25; + data.facet_size = 0.05; + data.facet_distance = 0.05; + + remesh_surface(out_mesh, data); + + // Now back to deal.II + cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out); + std::ofstream out_name_spheres("remeshed_union_spheres.vtk"); + go.write_vtk(tria_out, out_name_spheres); + deallog << "OK" << std::endl; +} + +int +main() +{ + initlog(); + test<3, 3>(); +} diff --git a/tests/cgal/cgal_remesh_surface.output b/tests/cgal/cgal_remesh_surface.output new file mode 100644 index 0000000000..3e5e6a7d8a --- /dev/null +++ b/tests/cgal/cgal_remesh_surface.output @@ -0,0 +1,3 @@ + +DEAL::dim= 3, spacedim= 3 +DEAL::OK diff --git a/tests/cgal/cgal_surface_mesh_05.cc b/tests/cgal/cgal_surface_mesh_05.cc index a542940233..371ca2c2e8 100644 --- a/tests/cgal/cgal_surface_mesh_05.cc +++ b/tests/cgal/cgal_surface_mesh_05.cc @@ -70,10 +70,10 @@ test() out_mesh); // Now back to deal.II cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out); - std::ofstream out_name_spheres("boolean_triangulation_spheres.vtk"); + std::ofstream out_name_spheres("boolean_intersection_hyper_spheres.vtk"); go.write_vtk(tria_out, out_name_spheres); deallog << "OK" << std::endl; - remove("tria_out.vtk"); + remove("boolean_intersection_hyper_spheres.vtk"); // Clear everything tria0.clear(); @@ -106,10 +106,10 @@ test() out_mesh); // Now back to deal.II cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out); - std::ofstream out_name_cubes("boolean_triangulation_cubes.vtk"); + std::ofstream out_name_cubes("boolean_intersection_cubes.vtk"); go.write_vtk(tria_out, out_name_cubes); - remove("tria_out.vtk"); deallog << "OK" << std::endl; + remove("boolean_intersection_cubes.vtk"); } int