Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
surface_mesh.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2022 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
18
19#ifdef DEAL_II_WITH_CGAL
20
21# include <CGAL/Exact_predicates_exact_constructions_kernel.h>
22# include <CGAL/Simple_cartesian.h>
23# include <CGAL/Surface_mesh/Surface_mesh.h>
24
25
27
28namespace
29{
30 template <typename dealiiFace, typename CGAL_Mesh>
31 void
32 add_facet(
33 const dealiiFace &face,
34 const std::map<unsigned int, typename CGAL_Mesh::Vertex_index> &deal2cgal,
35 CGAL_Mesh &mesh,
36 const bool clockwise_ordering = true)
37 {
38 const auto reference_cell_type = face->reference_cell();
39 std::vector<typename CGAL_Mesh::Vertex_index> indices;
40
41 if (reference_cell_type == ReferenceCells::Line)
42 {
43 mesh.add_edge(deal2cgal.at(face->vertex_index(0)),
44 deal2cgal.at(face->vertex_index(1)));
45 }
46 else if (reference_cell_type == ReferenceCells::Triangle)
47 {
48 indices = {deal2cgal.at(face->vertex_index(0)),
49 deal2cgal.at(face->vertex_index(1)),
50 deal2cgal.at(face->vertex_index(2))};
51 }
52
53 else if (reference_cell_type == ReferenceCells::Quadrilateral)
54 {
55 indices = {deal2cgal.at(face->vertex_index(0)),
56 deal2cgal.at(face->vertex_index(1)),
57 deal2cgal.at(face->vertex_index(3)),
58 deal2cgal.at(face->vertex_index(2))};
59 }
60 else
62
63 if (clockwise_ordering == true)
64 std::reverse(indices.begin(), indices.end());
65
66 [[maybe_unused]] const auto new_face = mesh.add_face(indices);
67 Assert(new_face != mesh.null_face(),
68 ExcInternalError("While trying to build a CGAL facet, "
69 "CGAL encountered a orientation problem that it "
70 "was not able to solve."));
71 }
72
73
74
75 template <typename dealiiFace, typename CGAL_Mesh>
76 void
77 map_vertices(
78 const dealiiFace &cell,
79 std::map<unsigned int, typename CGAL_Mesh::Vertex_index> &deal2cgal,
80 CGAL_Mesh &mesh)
81 {
82 for (const auto i : cell->vertex_indices())
83 {
84 deal2cgal[cell->vertex_index(i)] = mesh.add_vertex(
85 CGALWrappers::dealii_point_to_cgal_point<typename CGAL_Mesh::Point>(
86 cell->vertex(i)));
87 }
88 }
89} // namespace
90
91
92
93# ifndef DOXYGEN
94// Template implementations
95namespace CGALWrappers
96{
97 template <typename CGALPointType, int dim, int spacedim>
98 void
101 const Mapping<dim, spacedim> &mapping,
102 CGAL::Surface_mesh<CGALPointType> &mesh)
103 {
104 Assert(dim > 1, ExcImpossibleInDim(dim));
105 using Mesh = CGAL::Surface_mesh<CGALPointType>;
106 const auto &vertices = mapping.get_vertices(cell);
107 std::map<unsigned int, typename Mesh::Vertex_index> deal2cgal;
108
109 // Add all vertices to the mesh
110 // Store CGAL ordering
111 for (const auto &i : cell->vertex_indices())
112 deal2cgal[cell->vertex_index(i)] = mesh.add_vertex(
113 CGALWrappers::dealii_point_to_cgal_point<CGALPointType>(vertices[i]));
114
115 // Add faces
116 if (dim < 3)
117 // simplices and quads are allowable faces for CGAL
118 add_facet(cell, deal2cgal, mesh);
119 else
120 // in 3d, we build a surface mesh containing all the faces of the 3d cell.
121 // Simplices, Tetrahedrons, and Pyramids have their faces numbered in the
122 // same way as CGAL does (all faces of a bounding polyhedron are numbered
123 // counter-clockwise, so that their normal points outwards). Hexahedrons
124 // in deal.II have their faces numbered lexicographically, and one cannot
125 // deduce the direction of the normals by just looking at the vertices.
126 //
127 // In order for CGAL to be able to produce the right orientation, we need
128 // to reverse the order of the vertices for faces with even index.
129 // However, in order to allow for all kinds of meshes in 3d, including
130 // Moebius-loops, a deal.II face might even be rotated looking from one
131 // cell, whereas it is according to the standard when looking at it from
132 // the neighboring cell sharing that particular face. Therefore, when
133 // building a cgal face we must also take into account the fact that a
134 // face may have a non-standard orientation.
135 for (const auto &f : cell->face_indices())
136 {
137 // Check for standard orientation of faces
138 bool face_is_clockwise_oriented =
139 cell->reference_cell() != ReferenceCells::Hexahedron ||
140 (f % 2 == 0);
141 // Make sure that we revert the orientation if required
142 if (cell->face_orientation(f) == false)
143 face_is_clockwise_oriented = !face_is_clockwise_oriented;
144 add_facet(cell->face(f), deal2cgal, mesh, face_is_clockwise_oriented);
145 }
146 }
147
148
149
150 template <typename CGALPointType, int dim, int spacedim>
151 void
153 const ::Triangulation<dim, spacedim> &tria,
154 CGAL::Surface_mesh<CGALPointType> &mesh)
155 {
156 Assert(tria.n_cells() > 0,
158 "Triangulation cannot be empty upon calling this function."));
159 Assert(mesh.is_empty(),
161 "The surface mesh must be empty upon calling this function."));
162
163 Assert(dim > 1, ExcImpossibleInDim(dim));
164 using Mesh = CGAL::Surface_mesh<CGALPointType>;
165 using Vertex_index = typename Mesh::Vertex_index;
166
167 std::map<unsigned int, Vertex_index> deal2cgal;
168 if constexpr (dim == 2)
169 {
170 for (const auto &cell : tria.active_cell_iterators())
171 {
172 map_vertices(cell, deal2cgal, mesh);
173 add_facet(cell, deal2cgal, mesh);
174 }
175 }
176 else if constexpr (dim == 3 && spacedim == 3)
177 {
178 for (const auto &cell : tria.active_cell_iterators())
179 {
180 for (const auto &f : cell->face_indices())
181
182 if (cell->face(f)->at_boundary())
183 {
184 map_vertices(cell->face(f), deal2cgal, mesh);
185 add_facet(cell->face(f),
186 deal2cgal,
187 mesh,
188 (f % 2 == 0 || cell->n_vertices() != 8));
189 }
190 }
191 }
192 else
193 {
194 Assert(false, ExcImpossibleInDimSpacedim(dim, spacedim));
195 }
196 }
197
198// Explicit instantiations.
199//
200// We don't build the instantiations.inst file if deal.II isn't
201// configured with CGAL, but doxygen doesn't know that and tries to
202// find that file anyway for parsing -- which then of course it fails
203// on. So exclude the following from doxygen consideration.
204# ifndef DOXYGEN
205# include "cgal/surface_mesh.inst"
206# endif
207
208} // namespace CGALWrappers
209# endif
210
211
213
214#endif
Abstract base class for mapping classes.
Definition mapping.h:320
virtual boost::container::small_vector< Point< spacedim >, ReferenceCells::max_n_vertices< dim >() > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
unsigned int vertex_indices[2]
static ::ExceptionBase & ExcImpossibleInDimSpacedim(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void dealii_tria_to_cgal_surface_mesh(const ::Triangulation< dim, spacedim > &triangulation, CGAL::Surface_mesh< CGALPointType > &mesh)
void dealii_cell_to_cgal_surface_mesh(const typename ::Triangulation< dim, spacedim >::cell_iterator &cell, const ::Mapping< dim, spacedim > &mapping, CGAL::Surface_mesh< CGALPointType > &mesh)
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Line