Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
22
23namespace
24{
25 template <typename dealiiFace, typename CGAL_Mesh>
26 void
27 add_facet(
28 const dealiiFace &face,
29 const std::map<unsigned int, typename CGAL_Mesh::Vertex_index> &deal2cgal,
30 CGAL_Mesh &mesh,
31 const bool clockwise_ordering = true)
32 {
33 const auto reference_cell_type = face->reference_cell();
34 std::vector<typename CGAL_Mesh::Vertex_index> indices;
35
36 if (reference_cell_type == ReferenceCells::Line)
37 {
38 mesh.add_edge(deal2cgal.at(face->vertex_index(0)),
39 deal2cgal.at(face->vertex_index(1)));
40 }
41 else if (reference_cell_type == ReferenceCells::Triangle)
42 {
43 indices = {deal2cgal.at(face->vertex_index(0)),
44 deal2cgal.at(face->vertex_index(1)),
45 deal2cgal.at(face->vertex_index(2))};
46 }
47
48 else if (reference_cell_type == ReferenceCells::Quadrilateral)
49 {
50 indices = {deal2cgal.at(face->vertex_index(0)),
51 deal2cgal.at(face->vertex_index(1)),
52 deal2cgal.at(face->vertex_index(3)),
53 deal2cgal.at(face->vertex_index(2))};
54 }
55 else
57
58 if (clockwise_ordering == true)
59 std::reverse(indices.begin(), indices.end());
60
61 [[maybe_unused]] const auto new_face = mesh.add_face(indices);
62 Assert(new_face != mesh.null_face(),
63 ExcInternalError("While trying to build a CGAL facet, "
64 "CGAL encountered a orientation problem that it "
65 "was not able to solve."));
66 }
67
68
69
70 template <typename dealiiFace, typename CGAL_Mesh>
71 void
72 map_vertices(
73 const dealiiFace &cell,
74 std::map<unsigned int, typename CGAL_Mesh::Vertex_index> &deal2cgal,
75 CGAL_Mesh &mesh)
76 {
77 for (const auto i : cell->vertex_indices())
78 {
79 deal2cgal[cell->vertex_index(i)] = mesh.add_vertex(
80 CGALWrappers::dealii_point_to_cgal_point<typename CGAL_Mesh::Point>(
81 cell->vertex(i)));
82 }
83 }
84} // namespace
85
86
87
88# ifndef DOXYGEN
89// Template implementations
90namespace CGALWrappers
91{
92 template <typename CGALPointType, int dim, int spacedim>
93 void
96 const Mapping<dim, spacedim> &mapping,
97 CGAL::Surface_mesh<CGALPointType> &mesh)
98 {
99 Assert(dim > 1, ExcImpossibleInDim(dim));
100 using Mesh = CGAL::Surface_mesh<CGALPointType>;
101 const auto &vertices = mapping.get_vertices(cell);
102 std::map<unsigned int, typename Mesh::Vertex_index> deal2cgal;
103
104 // Add all vertices to the mesh
105 // Store CGAL ordering
106 for (const auto &i : cell->vertex_indices())
107 deal2cgal[cell->vertex_index(i)] = mesh.add_vertex(
108 CGALWrappers::dealii_point_to_cgal_point<CGALPointType>(vertices[i]));
109
110 // Add faces
111 if (dim < 3)
112 // simplices and quads are allowable faces for CGAL
113 add_facet(cell, deal2cgal, mesh);
114 else
115 // in 3d, we build a surface mesh containing all the faces of the 3d cell.
116 // Simplices, Tetrahedrons, and Pyramids have their faces numbered in the
117 // same way as CGAL does (all faces of a bounding polyhedron are numbered
118 // counter-clockwise, so that their normal points outwards). Hexahedrons
119 // in deal.II have their faces numbered lexicographically, and one cannot
120 // deduce the direction of the normals by just looking at the vertices.
121 //
122 // In order for CGAL to be able to produce the right orientation, we need
123 // to reverse the order of the vertices for faces with even index.
124 // However, in order to allow for all kinds of meshes in 3d, including
125 // Moebius-loops, a deal.II face might even be rotated looking from one
126 // cell, whereas it is according to the standard when looking at it from
127 // the neighboring cell sharing that particular face. Therefore, when
128 // building a cgal face we must also take into account the fact that a
129 // face may have a non-standard orientation.
130 for (const auto &f : cell->face_indices())
131 {
132 // Check for standard orientation of faces
133 bool face_is_clockwise_oriented =
134 cell->reference_cell() != ReferenceCells::Hexahedron ||
135 (f % 2 == 0);
136 // Make sure that we revert the orientation if required
137 if (cell->face_orientation(f) == false)
138 face_is_clockwise_oriented = !face_is_clockwise_oriented;
139 add_facet(cell->face(f), deal2cgal, mesh, face_is_clockwise_oriented);
140 }
141 }
142
143
144
145 template <typename CGALPointType, int dim, int spacedim>
146 void
148 const ::Triangulation<dim, spacedim> &tria,
149 CGAL::Surface_mesh<CGALPointType> &mesh)
150 {
151 Assert(tria.n_cells() > 0,
153 "Triangulation cannot be empty upon calling this function."));
154 Assert(mesh.is_empty(),
156 "The surface mesh must be empty upon calling this function."));
157
158 Assert(dim > 1, ExcImpossibleInDim(dim));
159 using Mesh = CGAL::Surface_mesh<CGALPointType>;
160 using Vertex_index = typename Mesh::Vertex_index;
161
162 std::map<unsigned int, Vertex_index> deal2cgal;
163 if constexpr (dim == 2)
164 {
165 for (const auto &cell : tria.active_cell_iterators())
166 {
167 map_vertices(cell, deal2cgal, mesh);
168 add_facet(cell, deal2cgal, mesh);
169 }
170 }
171 else if constexpr (dim == 3 && spacedim == 3)
172 {
173 for (const auto &cell : tria.active_cell_iterators())
174 {
175 for (const auto &f : cell->face_indices())
176
177 if (cell->face(f)->at_boundary())
178 {
179 map_vertices(cell->face(f), deal2cgal, mesh);
180 add_facet(cell->face(f),
181 deal2cgal,
182 mesh,
183 (f % 2 == 0 || cell->n_vertices() != 8));
184 }
185 }
186 }
187 else
188 {
189 Assert(false, ExcImpossibleInDimSpacedim(dim, spacedim));
190 }
191 } // explicit instantiations
192# include "surface_mesh.inst"
193
194} // namespace CGALWrappers
195# endif
196
197
199
200#endif
Abstract base class for mapping classes.
Definition mapping.h:318
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
Point< 3 > vertices[4]
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)
#define DEAL_II_ASSERT_UNREACHABLE()
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 const ReferenceCell Quadrilateral
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Line