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
utilities.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2022 - 2023 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_utilities_h
17#define dealii_cgal_utilities_h
18
19#include <deal.II/base/config.h>
20
22
24
25#ifdef DEAL_II_WITH_CGAL
27
28# include <deal.II/grid/tria.h>
29
30# include <boost/hana.hpp>
31
32# include <CGAL/Cartesian.h>
33# include <CGAL/Complex_2_in_triangulation_3.h>
34# include <CGAL/Exact_predicates_exact_constructions_kernel.h>
35# include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
36# include <CGAL/Kernel_traits.h>
37# include <CGAL/Mesh_complex_3_in_triangulation_3.h>
38# include <CGAL/Mesh_criteria_3.h>
39# include <CGAL/Mesh_triangulation_3.h>
40// Disable a warnung that we get with gcc-13 about a potential unitialized
41// usage of an <anonymous> lambda function in this external CGAL header.
43# include <CGAL/Polygon_mesh_processing/corefinement.h>
45# include <CGAL/Polygon_mesh_processing/measure.h>
46# include <CGAL/Polygon_mesh_processing/remesh.h>
47# include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
48# include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
49# include <CGAL/Simple_cartesian.h>
50# include <CGAL/Surface_mesh.h>
51# include <CGAL/Triangulation_3.h>
52# include <CGAL/boost/graph/copy_face_graph.h>
53# include <CGAL/boost/graph/selection.h>
54# include <CGAL/convex_hull_3.h>
55# include <CGAL/make_mesh_3.h>
56# include <CGAL/make_surface_mesh.h>
58
59# include <fstream>
60# include <limits>
61# include <type_traits>
62
63
64
78namespace CGALWrappers
79{
80# ifdef CGAL_CONCURRENT_MESH_3
81 using ConcurrencyTag = CGAL::Parallel_tag;
82# else
83 using ConcurrencyTag = CGAL::Sequential_tag;
84# endif
85
96 {
104 compute_corefinement = 1 << 0,
105
112 compute_difference = 1 << 1,
113
119 compute_intersection = 1 << 2,
120
126 compute_union = 1 << 3,
127 };
128
129
130
141 template <typename C3t3>
142 void
144 CGAL::Surface_mesh<typename C3t3::Point::Point> &surface_mesh,
145 C3t3 & triangulation,
146 const AdditionalData<3> &data = AdditionalData<3>{});
147
191 template <typename CGALPointType>
192 void
194 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_1,
195 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_2,
196 const BooleanOperation & boolean_operation,
197 CGAL::Surface_mesh<CGALPointType> & output_surface_mesh);
198
210 template <typename CGALTriangulationType>
212 compute_quadrature(const CGALTriangulationType &tria,
213 const unsigned int degree);
214
230 template <int dim0, int dim1, int spacedim>
233 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
234 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
235 const unsigned int degree,
236 const BooleanOperation & bool_op,
237 const Mapping<dim0, spacedim> &mapping0 =
238 (ReferenceCells::get_hypercube<dim0>()
239 .template get_default_linear_mapping<dim0, spacedim>()),
240 const Mapping<dim1, spacedim> &mapping1 =
241 (ReferenceCells::get_hypercube<dim1>()
242 .template get_default_linear_mapping<dim1, spacedim>()));
243
259 template <int dim0, int dim1, int spacedim>
262 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
263 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
264 const unsigned int degree,
265 const Mapping<dim0, spacedim> &mapping0,
266 const Mapping<dim1, spacedim> &mapping1);
267
296 template <typename CGALPointType>
297 void
298 remesh_surface(CGAL::Surface_mesh<CGALPointType> &surface_mesh,
299 const AdditionalData<3> & data = AdditionalData<3>{});
300} // namespace CGALWrappers
301
302# ifndef DOXYGEN
303// Template implementations
304namespace CGALWrappers
305{
306 template <typename C3t3>
307 void
309 CGAL::Surface_mesh<typename C3t3::Point::Point> &surface_mesh,
310 C3t3 & triangulation,
311 const AdditionalData<3> & data)
312 {
313 using CGALPointType = typename C3t3::Point::Point;
314 Assert(CGAL::is_closed(surface_mesh),
315 ExcMessage("The surface mesh must be closed."));
316
317 using K = typename CGAL::Kernel_traits<CGALPointType>::Kernel;
318 using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<
319 K,
320 CGAL::Surface_mesh<CGALPointType>>;
321 using Tr = typename CGAL::
322 Mesh_triangulation_3<Mesh_domain, CGAL::Default, ConcurrencyTag>::type;
323 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
324
325 CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh);
326 Mesh_domain domain(surface_mesh);
327 domain.detect_features();
328 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
329 CGAL::parameters::facet_angle = data.facet_angle,
330 CGAL::parameters::facet_distance =
331 data.facet_distance,
332 CGAL::parameters::cell_radius_edge_ratio =
333 data.cell_radius_edge_ratio,
334 CGAL::parameters::cell_size = data.cell_size);
335 // Mesh generation
336 triangulation = CGAL::make_mesh_3<C3t3>(domain,
337 criteria,
338 CGAL::parameters::no_perturb(),
339 CGAL::parameters::no_exude());
340 }
341
342
343
344 template <typename CGALPointType>
345 void
347 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_1,
348 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_2,
349 const BooleanOperation & boolean_operation,
350 CGAL::Surface_mesh<CGALPointType> & output_surface_mesh)
351 {
352 Assert(output_surface_mesh.is_empty(),
354 "output_surface_mesh must be empty upon calling this function"));
355 Assert(CGAL::is_closed(surface_mesh_1),
357 "The input surface_mesh_1 must be a closed surface mesh."));
358 Assert(CGAL::is_closed(surface_mesh_2),
360 "The input surface_mesh_2 must be a closed surface mesh."));
361 Assert(CGAL::is_triangle_mesh(surface_mesh_1),
362 ExcMessage("The first CGAL mesh must be triangulated."));
363 Assert(CGAL::is_triangle_mesh(surface_mesh_2),
364 ExcMessage("The second CGAL mesh must be triangulated."));
365
366 bool res = false;
367 auto surf_1 = surface_mesh_1;
368 auto surf_2 = surface_mesh_2;
369 namespace PMP = CGAL::Polygon_mesh_processing;
370 switch (boolean_operation)
371 {
373 res = PMP::corefine_and_compute_union(surf_1,
374 surf_2,
375 output_surface_mesh);
376 break;
378 res = PMP::corefine_and_compute_intersection(surf_1,
379 surf_2,
380 output_surface_mesh);
381 break;
383 res = PMP::corefine_and_compute_difference(surf_1,
384 surf_2,
385 output_surface_mesh);
386 break;
388 PMP::corefine(surf_1,
389 surf_2); // both surfaces are corefined
390 output_surface_mesh = std::move(surf_1);
391 res = true;
392 break;
393 default:
394 output_surface_mesh.clear();
395 break;
396 }
397 (void)res;
398 Assert(res,
399 ExcMessage("The boolean operation was not successfully computed."));
400 }
401
402
403
404 template <typename CGALTriangulationType>
406 compute_quadrature(const CGALTriangulationType &tria,
407 const unsigned int degree)
408 {
409 Assert(tria.is_valid(), ExcMessage("The triangulation is not valid."));
410 Assert(CGALTriangulationType::Point::Ambient_dimension::value == 3,
412 Assert(degree > 0,
413 ExcMessage("The degree of the Quadrature formula is not positive."));
414
415 constexpr int spacedim =
416 CGALTriangulationType::Point::Ambient_dimension::value;
417 std::vector<std::array<::Point<spacedim>, spacedim + 1>>
418 vec_of_simplices; // tets
419
420 std::array<::Point<spacedim>, spacedim + 1> simplex;
421 for (const auto &f : tria.finite_cell_handles())
422 {
423 for (unsigned int i = 0; i < (spacedim + 1); ++i)
424 {
425 simplex[i] =
426 cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
427 }
428
429 vec_of_simplices.push_back(simplex);
430 }
431
432 return QGaussSimplex<spacedim>(degree).mapped_quadrature(vec_of_simplices);
433 }
434
435
436
437 template <int dim0, int dim1, int spacedim>
440 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
441 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
442 const unsigned int degree,
443 const BooleanOperation & bool_op,
444 const Mapping<dim0, spacedim> &mapping0,
445 const Mapping<dim1, spacedim> &mapping1)
446 {
447 Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
448 ExcNotImplemented("2d geometries are not yet supported."));
449 if (dim1 > dim0)
450 {
452 cell1,
453 cell0,
454 degree,
455 bool_op,
456 mapping1,
457 mapping0); // This function works for dim1<=dim0, so swap them if this
458 // is not the case.
459 }
461 {
463 cell0, cell1, degree, mapping0, mapping1);
464 }
465 else
466 {
467 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
468 using CGALPoint = CGAL::Point_3<K>;
469 using CGALTriangulation = CGAL::Triangulation_3<K>;
470 CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
471 dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
472 dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
473 // They have to be triangle meshes
474 CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
475 CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
476 Assert(CGAL::is_triangle_mesh(surface_1) &&
477 CGAL::is_triangle_mesh(surface_2),
478 ExcMessage("The surface must be a triangle mesh."));
479 compute_boolean_operation(surface_1, surface_2, bool_op, out_surface);
480
481 CGALTriangulation tria;
482 tria.insert(out_surface.points().begin(), out_surface.points().end());
483 return compute_quadrature(tria, degree);
484 }
485 }
486
487
488
489 template <int dim0, int dim1, int spacedim>
492 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
493 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
494 const unsigned int degree,
495 const Mapping<dim0, spacedim> &mapping0,
496 const Mapping<dim1, spacedim> &mapping1)
497 {
498 Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
499 ExcNotImplemented("2d geometries are not yet supported."));
500 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
501 using CGALPoint = CGAL::Point_3<K>;
502 using CGALTriangulation = CGAL::Triangulation_3<K>;
503
504 CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
505 dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
506 dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
507 // They have to be triangle meshes
508 CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
509 CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
510
512 surface_2,
514 out_surface);
515 CGAL::Surface_mesh<CGALPoint> dummy;
516 CGALTriangulation tr;
517 CGAL::convex_hull_3(out_surface.points().begin(),
518 out_surface.points().end(),
519 dummy);
520 tr.insert(dummy.points().begin(), dummy.points().end());
521 return compute_quadrature(tr, degree);
522 }
523
524
525
526 template <typename CGALPointType>
527 void
528 remesh_surface(CGAL::Surface_mesh<CGALPointType> &cgal_mesh,
529 const AdditionalData<3> & data)
530 {
531 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
532 using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
533 // Polyhedron type
534 using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
535 // Triangulation
536 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
537 using C3t3 =
538 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
539 Mesh_domain::Corner_index,
540 Mesh_domain::Curve_index>;
541 // Criteria
542 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
543
544 Polyhedron poly;
545 CGAL::copy_face_graph(cgal_mesh, poly);
546 // Create a vector with only one element: the pointer to the
547 // polyhedron.
548 std::vector<Polyhedron *> poly_ptrs_vector(1, &poly);
549 // Create a polyhedral domain, with only one polyhedron,
550 // and no "bounding polyhedron", so the volumetric part of the
551 // domain will be empty.
552 Mesh_domain domain(poly_ptrs_vector.begin(), poly_ptrs_vector.end());
553 // Get sharp features
554 domain.detect_features(); // includes detection of borders
555
556 // Mesh criteria
557 const double default_value_edge_size = std::numeric_limits<double>::max();
558 if (data.edge_size > 0 &&
559 std::abs(data.edge_size - default_value_edge_size) > 1e-12)
560 {
561 Mesh_criteria criteria(CGAL::parameters::edge_size = data.edge_size,
562 CGAL::parameters::facet_size = data.facet_size,
563 CGAL::parameters::facet_angle = data.facet_angle,
564 CGAL::parameters::facet_distance =
565 data.facet_distance,
566 CGAL::parameters::cell_radius_edge_ratio =
567 data.cell_radius_edge_ratio,
568 CGAL::parameters::cell_size = data.cell_size);
569 // Mesh generation
570 C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
571 criteria,
572 CGAL::parameters::no_perturb(),
573 CGAL::parameters::no_exude());
574 cgal_mesh.clear();
575 CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
576 }
577 else if (std::abs(data.edge_size - default_value_edge_size) <= 1e-12)
578 {
579 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
580 CGAL::parameters::facet_angle = data.facet_angle,
581 CGAL::parameters::facet_distance =
582 data.facet_distance,
583 CGAL::parameters::cell_radius_edge_ratio =
584 data.cell_radius_edge_ratio,
585 CGAL::parameters::cell_size = data.cell_size);
586 // Mesh generation
587 C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
588 criteria,
589 CGAL::parameters::no_perturb(),
590 CGAL::parameters::no_exude());
591 cgal_mesh.clear();
592 CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
593 }
594 else
595 {
596 Assert(false, ExcInternalError());
597 }
598 }
599
600
601
608 template <int spacedim>
609 void
610 resort_dealii_vertices_to_cgal_order(const unsigned int structdim,
611 std::vector<Point<spacedim>> &vertices)
612 {
613 if (ReferenceCell::n_vertices_to_type(structdim, vertices.size()) ==
615 std::swap(vertices[2], vertices[3]);
616 }
617
618
619
627 template <int dim, int spacedim>
628 std::vector<Point<spacedim>>
629 get_vertices_in_cgal_order(
630 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
631 const Mapping<dim, spacedim> & mapping)
632 {
633 // Elements have to be rectangular or simplices
634 const unsigned int n_vertices = cell->n_vertices();
635 Assert((n_vertices == ReferenceCells::get_hypercube<dim>().n_vertices()) ||
636 (n_vertices == ReferenceCells::get_simplex<dim>().n_vertices()),
638
639 std::vector<Point<spacedim>> ordered_vertices(n_vertices);
640 std::copy_n(mapping.get_vertices(cell).begin(),
641 n_vertices,
642 ordered_vertices.begin());
643
644 resort_dealii_vertices_to_cgal_order(dim, ordered_vertices);
645
646 return ordered_vertices;
647 }
648
649
650
659 template <int dim, int spacedim>
660 std::vector<Point<spacedim>>
661 get_vertices_in_cgal_order(
662 const typename ::Triangulation<dim, spacedim>::cell_iterator &cell,
663 const unsigned int face_no,
664 const Mapping<dim, spacedim> & mapping)
665 {
666 // Elements have to be rectangular or simplices
667 const unsigned int n_vertices = cell->face(face_no)->n_vertices();
668 Assert(
669 (n_vertices == ReferenceCells::get_hypercube<dim - 1>().n_vertices()) ||
670 (n_vertices == ReferenceCells::get_simplex<dim - 1>().n_vertices()),
672
673 std::vector<Point<spacedim>> ordered_vertices(n_vertices);
674 std::copy_n(mapping.get_vertices(cell, face_no).begin(),
675 n_vertices,
676 ordered_vertices.begin());
677
678 resort_dealii_vertices_to_cgal_order(dim - 1, ordered_vertices);
679
680 return ordered_vertices;
681 }
682
683
684
685} // namespace CGALWrappers
686# endif
687
689
690#endif
691#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
Definition point.h:112
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 > > &simplices) const
static ReferenceCell n_vertices_to_type(const int dim, const unsigned int n_vertices)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition config.h:486
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition config.h:529
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void remesh_surface(CGAL::Surface_mesh< CGALPointType > &surface_mesh, const AdditionalData< 3 > &data=AdditionalData< 3 >{})
::Quadrature< CGALTriangulationType::Point::Ambient_dimension::value > compute_quadrature(const CGALTriangulationType &tria, const unsigned int degree)
CGAL::Sequential_tag ConcurrencyTag
Definition utilities.h:83
::Quadrature< spacedim > compute_quadrature_on_boolean_operation(const typename ::Triangulation< dim0, spacedim >::cell_iterator &cell0, const typename ::Triangulation< dim1, spacedim >::cell_iterator &cell1, const unsigned int degree, const BooleanOperation &bool_op, const Mapping< dim0, spacedim > &mapping0=(ReferenceCells::get_hypercube< dim0 >() .template get_default_linear_mapping< dim0, spacedim >()), const Mapping< dim1, spacedim > &mapping1=(ReferenceCells::get_hypercube< dim1 >() .template get_default_linear_mapping< dim1, spacedim >()))
void cgal_surface_mesh_to_cgal_triangulation(CGAL::Surface_mesh< typename C3t3::Point::Point > &surface_mesh, C3t3 &triangulation, const AdditionalData< 3 > &data=AdditionalData< 3 >{})
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
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)
void compute_boolean_operation(const CGAL::Surface_mesh< CGALPointType > &surface_mesh_1, const CGAL::Surface_mesh< CGALPointType > &surface_mesh_2, const BooleanOperation &boolean_operation, CGAL::Surface_mesh< CGALPointType > &output_surface_mesh)
::Quadrature< spacedim > compute_quadrature_on_intersection(const typename ::Triangulation< dim0, spacedim >::cell_iterator &cell0, const typename ::Triangulation< dim1, spacedim >::cell_iterator &cell1, const unsigned int degree, const Mapping< dim0, spacedim > &mapping0, const Mapping< dim1, spacedim > &mapping1)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
constexpr const ReferenceCell Quadrilateral
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria