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