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
triangulation.h
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#ifndef dealii_cgal_triangulation_h
17#define dealii_cgal_triangulation_h
18
19#include <deal.II/base/config.h>
20
23
24#include <deal.II/grid/tria.h>
25
27
28#ifdef DEAL_II_WITH_CGAL
29
30# include <boost/hana.hpp>
31
32# include <CGAL/Complex_2_in_triangulation_3.h>
33# include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
34# include <CGAL/Implicit_surface_3.h>
35# include <CGAL/Labeled_mesh_domain_3.h>
36# include <CGAL/Mesh_complex_3_in_triangulation_3.h>
37# include <CGAL/Mesh_criteria_3.h>
38# include <CGAL/Mesh_triangulation_3.h>
39# include <CGAL/Polyhedron_3.h>
40# include <CGAL/Surface_mesh.h>
41# include <CGAL/Surface_mesh_default_triangulation_3.h>
42# include <CGAL/Triangulation_2.h>
43# include <CGAL/Triangulation_3.h>
44# include <CGAL/make_mesh_3.h>
45# include <CGAL/make_surface_mesh.h>
47
49
50namespace CGALWrappers
51{
82 template <int spacedim, typename CGALTriangulation>
83 void
84 add_points_to_cgal_triangulation(const std::vector<Point<spacedim>> &points,
85 CGALTriangulation &triangulation);
86
118 template <typename CGALTriangulation, int dim, int spacedim>
119 void
120 cgal_triangulation_to_dealii_triangulation(
121 const CGALTriangulation & cgal_triangulation,
122 Triangulation<dim, spacedim> &dealii_triangulation);
123
147 template <typename CGALTriangulationType,
148 typename CornerIndexType,
149 typename CurveIndexType>
150 void
151 cgal_triangulation_to_dealii_triangulation(
152 const CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
153 CornerIndexType,
154 CurveIndexType>
155 & cgal_triangulation,
156 Triangulation<3> &dealii_triangulation);
157
179 template <typename CGAL_MeshType>
180 void
181 cgal_surface_mesh_to_dealii_triangulation(const CGAL_MeshType &cgal_mesh,
183
184
185
186# ifndef DOXYGEN
187 // Template implementation
188 template <int spacedim, typename CGALTriangulation>
189 void
190 add_points_to_cgal_triangulation(const std::vector<Point<spacedim>> &points,
191 CGALTriangulation &triangulation)
192 {
193 Assert(triangulation.is_valid(),
195 "The triangulation you pass to this function should be a valid "
196 "CGAL triangulation."));
197
198 std::vector<typename CGALTriangulation::Point> cgal_points(points.size());
199 std::transform(points.begin(),
200 points.end(),
201 cgal_points.begin(),
202 [](const auto &p) {
203 return CGALWrappers::dealii_point_to_cgal_point<
204 typename CGALTriangulation::Point>(p);
205 });
206
207 triangulation.insert(cgal_points.begin(), cgal_points.end());
208 Assert(triangulation.is_valid(),
210 "The Triangulation is no longer valid after inserting the points. "
211 "Bailing out."));
212 }
213
214
215
216 template <typename CGALTriangulationType,
217 typename CornerIndexType,
218 typename CurveIndexType>
219 void
220 cgal_triangulation_to_dealii_triangulation(
221 const CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
222 CornerIndexType,
223 CurveIndexType>
224 & cgal_triangulation,
225 Triangulation<3> &dealii_triangulation)
226 {
227 using C3T3 = CGAL::Mesh_complex_3_in_triangulation_3<CGALTriangulationType,
228 CornerIndexType,
229 CurveIndexType>;
230
231 // Extract all vertices first
232 std::vector<Point<3>> dealii_vertices;
233 std::map<typename C3T3::Vertex_handle, unsigned int>
234 cgal_to_dealii_vertex_map;
235
236 std::size_t inum = 0;
237 for (auto vit = cgal_triangulation.triangulation().finite_vertices_begin();
238 vit != cgal_triangulation.triangulation().finite_vertices_end();
239 ++vit)
240 {
241 if (vit->in_dimension() <= -1)
242 continue;
243 cgal_to_dealii_vertex_map[vit] = inum++;
244 dealii_vertices.emplace_back(
245 CGALWrappers::cgal_point_to_dealii_point<3>(vit->point()));
246 }
247
248 // Now build cell connectivity
249 std::vector<CellData<3>> cells;
250 for (auto cgal_cell = cgal_triangulation.cells_in_complex_begin();
251 cgal_cell != cgal_triangulation.cells_in_complex_end();
252 ++cgal_cell)
253 {
254 CellData<3> cell(ReferenceCells::Tetrahedron.n_vertices());
255 for (unsigned int i = 0; i < 4; ++i)
256 cell.vertices[i] = cgal_to_dealii_vertex_map[cgal_cell->vertex(i)];
257 cell.manifold_id = cgal_triangulation.subdomain_index(cgal_cell);
258 cells.push_back(cell);
259 }
260
261 // Do the same with surface patches, if possible
262 SubCellData subcell_data;
263 if constexpr (std::is_integral<typename C3T3::Surface_patch_index>::value)
264 {
265 for (auto face = cgal_triangulation.facets_in_complex_begin();
266 face != cgal_triangulation.facets_in_complex_end();
267 ++face)
268 {
269 const auto &[cgal_cell, cgal_vertex_face_index] = *face;
270 CellData<2> dealii_face(ReferenceCells::Triangle.n_vertices());
271 // A face is identified by a cell and by the index within the cell
272 // of the opposite vertex. Loop over vertices, and retain only those
273 // that belong to this face.
274 int j = 0;
275 for (int i = 0; i < 4; ++i)
276 if (i != cgal_vertex_face_index)
277 dealii_face.vertices[j++] =
278 cgal_to_dealii_vertex_map[cgal_cell->vertex(i)];
279 dealii_face.manifold_id =
280 cgal_triangulation.surface_patch_index(cgal_cell,
281 cgal_vertex_face_index);
282 subcell_data.boundary_quads.emplace_back(dealii_face);
283 }
284 }
285 // and curves
286 if constexpr (std::is_integral<typename C3T3::Curve_index>::value)
287 {
288 for (auto edge = cgal_triangulation.edges_in_complex_begin();
289 edge != cgal_triangulation.edges_in_complex_end();
290 ++edge)
291 {
292 const auto &[cgal_cell, v1, v2] = *edge;
293 CellData<1> dealii_edge(ReferenceCells::Line.n_vertices());
294 dealii_edge.vertices[0] =
295 cgal_to_dealii_vertex_map[cgal_cell->vertex(v1)];
296 dealii_edge.vertices[1] =
297 cgal_to_dealii_vertex_map[cgal_cell->vertex(v2)];
298 dealii_edge.manifold_id =
299 cgal_triangulation.curve_index(cgal_cell->vertex(v1),
300 cgal_cell->vertex(v2));
301 subcell_data.boundary_lines.emplace_back(dealii_edge);
302 }
303 }
304
305 dealii_triangulation.create_triangulation(dealii_vertices,
306 cells,
307 subcell_data);
308 }
309
310
311
312 template <typename CGALTriangulation, int dim, int spacedim>
313 void
314 cgal_triangulation_to_dealii_triangulation(
315 const CGALTriangulation & cgal_triangulation,
316 Triangulation<dim, spacedim> &dealii_triangulation)
317 {
318 AssertThrow(cgal_triangulation.dimension() == dim,
319 ExcMessage("The dimension of the input CGAL triangulation (" +
320 std::to_string(cgal_triangulation.dimension()) +
321 ") does not match the dimension of the output "
322 "deal.II triangulation (" +
323 std::to_string(dim) + ")."));
324
325 Assert(dealii_triangulation.n_cells() == 0,
326 ExcMessage("The output triangulation object needs to be empty."));
327
328 // deal.II storage data structures
329 std::vector<Point<spacedim>> vertices(
330 cgal_triangulation.number_of_vertices());
331 std::vector<CellData<dim>> cells;
332 SubCellData subcell_data;
333
334 // CGAL storage data structures
335 std::map<typename CGALTriangulation::Vertex_handle, unsigned int>
336 vertex_map;
337 {
338 unsigned int i = 0;
339 for (const auto &v : cgal_triangulation.finite_vertex_handles())
340 {
341 vertices[i] =
342 CGALWrappers::cgal_point_to_dealii_point<spacedim>(v->point());
343 vertex_map[v] = i++;
344 }
345 }
346
347 const auto has_faces = boost::hana::is_valid(
348 [](auto &&obj) -> decltype(obj.finite_face_handles()) {});
349
350 const auto has_cells = boost::hana::is_valid(
351 [](auto &&obj) -> decltype(obj.finite_cell_handles()) {});
352
353 // Different loops for Triangulation_2 and Triangulation_3 types.
354 if constexpr (decltype(has_faces(cgal_triangulation)){})
355 {
356 // This is a non-degenerate Triangulation_2
357 if (cgal_triangulation.dimension() == 2)
358 for (const auto &f : cgal_triangulation.finite_face_handles())
359 {
360 CellData<dim> cell(ReferenceCells::Triangle.n_vertices());
361 for (unsigned int i = 0;
363 ++i)
364 cell.vertices[i] = vertex_map[f->vertex(i)];
365 cells.push_back(cell);
366 }
367 else if (cgal_triangulation.dimension() == 1)
368 // This is a degenerate Triangulation_2, made of edges
369 for (const auto &e : cgal_triangulation.finite_edges())
370 {
371 // An edge is identified by a face and a vertex index in the
372 // face
373 const auto & f = e.first;
374 const auto & i = e.second;
375 CellData<dim> cell(ReferenceCells::Line.n_vertices());
376 unsigned int id = 0;
377 // Since an edge is identified by a face (a triangle) and the
378 // index of the opposite vertex to the edge, we can use this
379 // logic to infer the indices of the vertices of the edge: loop
380 // over all vertices, and keep only those that are not the
381 // opposite vertex of the edge.
382 for (unsigned int j = 0;
384 ++j)
385 if (j != i)
386 cell.vertices[id++] = vertex_map[f->vertex(j)];
387 cells.push_back(cell);
388 }
389 else
390 {
391 Assert(false, ExcInternalError());
392 }
393 }
394 else if constexpr (decltype(has_cells(cgal_triangulation)){})
395 {
396 // This is a non-degenerate Triangulation_3
397 if (cgal_triangulation.dimension() == 3)
398 for (const auto &c : cgal_triangulation.finite_cell_handles())
399 {
401 for (unsigned int i = 0;
403 ++i)
404 cell.vertices[i] = vertex_map[c->vertex(i)];
405 cells.push_back(cell);
406 }
407 else if (cgal_triangulation.dimension() == 2)
408 // This is a degenerate Triangulation_3, made of triangles
409 for (const auto &facet : cgal_triangulation.finite_facets())
410 {
411 // A facet is identified by a cell and the opposite vertex index
412 // in the face
413 const auto & c = facet.first;
414 const auto & i = facet.second;
415 CellData<dim> cell(ReferenceCells::Triangle.n_vertices());
416 unsigned int id = 0;
417 // Since a face is identified by a cell (a tetrahedron) and the
418 // index of the opposite vertex to the face, we can use this
419 // logic to infer the indices of the vertices of the face: loop
420 // over all vertices, and keep only those that are not the
421 // opposite vertex of the face.
422 for (unsigned int j = 0;
424 ++j)
425 if (j != i)
426 cell.vertices[id++] = vertex_map[c->vertex(j)];
427 cells.push_back(cell);
428 }
429 else if (cgal_triangulation.dimension() == 1)
430 // This is a degenerate Triangulation_3, made of edges
431 for (const auto &edge : cgal_triangulation.finite_edges())
432 {
433 // An edge is identified by a cell and its two vertices
434 const auto &[c, i, j] = edge;
435 CellData<dim> cell(ReferenceCells::Line.n_vertices());
436 cell.vertices[0] = vertex_map[c->vertex(i)];
437 cell.vertices[1] = vertex_map[c->vertex(j)];
438 cells.push_back(cell);
439 }
440 else
441 {
442 Assert(false, ExcInternalError());
443 }
444 }
445 dealii_triangulation.create_triangulation(vertices, cells, subcell_data);
446 }
447
448
449
450 template <typename CGAL_MeshType>
451 void
452 cgal_surface_mesh_to_dealii_triangulation(const CGAL_MeshType &cgal_mesh,
454 {
455 Assert(triangulation.n_cells() == 0,
457 "Triangulation must be empty upon calling this function."));
458
459 const auto is_surface_mesh =
460 boost::hana::is_valid([](auto &&obj) -> decltype(obj.faces()) {});
461
462 const auto is_polyhedral =
463 boost::hana::is_valid([](auto &&obj) -> decltype(obj.facets_begin()) {});
464
465 // Collect Vertices and cells
466 std::vector<::Point<3>> vertices;
467 std::vector<CellData<2>> cells;
468 SubCellData subcell_data;
469
470 // Different loops for Polyhedron or Surface_mesh types
471 if constexpr (decltype(is_surface_mesh(cgal_mesh)){})
472 {
473 AssertThrow(cgal_mesh.num_vertices() > 0,
474 ExcMessage("CGAL surface mesh is empty."));
475 vertices.reserve(cgal_mesh.num_vertices());
476 std::map<typename CGAL_MeshType::Vertex_index, unsigned int> vertex_map;
477 {
478 unsigned int i = 0;
479 for (const auto &v : cgal_mesh.vertices())
480 {
481 vertices.emplace_back(CGALWrappers::cgal_point_to_dealii_point<3>(
482 cgal_mesh.point(v)));
483 vertex_map[v] = i++;
484 }
485 }
486
487 // Collect CellData
488 for (const auto &face : cgal_mesh.faces())
489 {
490 const auto face_vertices =
491 CGAL::vertices_around_face(cgal_mesh.halfedge(face), cgal_mesh);
492
493 AssertThrow(face_vertices.size() == 3 || face_vertices.size() == 4,
494 ExcMessage("Only triangle or quadrilateral surface "
495 "meshes are supported in deal.II"));
496
497 CellData<2> c(face_vertices.size());
498 auto it_vertex = c.vertices.begin();
499 for (const auto &v : face_vertices)
500 *(it_vertex++) = vertex_map[v];
501
502 // Make sure the numberfing is consistent with the one in deal.II
503 if (face_vertices.size() == 4)
504 std::swap(c.vertices[3], c.vertices[2]);
505
506 cells.emplace_back(c);
507 }
508 }
509 else if constexpr (decltype(is_polyhedral(cgal_mesh)){})
510 {
511 AssertThrow(cgal_mesh.size_of_vertices() > 0,
512 ExcMessage("CGAL surface mesh is empty."));
513 vertices.reserve(cgal_mesh.size_of_vertices());
514 std::map<decltype(cgal_mesh.vertices_begin()), unsigned int> vertex_map;
515 {
516 unsigned int i = 0;
517 for (auto it = cgal_mesh.vertices_begin();
518 it != cgal_mesh.vertices_end();
519 ++it)
520 {
521 vertices.emplace_back(
522 CGALWrappers::cgal_point_to_dealii_point<3>(it->point()));
523 vertex_map[it] = i++;
524 }
525 }
526
527 // Loop over faces of Polyhedron, fill CellData
528 for (auto face = cgal_mesh.facets_begin();
529 face != cgal_mesh.facets_end();
530 ++face)
531 {
532 auto j = face->facet_begin();
533 const unsigned int vertices_per_face = CGAL::circulator_size(j);
534 AssertThrow(vertices_per_face == 3 || vertices_per_face == 4,
535 ExcMessage("Only triangle or quadrilateral surface "
536 "meshes are supported in deal.II. You "
537 "tried to read a mesh where a face has " +
538 std::to_string(vertices_per_face) +
539 " vertices per face."));
540
541 CellData<2> c(vertices_per_face);
542 auto it = c.vertices.begin();
543 for (unsigned int i = 0; i < vertices_per_face; ++i)
544 {
545 *(it++) = vertex_map[j->vertex()];
546 ++j;
547 }
548
549 if (vertices_per_face == 4)
550 std::swap(c.vertices[3], c.vertices[2]);
551
552 cells.emplace_back(c);
553 }
554 }
555 else
556 {
557 AssertThrow(false,
559 "Unsupported CGAL surface triangulation type."));
560 }
561 triangulation.create_triangulation(vertices, cells, subcell_data);
562 }
563} // namespace CGALWrappers
564# endif // doxygen
565
567
568#endif
569#endif
Definition point.h:112
unsigned int n_vertices() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
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]
const unsigned int v1
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Line
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< CellData< 2 > > boundary_quads
std::vector< CellData< 1 > > boundary_lines