#include <deal.II/base/config.h>
#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/tria.h>
#include <deal.II/cgal/utilities.h>
#ifdef DEAL_II_WITH_CGAL
+# include <CGAL/Polygon_mesh_processing/stitch_borders.h>
# include <CGAL/Surface_mesh.h>
+
DEAL_II_NAMESPACE_OPEN
namespace CGALWrappers
const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
const dealii::Mapping<dim, spacedim> & mapping,
CGAL::Surface_mesh<CGALPointType> & 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 <typename CGALPointType, int dim, int spacedim>
+ void
+ dealii_tria_to_cgal_surface_mesh(
+ const dealii::Triangulation<dim, spacedim> &triangulation,
+ CGAL::Surface_mesh<CGALPointType> & mesh);
} // 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;
"CGAL encountered a orientation problem that it "
"was not able to solve."));
}
+
+
+
+ template <typename dealiiFace, typename CGAL_Mesh>
+ void
+ map_vertices(
+ const dealiiFace & cell,
+ std::map<unsigned int, typename CGAL_Mesh::Vertex_index> &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<typename CGAL_Mesh::Point>(
+ cell->vertex(i)));
+ }
+ }
} // namespace
CGAL::Surface_mesh<CGALPointType> & 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<CGALPointType>;
const auto &vertices = mapping.get_vertices(cell);
std::map<unsigned int, typename Mesh::Vertex_index> deal2cgal;
(f % 2 == 0));
}
- // explicit instantiations
-# include "surface_mesh.inst"
+ template <typename CGALPointType, int dim, int spacedim>
+ void
+ dealii_tria_to_cgal_surface_mesh(
+ const dealii::Triangulation<dim, spacedim> &tria,
+ CGAL::Surface_mesh<CGALPointType> & 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<CGALPointType>;
+ using Vertex_index = typename Mesh::Vertex_index;
+
+ std::map<unsigned int, Vertex_index> 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
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Mapping<dim, spacedim> & mapping,
CGAL::Surface_mesh<typename cgal_kernel::Point_3> & mesh);
+
+ template void dealii_tria_to_cgal_surface_mesh<
+ typename cgal_kernel::Point_3,
+ dim,
+ spacedim>(const Triangulation<dim, spacedim> & cell,
+ CGAL::Surface_mesh<typename cgal_kernel::Point_3> &mesh);
#endif
}
4 3 7 6 2
4 0 1 3 2
4 6 7 5 4
-
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/surface_mesh.h>
+#include <deal.II/cgal/triangulation.h>
+#include <string.h>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_3<K>;
+using CGALMesh = CGAL::Surface_mesh<CGALPoint>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+ deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl;
+
+ Triangulation<spacedim> tria_in;
+ Triangulation<2, 3> tria_out;
+ GridOut go;
+ CGALMesh surface_mesh;
+
+ std::vector<std::pair<std::string, std::string>> 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>();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/surface_mesh.h>
+#include <deal.II/cgal/triangulation.h>
+#include <deal.II/cgal/utilities.h>
+#include <string.h>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_3<K>;
+using CGALMesh = CGAL::Surface_mesh<CGALPoint>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+ deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl;
+
+ Triangulation<spacedim> 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>();
+}
--- /dev/null
+
+DEAL::dim= 3, spacedim= 3
+DEAL::OK
+DEAL::OK
projects / dealii.git / commitdiff