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