]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add 3D/3D cuts with vertices only and minor fixes 14166/head
authorMarco Feder <marco.feder@sissa.it>
Fri, 9 Sep 2022 07:57:01 +0000 (09:57 +0200)
committerMarco Feder <marco.feder@sissa.it>
Tue, 20 Sep 2022 15:32:14 +0000 (17:32 +0200)
cmake/modules/FindCGAL.cmake
doc/news/changes/minor/20220920Feder [new file with mode: 0644]
include/deal.II/base/quadrature_lib.h
include/deal.II/cgal/intersections.h
include/deal.II/cgal/utilities.h
source/cgal/intersections.cc

index 34692bee59b627d099078679c39723e85ab94cf8..7071d5b31ea9d339dd71113892b65258183d3957 100644 (file)
@@ -1,17 +1,17 @@
-# # ---------------------------------------------------------------------
-# #
-# # Copyright (C) 2022 by the deal.II authors
-# #
-# # This file is part of the deal.II library.
-# #
-# # The deal.II library is free software; you can use it, redistribute
-# # it, and/or modify it under the terms of the GNU Lesser General
-# # Public License as published by the Free Software Foundation; either
-# # version 2.1 of the License, or (at your option) any later version.
-# # The full text of the license can be found in the file LICENSE.md at
-# # the top level directory of deal.II.
-# #
-# # ---------------------------------------------------------------------
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2022 by the deal.II authors
+##
+## This file is part of the deal.II library.
+##
+## The deal.II library is free software; you can use it, redistribute
+## it, and/or modify it under the terms of the GNU Lesser General
+## Public License as published by the Free Software Foundation; either
+## version 2.1 of the License, or (at your option) any later version.
+## The full text of the license can be found in the file LICENSE.md at
+## the top level directory of deal.II.
+##
+## ---------------------------------------------------------------------
 
 #
 # Try to find the CGAL libraries
diff --git a/doc/news/changes/minor/20220920Feder b/doc/news/changes/minor/20220920Feder
new file mode 100644 (file)
index 0000000..aaa411a
--- /dev/null
@@ -0,0 +1,6 @@
+New: The new function CGALWrappers::compute_intersection_of_cells computes
+the intersection between two (affine) cells starting from the location of the
+vertices. The intersection is described by a vector of arrays, where each array 
+identifies a simplex of the sub-tassellation of the intersection.
+<br>
+(Marco Feder, Johannes Heinz, 2022/09/20)
index 5284bf849c320ba26fac4f195742c44925c6de27..0cfe1fa4ee26fe19cff85b7da471a7f91e00fe4d 100644 (file)
@@ -644,12 +644,12 @@ public:
     const std::array<Point<spacedim>, dim + 1> &vertices) const;
 
   /**
-   *
-   * Given a partition of a cell into simplices, this function creates a
-   * quadrature rule on the cell by collecting Quadrature objects on each
-   * simplex. A simplex is identified by its vertices, which are stored into an
-   * array of Points. Hence, this function can provide quadrature rules on
-   * polygons (or polyhedra), as they can be split into simplices.
+   * Given a collection of simplices, this function creates a global
+   * quadrature rule on the area covered by the simplices, by mapping the
+   * current quadrature on each simplex. A simplex is identified by its
+   * vertices, which are stored into an array of Points. Hence, this function
+   * can provide quadrature rules on polygons (or polyhedra), as they can be
+   * split into simplices.
    *
    *
    * @param simplices A std::vector where each entry is an array of `dim+1` points, which identifies the vertices of a simplex.
index 8254d10e35c74629830b033fed2c60fc3dd21c53..acf05fae360a2d35902b07df5e5e25b307b46a02 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 namespace CGALWrappers
 {
-  template <int dim0,
-            int dim1,
-            int spacedim,
-            int n_components0,
-            int n_components1>
-  std::vector<std::array<Point<spacedim>, dim1 + 1>>
-  compute_intersection_of_cells(
-    const std::array<Point<spacedim>, n_components0> &vertices0,
-    const std::array<Point<spacedim>, n_components1> &vertices1,
-    const double                                      tol)
-  {
-    (void)vertices0;
-    (void)vertices1;
-    (void)tol;
-    Assert(false, ExcMessage("No explicit template instantiation available"));
-    return {};
-  }
-
   /**
-   * Given two deal.II affine cells, compute the intersection and return a
+   * Given two deal.II affine cells, compute their intersection and return a
    * vector of simplices, each one identified by an array of deal.II Points. All
    * the simplices together are a subdivision of the intersection. If cells are
    * non-affine, a geometrical error is introduced. If the
@@ -75,6 +57,24 @@ namespace CGALWrappers
     const Mapping<dim0, spacedim> &                              mapping0,
     const Mapping<dim1, spacedim> &                              mapping1,
     const double                                                 tol = 1e-9);
+
+
+  /**
+   * Same function as above, but working directly with vertices.
+   */
+  template <int dim0, int dim1, int spacedim, int n_vertices0, int n_vertices1>
+  std::vector<std::array<Point<spacedim>, dim1 + 1>>
+  compute_intersection_of_cells(
+    const std::array<Point<spacedim>, n_vertices0> &vertices0,
+    const std::array<Point<spacedim>, n_vertices1> &vertices1,
+    const double                                    tol = 1e-9)
+  {
+    (void)vertices0;
+    (void)vertices1;
+    (void)tol;
+    Assert(false, ExcMessage("No explicit template instantiation available"));
+    return {};
+  }
 } // namespace CGALWrappers
 
 DEAL_II_NAMESPACE_CLOSE
index 4e8bc23a80ab1b382a90e032eba6e76bed39f436..10fe36523fd8c2b3cf7eac85513ae1486c473f31 100644 (file)
@@ -595,11 +595,11 @@ namespace CGALWrappers
 
 
   /**
-   * Get vertices of cell in CGAL ordering
+   * Get vertices of cell in CGAL ordering.
    *
-   * @param [in] cell A cell_iterator to a deal.II cell.
-   * @param [in] mapping Mapping object for the cell.
-   * @return [out] Array of vertices in CGAL order.
+   * @param cell A cell_iterator to a deal.II cell.
+   * @param mapping Mapping object for the cell.
+   * @return  Array of vertices in CGAL order.
    */
   template <int n_vertices, int dim, int spacedim>
   std::array<Point<spacedim>, n_vertices>
index 7a859849b92fb583545cb907eb6370a1b4d7f53f..fbf5dd76b17c9a4240fca225920107182966f0c3 100644 (file)
@@ -70,6 +70,7 @@ namespace CGALWrappers
   using CGALPoint2             = K::Point_2;
   using CGALPoint3             = K::Point_3;
   using CGALPoint3_exact       = K_exact::Point_3;
+  using CGALPoint3_inexact     = K_inexact::Point_3;
   using CGALSegment2           = K::Segment_2;
   using Surface_mesh           = CGAL::Surface_mesh<K_inexact::Point_3>;
   using CGALSegment3           = K::Segment_3;
@@ -141,6 +142,8 @@ namespace CGALWrappers
         }
     }
 
+
+
     void
     mark_domains(CDT &cdt)
     {
@@ -200,6 +203,8 @@ namespace CGALWrappers
       return CGAL::intersection(triangle1, triangle2);
     }
 
+
+
     boost::optional<boost::variant<CGALPoint2, CGALSegment2>>
     compute_intersection(const std::array<Point<2>, 3> &first_simplex,
                          const std::array<Point<2>, 2> &second_simplex)
@@ -227,6 +232,8 @@ namespace CGALWrappers
       return CGAL::intersection(segm, triangle);
     }
 
+
+
     // rectangle-rectangle
     std::vector<Polygon_with_holes_2>
     compute_intersection(const std::array<Point<2>, 4> &first_simplex,
@@ -255,6 +262,8 @@ namespace CGALWrappers
       return poly_list;
     }
 
+
+
     boost::optional<boost::variant<CGALPoint3, CGALSegment3>>
     compute_intersection(const std::array<Point<3>, 2> &first_simplex,
                          const std::array<Point<3>, 4> &second_simplex)
@@ -335,6 +344,8 @@ namespace CGALWrappers
     }
   } // namespace internal
 
+
+
   // Specialization for quads
   template <>
   std::vector<std::array<Point<2>, 3>>
@@ -400,6 +411,8 @@ namespace CGALWrappers
       }
   }
 
+
+
   // Specialization for quad \cap line
   template <>
   std::vector<std::array<Point<2>, 2>>
@@ -601,51 +614,93 @@ namespace CGALWrappers
 #  endif
   }
 
+
+
   template <>
   std::vector<std::array<Point<3>, 4>>
   compute_intersection_of_cells<3, 3, 3, 8, 8>(
     const std::array<Point<3>, 8> &vertices0,
     const std::array<Point<3>, 8> &vertices1,
-    const double                   tol)
+
+    const double tol)
   {
-    // Surface_mesh surf0, surf1, sm;
-    // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell0, mapping0, surf0);
-    // CGALWrappers::dealii_cell_to_cgal_surface_mesh(cell1, mapping1, surf1);
-    // CGAL::Polygon_mesh_processing::triangulate_faces(surf0);
-    // CGAL::Polygon_mesh_processing::triangulate_faces(surf1);
-    // CGALWrappers::compute_boolean_operation(
-    //   surf0, surf1, CGALWrappers::BooleanOperation::compute_intersection,
-    //   sm);
-    // std::vector<std::array<Point<3>, 4>> vertices;
-    // if (CGAL::Polygon_mesh_processing::volume(sm) > tol)
-    //   {
-    //     // Collect tetrahedrons
-    //     Triangulation3_inexact tria;
-    //     tria.insert(sm.points().begin(), sm.points().end());
-    //     for (const auto &c : tria.finite_cell_handles())
-    //       {
-    //         const auto &tet = tria.tetrahedron(c);
-    //         vertices.push_back(
-    //           {{CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(0)),
-    //             CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(1)),
-    //             CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(2)),
-    //             CGALWrappers::cgal_point_to_dealii_point<3>(tet.vertex(3))}});
-    //       }
-    //     return vertices;
-    //   }
-    // else
-    //   {
-    //     return vertices;
-    //   }
-
-    // TODO: implememnt 3d/3d cut
-    AssertThrow(false, ExcMessage("Not yet implemented"));
-    (void)vertices0;
-    (void)vertices1;
-    (void)tol;
-    return {};
+    std::array<CGALPoint3_inexact, 8> pts_hex0;
+    std::array<CGALPoint3_inexact, 8> pts_hex1;
+    std::transform(
+      vertices0.begin(),
+      vertices0.end(),
+      pts_hex0.begin(),
+      [&](const Point<3> &p) {
+        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
+      });
+
+    std::transform(
+      vertices1.begin(),
+      vertices1.end(),
+      pts_hex1.begin(),
+      [&](const Point<3> &p) {
+        return CGALWrappers::dealii_point_to_cgal_point<CGALPoint3_inexact>(p);
+      });
+
+
+    Surface_mesh surf0, surf1, sm;
+    // Subdivide hex into tetrahedrons
+    std::vector<std::array<Point<3>, 4>> vertices;
+    Triangulation3_inexact               tria0, tria1;
+
+    tria0.insert(pts_hex0.begin(), pts_hex0.end());
+    tria1.insert(pts_hex1.begin(), pts_hex1.end());
+
+    for (const auto &c0 : tria0.finite_cell_handles())
+      {
+        const auto &                 tet0 = tria1.tetrahedron(c0);
+        [[maybe_unused]] const auto &tetg0 =
+          CGAL::make_tetrahedron(tet0.vertex(0),
+                                 tet0.vertex(1),
+                                 tet0.vertex(2),
+                                 tet0.vertex(3),
+                                 surf0);
+        for (const auto &c1 : tria1.finite_cell_handles())
+          {
+            const auto &                 tet1 = tria1.tetrahedron(c1);
+            [[maybe_unused]] const auto &tetg1 =
+              CGAL::make_tetrahedron(tet1.vertex(0),
+                                     tet1.vertex(1),
+                                     tet1.vertex(2),
+                                     tet1.vertex(3),
+                                     surf1);
+            const bool test_intersection =
+              CGAL::PMP::corefine_and_compute_intersection(surf0, surf1, sm);
+            if (CGAL::Polygon_mesh_processing::volume(sm) > tol &&
+                test_intersection)
+              {
+                // Collect tetrahedrons
+                Triangulation3_inexact tria;
+                tria.insert(sm.points().begin(), sm.points().end());
+                for (const auto &c : tria.finite_cell_handles())
+                  {
+                    const auto &tet = tria.tetrahedron(c);
+                    vertices.push_back(
+                      {{CGALWrappers::cgal_point_to_dealii_point<3>(
+                          tet.vertex(0)),
+                        CGALWrappers::cgal_point_to_dealii_point<3>(
+                          tet.vertex(1)),
+                        CGALWrappers::cgal_point_to_dealii_point<3>(
+                          tet.vertex(2)),
+                        CGALWrappers::cgal_point_to_dealii_point<3>(
+                          tet.vertex(3))}});
+                  }
+              }
+            surf1.clear();
+            sm.clear();
+          }
+        surf0.clear();
+      }
+    return vertices;
   }
 
+
+
   template <int dim0, int dim1, int spacedim>
   std::vector<std::array<Point<spacedim>, dim1 + 1>>
   compute_intersection_of_cells(
@@ -663,7 +718,6 @@ namespace CGALWrappers
     const auto vertices0 =
       CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim0))>(
         cell0, mapping0);
-
     const auto vertices1 =
       CGALWrappers::get_vertices_in_cgal_order<int(std::pow(2, dim1))>(
         cell1, mapping1);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.