Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
utilities.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_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# include <CGAL/Polygon_mesh_processing/corefinement.h>
41# include <CGAL/Polygon_mesh_processing/measure.h>
42# include <CGAL/Polygon_mesh_processing/remesh.h>
43# include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
44# include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
45# include <CGAL/Simple_cartesian.h>
46# include <CGAL/Surface_mesh.h>
47# include <CGAL/Triangulation_3.h>
48# include <CGAL/boost/graph/copy_face_graph.h>
49# include <CGAL/boost/graph/selection.h>
50# include <CGAL/convex_hull_3.h>
51# include <CGAL/make_mesh_3.h>
52# include <CGAL/make_surface_mesh.h>
54//# include <CGAL/tetrahedral_remeshing.h> REQUIRES CGAL_VERSION>=5.1.5
55
56# include <fstream>
57# include <type_traits>
58
59
60
74namespace CGALWrappers
75{
76# ifdef CGAL_CONCURRENT_MESH_3
77 using ConcurrencyTag = CGAL::Parallel_tag;
78# else
79 using ConcurrencyTag = CGAL::Sequential_tag;
80# endif
81
92 {
100 compute_corefinement = 1 << 0,
101
108 compute_difference = 1 << 1,
109
115 compute_intersection = 1 << 2,
116
122 compute_union = 1 << 3,
123 };
124
125
126
137 template <typename C3t3>
138 void
140 CGAL::Surface_mesh<typename C3t3::Point::Point> &surface_mesh,
141 C3t3 & triangulation,
142 const AdditionalData<3> &data = AdditionalData<3>{});
143
187 template <typename CGALPointType>
188 void
190 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_1,
191 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_2,
192 const BooleanOperation & boolean_operation,
193 CGAL::Surface_mesh<CGALPointType> & output_surface_mesh);
194
206 template <typename CGALTriangulationType>
208 compute_quadrature(const CGALTriangulationType &tria,
209 const unsigned int degree);
210
226 template <int dim0, int dim1, int spacedim>
229 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
230 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
231 const unsigned int degree,
232 const BooleanOperation & bool_op,
233 const Mapping<dim0, spacedim> &mapping0 =
234 (::ReferenceCells::get_hypercube<dim0>()
235 .template get_default_linear_mapping<dim0, spacedim>()),
236 const Mapping<dim1, spacedim> &mapping1 =
237 (::ReferenceCells::get_hypercube<dim1>()
238 .template get_default_linear_mapping<dim1, spacedim>()));
239
255 template <int dim0, int dim1, int spacedim>
258 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
259 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
260 const unsigned int degree,
261 const Mapping<dim0, spacedim> &mapping0,
262 const Mapping<dim1, spacedim> &mapping1);
263
292 template <typename CGALPointType>
293 void
294 remesh_surface(CGAL::Surface_mesh<CGALPointType> &surface_mesh,
295 const AdditionalData<3> & data = AdditionalData<3>{});
296} // namespace CGALWrappers
297
298# ifndef DOXYGEN
299// Template implementations
300namespace CGALWrappers
301{
302 template <typename C3t3>
303 void
305 CGAL::Surface_mesh<typename C3t3::Point::Point> &surface_mesh,
306 C3t3 & triangulation,
307 const AdditionalData<3> & data)
308 {
309 using CGALPointType = typename C3t3::Point::Point;
310 Assert(CGAL::is_closed(surface_mesh),
311 ExcMessage("The surface mesh must be closed."));
312
313 using K = typename CGAL::Kernel_traits<CGALPointType>::Kernel;
314 using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<
315 K,
316 CGAL::Surface_mesh<CGALPointType>>;
317 using Tr = typename CGAL::
318 Mesh_triangulation_3<Mesh_domain, CGAL::Default, ConcurrencyTag>::type;
319 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
320
321 CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh);
322 Mesh_domain domain(surface_mesh);
323 domain.detect_features();
324 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
325 CGAL::parameters::facet_angle = data.facet_angle,
326 CGAL::parameters::facet_distance =
327 data.facet_distance,
328 CGAL::parameters::cell_radius_edge_ratio =
329 data.cell_radius_edge_ratio,
330 CGAL::parameters::cell_size = data.cell_size);
331 // Mesh generation
332 triangulation = CGAL::make_mesh_3<C3t3>(domain,
333 criteria,
334 CGAL::parameters::no_perturb(),
335 CGAL::parameters::no_exude());
336 }
337
338
339
340 template <typename CGALPointType>
341 void
343 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_1,
344 const CGAL::Surface_mesh<CGALPointType> &surface_mesh_2,
345 const BooleanOperation & boolean_operation,
346 CGAL::Surface_mesh<CGALPointType> & output_surface_mesh)
347 {
348 Assert(output_surface_mesh.is_empty(),
350 "output_surface_mesh must be empty upon calling this function"));
351 Assert(CGAL::is_closed(surface_mesh_1),
353 "The input surface_mesh_1 must be a closed surface mesh."));
354 Assert(CGAL::is_closed(surface_mesh_2),
356 "The input surface_mesh_2 must be a closed surface mesh."));
357 Assert(CGAL::is_triangle_mesh(surface_mesh_1),
358 ExcMessage("The first CGAL mesh must be triangulated."));
359 Assert(CGAL::is_triangle_mesh(surface_mesh_2),
360 ExcMessage("The second CGAL mesh must be triangulated."));
361
362 bool res = false;
363 auto surf_1 = surface_mesh_1;
364 auto surf_2 = surface_mesh_2;
365 namespace PMP = CGAL::Polygon_mesh_processing;
366 switch (boolean_operation)
367 {
369 res = PMP::corefine_and_compute_union(surf_1,
370 surf_2,
371 output_surface_mesh);
372 break;
374 res = PMP::corefine_and_compute_intersection(surf_1,
375 surf_2,
376 output_surface_mesh);
377 break;
379 res = PMP::corefine_and_compute_difference(surf_1,
380 surf_2,
381 output_surface_mesh);
382 break;
384 PMP::corefine(surf_1,
385 surf_2); // both surfaces are corefined
386 output_surface_mesh = std::move(surf_1);
387 res = true;
388 break;
389 default:
390 output_surface_mesh.clear();
391 break;
392 }
393 (void)res;
394 Assert(res,
395 ExcMessage("The boolean operation was not successfully computed."));
396 }
397
398
399
400 template <typename CGALTriangulationType>
402 compute_quadrature(const CGALTriangulationType &tria,
403 const unsigned int degree)
404 {
405 Assert(tria.is_valid(), ExcMessage("The triangulation is not valid."));
406 Assert(CGALTriangulationType::Point::Ambient_dimension::value == 3,
408 Assert(degree > 0,
409 ExcMessage("The degree of the Quadrature formula is not positive."));
410
411 constexpr int spacedim =
412 CGALTriangulationType::Point::Ambient_dimension::value;
413 QGaussSimplex<spacedim> quad(degree);
414 std::vector<::Point<spacedim>> pts;
415 std::vector<double> wts;
416 std::array<::Point<spacedim>, spacedim + 1> vertices; // tets
417
418 const auto is_c3t3 = boost::hana::is_valid(
419 [](auto &&obj) -> decltype(obj.cells_in_complex_begin()) {});
420
421 const auto is_tria3 = boost::hana::is_valid(
422 [](auto &&obj) -> decltype(obj.finite_cell_handles()) {});
423
424 if constexpr (decltype(is_c3t3(tria)){})
425 {
426 for (auto iter = tria.cells_in_complex_begin();
427 iter != tria.cells_in_complex_end();
428 ++iter)
429 {
430 for (unsigned int i = 0; i < (spacedim + 1); ++i)
431 {
432 vertices[i] = cgal_point_to_dealii_point<spacedim>(
433 iter->vertex(i)->point());
434 }
435
436 auto local_quad = quad.compute_affine_transformation(vertices);
437 std::transform(local_quad.get_points().begin(),
438 local_quad.get_points().end(),
439 std::back_inserter(pts),
440 [&pts](const auto &p) { return p; });
441 std::transform(local_quad.get_weights().begin(),
442 local_quad.get_weights().end(),
443 std::back_inserter(wts),
444 [&wts](const double w) { return w; });
445 }
446 }
447 else if constexpr (decltype(is_tria3(tria)){})
448 {
449 for (const auto &f : tria.finite_cell_handles())
450 {
451 for (unsigned int i = 0; i < (spacedim + 1); ++i)
452 {
453 vertices[i] =
454 cgal_point_to_dealii_point<spacedim>(f->vertex(i)->point());
455 }
456
457 auto local_quad = quad.compute_affine_transformation(vertices);
458 std::transform(local_quad.get_points().begin(),
459 local_quad.get_points().end(),
460 std::back_inserter(pts),
461 [&pts](const auto &p) { return p; });
462 std::transform(local_quad.get_weights().begin(),
463 local_quad.get_weights().end(),
464 std::back_inserter(wts),
465 [&wts](const double w) { return w; });
466 }
467 }
468 else
469 {
470 Assert(false, ExcMessage("Not a valid CGAL Triangulation."));
471 }
472 return Quadrature<spacedim>(pts, wts);
473 }
474
475
476
477 template <int dim0, int dim1, int spacedim>
480 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
481 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
482 const unsigned int degree,
483 const BooleanOperation & bool_op,
484 const Mapping<dim0, spacedim> &mapping0,
485 const Mapping<dim1, spacedim> &mapping1)
486 {
487 Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
488 ExcNotImplemented("2D geometries are not yet supported."));
489 if (dim1 > dim0)
490 {
492 cell1,
493 cell0,
494 degree,
495 bool_op,
496 mapping1,
497 mapping0); // This function works for dim1<=dim0, so swap them if this
498 // is not the case.
499 }
501 {
503 cell0, cell1, degree, mapping0, mapping1);
504 }
505 else
506 {
507 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
508 using CGALPoint = CGAL::Point_3<K>;
509 CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
510 dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
511 dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
512 // They have to be triangle meshes
513 CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
514 CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
515 Assert(CGAL::is_triangle_mesh(surface_1),
516 ExcMessage("The surface must be a triangle mesh."));
517 compute_boolean_operation(surface_1, surface_2, bool_op, out_surface);
518
519 using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<
520 K,
521 CGAL::Surface_mesh<CGALPoint>>;
522 using Tr = typename CGAL::Mesh_triangulation_3<Mesh_domain,
523 CGAL::Default,
524 ConcurrencyTag>::type;
525 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
526 using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
527 C3t3 tria;
529 return compute_quadrature(tria, degree);
530 }
531 }
532
533
534
535 template <int dim0, int dim1, int spacedim>
538 const typename ::Triangulation<dim0, spacedim>::cell_iterator &cell0,
539 const typename ::Triangulation<dim1, spacedim>::cell_iterator &cell1,
540 const unsigned int degree,
541 const Mapping<dim0, spacedim> &mapping0,
542 const Mapping<dim1, spacedim> &mapping1)
543 {
544 Assert(dim0 == 3 && dim1 == 3 && spacedim == 3,
545 ExcNotImplemented("2D geometries are not yet supported."));
546 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
547 using CGALPoint = CGAL::Point_3<K>;
548 using CGALTriangulation = CGAL::Triangulation_3<K>;
549
550 CGAL::Surface_mesh<CGALPoint> surface_1, surface_2, out_surface;
551 dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surface_1);
552 dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surface_2);
553 // They have to be triangle meshes
554 CGAL::Polygon_mesh_processing::triangulate_faces(surface_1);
555 CGAL::Polygon_mesh_processing::triangulate_faces(surface_2);
556
558 surface_2,
560 out_surface);
561 CGAL::Surface_mesh<CGALPoint> dummy;
562 CGALTriangulation tr;
563 CGAL::convex_hull_3(out_surface.points().begin(),
564 out_surface.points().end(),
565 dummy);
566 tr.insert(dummy.points().begin(), dummy.points().end());
567 return compute_quadrature(tr, degree);
568 }
569
570
571
572 template <typename CGALPointType>
573 void
574 remesh_surface(CGAL::Surface_mesh<CGALPointType> &cgal_mesh,
575 const AdditionalData<3> & data)
576 {
577 using K = CGAL::Exact_predicates_inexact_constructions_kernel;
578 using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
579 // Polyhedron type
580 using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
581 // Triangulation
582 using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
583 using C3t3 =
584 CGAL::Mesh_complex_3_in_triangulation_3<Tr,
585 Mesh_domain::Corner_index,
586 Mesh_domain::Curve_index>;
587 // Criteria
588 using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
589
590 Polyhedron poly;
591 CGAL::copy_face_graph(cgal_mesh, poly);
592 // Create a vector with only one element: the pointer to the
593 // polyhedron.
594 std::vector<Polyhedron *> poly_ptrs_vector(1, &poly);
595 // Create a polyhedral domain, with only one polyhedron,
596 // and no "bounding polyhedron", so the volumetric part of the
597 // domain will be empty.
598 Mesh_domain domain(poly_ptrs_vector.begin(), poly_ptrs_vector.end());
599 // Get sharp features
600 domain.detect_features(); // includes detection of borders
601
602 // Mesh criteria
603 const double default_value_edge_size = std::numeric_limits<double>::max();
604 if (data.edge_size > 0 &&
605 std::abs(data.edge_size - default_value_edge_size) > 1e-12)
606 {
607 Mesh_criteria criteria(CGAL::parameters::edge_size = data.edge_size,
608 CGAL::parameters::facet_size = data.facet_size,
609 CGAL::parameters::facet_angle = data.facet_angle,
610 CGAL::parameters::facet_distance =
611 data.facet_distance,
612 CGAL::parameters::cell_radius_edge_ratio =
613 data.cell_radius_edge_ratio,
614 CGAL::parameters::cell_size = data.cell_size);
615 // Mesh generation
616 C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
617 criteria,
618 CGAL::parameters::no_perturb(),
619 CGAL::parameters::no_exude());
620 cgal_mesh.clear();
621 CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
622 }
623 else if (std::abs(data.edge_size - default_value_edge_size) <= 1e-12)
624 {
625 Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
626 CGAL::parameters::facet_angle = data.facet_angle,
627 CGAL::parameters::facet_distance =
628 data.facet_distance,
629 CGAL::parameters::cell_radius_edge_ratio =
630 data.cell_radius_edge_ratio,
631 CGAL::parameters::cell_size = data.cell_size);
632 // Mesh generation
633 C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
634 criteria,
635 CGAL::parameters::no_perturb(),
636 CGAL::parameters::no_exude());
637 cgal_mesh.clear();
638 CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
639 }
640 else
641 {
642 Assert(false, ExcInternalError());
643 }
644 }
645} // namespace CGALWrappers
646# endif
647
649
650#endif
651#endif
Abstract base class for mapping classes.
Definition: mapping.h:311
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
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)
::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 >()))
CGAL::Sequential_tag ConcurrencyTag
Definition: utilities.h:79
void cgal_surface_mesh_to_cgal_triangulation(CGAL::Surface_mesh< typename C3t3::Point::Point > &surface_mesh, C3t3 &triangulation, const AdditionalData< 3 > &data=AdditionalData< 3 >{})
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)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const ::Triangulation< dim, spacedim > & tria