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