From: Marco Feder Date: Sat, 14 May 2022 18:07:50 +0000 (+0200) Subject: Convert a deal.II tria to a CGAL::Surface_mesh X-Git-Tag: v9.4.0-rc1~188^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1e023eeb634386a79ab7f49c65e08f193d38f074;p=dealii.git Convert a deal.II tria to a CGAL::Surface_mesh --- diff --git a/include/deal.II/cgal/surface_mesh.h b/include/deal.II/cgal/surface_mesh.h index e10fc76df7..818be5c069 100644 --- a/include/deal.II/cgal/surface_mesh.h +++ b/include/deal.II/cgal/surface_mesh.h @@ -19,14 +19,17 @@ #include #include +#include #include #include #ifdef DEAL_II_WITH_CGAL +# include # include + DEAL_II_NAMESPACE_OPEN namespace CGALWrappers @@ -66,6 +69,26 @@ namespace CGALWrappers const typename dealii::Triangulation::cell_iterator &cell, const dealii::Mapping & mapping, CGAL::Surface_mesh & mesh); + + /** + * Convert a deal.II triangulation to a CGAL::Surface_mesh. The output depends + * on the intrinsic dimension of the input deal.II triangulation. + * + * In 2D, i.e. with a + * Triangulation<2> or a Triangulation<2,3>, the output is the + * CGAL::Surface_mesh describing the whole triangulation. + * + * In 3D, the boundary the of the deal.II Triangulation is converted to + * a CGAL::Surface_mesh by looping over all the boundary faces. + * + * @param[in] triangulation The input deal.II triangulation. + * @param[out] mesh The output CGAL::Surface_mesh. + */ + template + void + dealii_tria_to_cgal_surface_mesh( + const dealii::Triangulation &triangulation, + CGAL::Surface_mesh & mesh); } // namespace CGALWrappers diff --git a/include/deal.II/cgal/utilities.h b/include/deal.II/cgal/utilities.h index 23e0efeda7..4c37536268 100644 --- a/include/deal.II/cgal/utilities.h +++ b/include/deal.II/cgal/utilities.h @@ -291,6 +291,10 @@ namespace CGALWrappers Assert(CGAL::is_closed(surface_mesh_2), ExcMessage( "The input surface_mesh_2 must be a closed surface mesh.")); + Assert(CGAL::is_triangle_mesh(surface_mesh_1), + ExcMessage("The first CGAL mesh must be triangulated.")); + Assert(CGAL::is_triangle_mesh(surface_mesh_2), + ExcMessage("The second CGAL mesh must be triangulated.")); bool res = false; auto surf_1 = surface_mesh_1; diff --git a/source/cgal/surface_mesh.cc b/source/cgal/surface_mesh.cc index 8d4fcf49f7..49aa0839b6 100644 --- a/source/cgal/surface_mesh.cc +++ b/source/cgal/surface_mesh.cc @@ -64,6 +64,23 @@ namespace "CGAL encountered a orientation problem that it " "was not able to solve.")); } + + + + template + void + map_vertices( + const dealiiFace & cell, + std::map &deal2cgal, + CGAL_Mesh & mesh) + { + for (const auto i : cell->vertex_indices()) + { + deal2cgal[cell->vertex_index(i)] = mesh.add_vertex( + CGALWrappers::dealii_point_to_cgal_point( + cell->vertex(i))); + } + } } // namespace @@ -80,10 +97,6 @@ namespace CGALWrappers CGAL::Surface_mesh & mesh) { Assert(dim > 1, ExcImpossibleInDim(dim)); - Assert( - mesh.is_empty(), - ExcMessage( - "The CGAL::Surface_mesh object must be empty upon calling this function.")); using Mesh = CGAL::Surface_mesh; const auto &vertices = mapping.get_vertices(cell); std::map deal2cgal; @@ -114,10 +127,57 @@ namespace CGALWrappers (f % 2 == 0)); } - // explicit instantiations -# include "surface_mesh.inst" + template + void + dealii_tria_to_cgal_surface_mesh( + const dealii::Triangulation &tria, + CGAL::Surface_mesh & mesh) + { + Assert(tria.n_cells() > 0, + ExcMessage( + "Triangulation cannot be empty upon calling this function.")); + Assert(mesh.is_empty(), + ExcMessage( + "The surface mesh must be empty upon calling this function.")); + + Assert(dim > 1, ExcImpossibleInDim(dim)); + using Mesh = CGAL::Surface_mesh; + using Vertex_index = typename Mesh::Vertex_index; + + std::map deal2cgal; + if constexpr (dim == 2) + { + for (const auto &cell : tria.active_cell_iterators()) + { + map_vertices(cell, deal2cgal, mesh); + add_facet(cell, deal2cgal, mesh); + } + } + else if constexpr (dim == 3 && spacedim == 3) + { + for (const auto &cell : tria.active_cell_iterators()) + { + for (const auto &f : cell->face_indices()) + + if (cell->face(f)->at_boundary()) + { + map_vertices(cell->face(f), deal2cgal, mesh); + add_facet(cell->face(f), + deal2cgal, + mesh, + (f % 2 == 0 || cell->n_vertices() != 8)); + } + } + } + else + { + Assert(false, ExcImpossibleInDimSpacedim(dim, spacedim)); + } + } // explicit instantiations +# include "surface_mesh.inst" + } // namespace CGALWrappers # endif diff --git a/source/cgal/surface_mesh.inst.in b/source/cgal/surface_mesh.inst.in index 985bd35eb9..6b84c7a9fd 100644 --- a/source/cgal/surface_mesh.inst.in +++ b/source/cgal/surface_mesh.inst.in @@ -25,5 +25,11 @@ for (dim : DIMENSIONS; spacedim : SPACE_DIMENSIONS; cgal_kernel : CGAL_KERNELS) const typename Triangulation::cell_iterator &cell, const Mapping & mapping, CGAL::Surface_mesh & mesh); + + template void dealii_tria_to_cgal_surface_mesh< + typename cgal_kernel::Point_3, + dim, + spacedim>(const Triangulation & cell, + CGAL::Surface_mesh &mesh); #endif } diff --git a/tests/cgal/cgal_surface_mesh_01.output b/tests/cgal/cgal_surface_mesh_01.output index df27e769bd..b2df3eeeac 100644 --- a/tests/cgal/cgal_surface_mesh_01.output +++ b/tests/cgal/cgal_surface_mesh_01.output @@ -110,4 +110,3 @@ DEAL::OFF 4 3 7 6 2 4 0 1 3 2 4 6 7 5 4 - diff --git a/tests/cgal/cgal_surface_mesh_04.cc b/tests/cgal/cgal_surface_mesh_04.cc new file mode 100644 index 0000000000..d603012725 --- /dev/null +++ b/tests/cgal/cgal_surface_mesh_04.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// 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 deal.II triangulation to a CGAL surface mesh. In the 2D case, +// thw whole triangulation is a 2D surface mesh. In 3D, the surface mesh +// describes the boundary of the deal.II Triangulation. + +#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 tria_in; + Triangulation<2, 3> tria_out; + GridOut go; + CGALMesh surface_mesh; + + std::vector> names_and_args; + if constexpr (spacedim == 2) + { + names_and_args = {{"hyper_cube", "0.0 : 1.0 : false"}, + {"hyper_ball", "0.,0. : 1. : false"}, + {"hyper_L", "0.0 : 1.0 : false"}, + {"channel_with_cylinder", "0.03 : 2 : 2.0 : false"}}; + } + else + { + names_and_args = {{"hyper_cube", "0.0 : 1.0 : false"}, + {"hyper_ball", "0.,0.,0. : 1. : false"}, + {"hyper_L", "0.0 : 1.0 : false"}, + {"channel_with_cylinder", "0.03 : 2 : 2.0 : false"}}; + } + + + for (const auto &info_pair : names_and_args) + { + auto name = info_pair.first; + auto args = info_pair.second; + deallog << "dim = " << dim << ", spacedim = " << spacedim + << " name: " << name << std::endl; + GridGenerator::generate_from_name_and_arguments(tria_in, name, args); + tria_in.refine_global(2); + dealii_tria_to_cgal_surface_mesh(tria_in, surface_mesh); + Assert(surface_mesh.is_valid(), + ExcMessage("The CGAL surface mesh is not valid.")); + + // Now back to the original dealii tria. + cgal_surface_mesh_to_dealii_triangulation(surface_mesh, tria_out); + std::ofstream out_name(name + std::to_string(spacedim) + ".vtk"); + go.write_vtk(tria_out, out_name); + // If we got here, everything was ok, including writing the grid. + deallog << "OK" << std::endl; + tria_in.clear(); + tria_out.clear(); + surface_mesh.clear(); + } +} + +int +main() +{ + initlog(); + test<2, 2>(); + test<3, 3>(); +} diff --git a/tests/cgal/cgal_surface_mesh_04.output b/tests/cgal/cgal_surface_mesh_04.output new file mode 100644 index 0000000000..e6f5341b18 --- /dev/null +++ b/tests/cgal/cgal_surface_mesh_04.output @@ -0,0 +1,19 @@ + +DEAL::dim= 2, spacedim= 2 +DEAL::dim = 2, spacedim = 2 name: hyper_cube +DEAL::OK +DEAL::dim = 2, spacedim = 2 name: hyper_ball +DEAL::OK +DEAL::dim = 2, spacedim = 2 name: hyper_L +DEAL::OK +DEAL::dim = 2, spacedim = 2 name: channel_with_cylinder +DEAL::OK +DEAL::dim= 3, spacedim= 3 +DEAL::dim = 3, spacedim = 3 name: hyper_cube +DEAL::OK +DEAL::dim = 3, spacedim = 3 name: hyper_ball +DEAL::OK +DEAL::dim = 3, spacedim = 3 name: hyper_L +DEAL::OK +DEAL::dim = 3, spacedim = 3 name: channel_with_cylinder +DEAL::OK diff --git a/tests/cgal/cgal_surface_mesh_05.cc b/tests/cgal/cgal_surface_mesh_05.cc new file mode 100644 index 0000000000..a542940233 --- /dev/null +++ b/tests/cgal/cgal_surface_mesh_05.cc @@ -0,0 +1,120 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Performs the following: deal.II Triangulations -> CGAL::Surface_mesh(es) -> +// Perform boolean operation -> deal.II tria. + +#include + +#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.5); + GridGenerator::hyper_ball(tria1, {0.2, 0.2, 0.2}, 0.5); + 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_intersection, + 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"); + go.write_vtk(tria_out, out_name_spheres); + deallog << "OK" << std::endl; + remove("tria_out.vtk"); + + // Clear everything + tria0.clear(); + tria1.clear(); + tria_out.clear(); + surface_mesh0.clear(); + surface_mesh1.clear(); + out_mesh.clear(); + + GridGenerator::hyper_cube(tria0); + GridGenerator::hyper_cube(tria1, 0.5, 1.5); + tria0.refine_global(3); + tria1.refine_global(3); + GridTools::rotate(numbers::PI_4, 2, tria1); + + // 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_intersection, + 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"); + go.write_vtk(tria_out, out_name_cubes); + remove("tria_out.vtk"); + deallog << "OK" << std::endl; +} + +int +main() +{ + initlog(); + test<3, 3>(); +} diff --git a/tests/cgal/cgal_surface_mesh_05.output b/tests/cgal/cgal_surface_mesh_05.output new file mode 100644 index 0000000000..b2c6daa84c --- /dev/null +++ b/tests/cgal/cgal_surface_mesh_05.output @@ -0,0 +1,4 @@ + +DEAL::dim= 3, spacedim= 3 +DEAL::OK +DEAL::OK