]> https://gitweb.dealii.org/ - dealii.git/commitdiff
changes to CGAL functions 18465/head
authorSascha Hofstetter <hofssasc@gmail.com>
Fri, 16 May 2025 08:17:46 +0000 (10:17 +0200)
committerSascha Hofstetter <hofssasc@gmail.com>
Wed, 21 May 2025 12:18:16 +0000 (14:18 +0200)
Fix the problem of not closed and oriented meshes

clean up

fix indent

changes in test files

clean up as proposed by bergbauer

indent

give correct output for test case

patch credit added

forgot indent

doc/news/changes/minor/20250520Hofstetter [new file with mode: 0644]
source/cgal/surface_mesh.cc
tests/cgal/cgal_remesh_surface.cc
tests/cgal/cgal_surface_mesh_01.cc
tests/cgal/cgal_surface_mesh_01.output
tests/cgal/cgal_surface_mesh_01.output.cgal-5.4.1
tests/cgal/cgal_surface_mesh_04.cc
tests/cgal/cgal_surface_mesh_05.cc

diff --git a/doc/news/changes/minor/20250520Hofstetter b/doc/news/changes/minor/20250520Hofstetter
new file mode 100644 (file)
index 0000000..8673450
--- /dev/null
@@ -0,0 +1,6 @@
+Fixed: CGALWrappers dealii_cell_to_cgal_surface_mesh
+and dealii_tria_to_cgal_surface_mesh now consistently
+return closed and oriented surface meshes for tets
+and hex. This avoids errors in boolean operations.
+<br>
+(Sascha Hofstetter, 2025/05/20)
index aaee70f2c0460456b880dd82aafb4215f184c7c2..973643028d58f790d2fd5c4f79019d2194ba15c0 100644 (file)
@@ -42,27 +42,40 @@ namespace
       {
         mesh.add_edge(deal2cgal.at(face->vertex_index(0)),
                       deal2cgal.at(face->vertex_index(1)));
+
+        if (clockwise_ordering)
+          std::reverse(indices.begin(), indices.end());
       }
     else if (reference_cell_type == ReferenceCells::Triangle)
       {
         indices = {deal2cgal.at(face->vertex_index(0)),
                    deal2cgal.at(face->vertex_index(1)),
                    deal2cgal.at(face->vertex_index(2))};
+
+        if (clockwise_ordering)
+          std::reverse(indices.begin(), indices.end());
       }
 
     else if (reference_cell_type == ReferenceCells::Quadrilateral)
       {
-        indices = {deal2cgal.at(face->vertex_index(0)),
-                   deal2cgal.at(face->vertex_index(1)),
-                   deal2cgal.at(face->vertex_index(3)),
-                   deal2cgal.at(face->vertex_index(2))};
+        if (clockwise_ordering)
+          {
+            indices = {deal2cgal.at(face->vertex_index(0)),
+                       deal2cgal.at(face->vertex_index(2)),
+                       deal2cgal.at(face->vertex_index(3)),
+                       deal2cgal.at(face->vertex_index(1))};
+          }
+        else
+          {
+            indices = {deal2cgal.at(face->vertex_index(0)),
+                       deal2cgal.at(face->vertex_index(1)),
+                       deal2cgal.at(face->vertex_index(3)),
+                       deal2cgal.at(face->vertex_index(2))};
+          }
       }
     else
       DEAL_II_ASSERT_UNREACHABLE();
 
-    if (clockwise_ordering == true)
-      std::reverse(indices.begin(), indices.end());
-
     [[maybe_unused]] const auto new_face = mesh.add_face(indices);
     Assert(new_face != mesh.null_face(),
            ExcInternalError("While trying to build a CGAL facet, "
@@ -81,9 +94,12 @@ namespace
   {
     for (const auto i : cell->vertex_indices())
       {
-        deal2cgal[cell->vertex_index(i)] = mesh.add_vertex(
-          CGALWrappers::dealii_point_to_cgal_point<typename CGAL_Mesh::Point>(
-            cell->vertex(i)));
+        if (deal2cgal.find(cell->vertex_index(i)) == deal2cgal.end())
+          {
+            deal2cgal[cell->vertex_index(i)] =
+              mesh.add_vertex(CGALWrappers::dealii_point_to_cgal_point<
+                              typename CGAL_Mesh::Point>(cell->vertex(i)));
+          }
       }
   }
 } // namespace
@@ -178,15 +194,16 @@ namespace CGALWrappers
         for (const auto &cell : tria.active_cell_iterators())
           {
             for (const auto &f : cell->face_indices())
-
-              if (cell->face(f)->at_boundary())
-                {
-                  map_vertices(cell->face(f), deal2cgal, mesh);
-                  add_facet(cell->face(f),
-                            deal2cgal,
-                            mesh,
-                            (f % 2 == 0 || cell->n_vertices() != 8));
-                }
+              {
+                if (cell->face(f)->at_boundary())
+                  {
+                    map_vertices(cell->face(f), deal2cgal, mesh);
+                    add_facet(cell->face(f),
+                              deal2cgal,
+                              mesh,
+                              (f % 2 == 0 || cell->n_vertices() != 8));
+                  }
+              }
           }
       }
     else
index cc014f52e4d27d70fbf0d80014429258c2503004..e9878333901d34096b7247469536b6a133e4e1ed 100644 (file)
@@ -27,6 +27,7 @@
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>
 
+#include <CGAL/Polygon_mesh_processing/orientation.h>
 #include <string.h>
 
 #include "../tests.h"
@@ -66,9 +67,18 @@ test()
   dealii_tria_to_cgal_surface_mesh(tria0, surface_mesh0);
   dealii_tria_to_cgal_surface_mesh(tria1, surface_mesh1);
 
-  // close the surfaces
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh0);
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh1);
+  Assert(surface_mesh0.is_valid() && surface_mesh1.is_valid(),
+         ExcMessage("The CGAL surface mesh is not valid."));
+  if (dim == 3)
+    {
+      Assert(CGAL::is_closed(surface_mesh0) && CGAL::is_closed(surface_mesh1),
+             dealii::ExcMessage("The CGAL mesh is not closed"));
+      Assert(
+        CGAL::Polygon_mesh_processing::is_outward_oriented(surface_mesh0) &&
+          CGAL::Polygon_mesh_processing::is_outward_oriented(surface_mesh1),
+        dealii::ExcMessage(
+          "The normal vectors of the CGAL mesh are not oriented outwards"));
+    }
 
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0);
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1);
index 8e7624b6ec26cdbb461f878aa52f7c66a64b64e1..f05e16a0cce3f98bcb22f7db49f1587e7630c48d 100644 (file)
@@ -26,6 +26,7 @@
 #include <deal.II/grid/tria.h>
 
 #include <CGAL/IO/io.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
 
 #include "../tests.h"
 
@@ -56,7 +57,24 @@ test()
       const auto cell = tria.begin_active();
       dealii_cell_to_cgal_surface_mesh(cell, *mapping, mesh);
 
-      Assert(mesh.is_valid(), dealii::ExcMessage("The CGAL mesh is not valid"));
+      Assert(mesh.is_valid(), ExcMessage("The CGAL mesh is not valid"));
+
+      if (dim == 3)
+        {
+          // is closed/oriented only for 3 dimensional objects important
+          Assert(CGAL::is_closed(mesh),
+                 ExcMessage("The CGAL mesh is not closed"));
+
+          // orientation not supported for wedges, this needs special treatment
+          if (r_cell != ref_cells[3][2])
+            {
+              Assert(
+                CGAL::Polygon_mesh_processing::is_outward_oriented(mesh),
+                ExcMessage(
+                  "The normal vectors of the CGAL mesh are not oriented outwards"));
+            }
+        }
+
       deallog << "deal vertices: " << cell->n_vertices() << ", cgal vertices "
               << mesh.num_vertices() << std::endl;
       deallog << "deal faces: " << cell->n_faces() << ", cgal faces "
index 6a1b26a53fee7bd802828162a26b1b207c3669f5..d92dbd368c0f82a0ec7b027eba12194fc45d3980 100644 (file)
@@ -19,7 +19,7 @@ DEAL::OFF
 1.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-4 2 3 1 0
+4 0 2 3 1
 
 DEAL::dim= 2,   spacedim= 3
 DEAL::deal vertices: 3, cgal vertices 3
@@ -41,7 +41,7 @@ DEAL::OFF
 1.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-4 2 3 1 0
+4 0 2 3 1
 
 DEAL::dim= 3,   spacedim= 3
 DEAL::deal vertices: 4, cgal vertices 4
@@ -68,7 +68,7 @@ DEAL::OFF
 -1.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
 0.00000 0.00000 1.00000
-4 2 3 1 0
+4 0 2 3 1
 3 4 2 0
 3 4 1 3
 3 4 0 1
@@ -87,9 +87,9 @@ DEAL::OFF
 0.00000 1.00000 1.00000
 3 2 0 1
 3 5 4 3
-4 3 4 1 0
-4 4 5 2 1
-4 5 3 0 2
+4 0 3 4 1
+4 1 4 5 2
+4 2 5 3 0
 
 DEAL::deal vertices: 8, cgal vertices 8
 DEAL::deal faces: 6, cgal faces 6
@@ -104,10 +104,10 @@ DEAL::OFF
 1.00000 0.00000 1.00000
 0.00000 1.00000 1.00000
 1.00000 1.00000 1.00000
-4 4 6 2 0
+4 0 4 6 2
 4 1 3 7 5
-4 1 5 4 0
+4 0 1 5 4
 4 2 6 7 3
-4 2 3 1 0
+4 0 2 3 1
 4 4 5 7 6
 
index 2bf9a298009e9b583e6736d5577cb891e20ea3c6..a77fd3f13f9bc18aea419d086a533283ad81d47f 100644 (file)
@@ -22,7 +22,7 @@ DEAL::OFF
 1.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-4  2 3 1 0
+4  0 2 3 1
 
 
 DEAL::dim= 2,   spacedim= 3
@@ -48,7 +48,7 @@ DEAL::OFF
 1.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-4  2 3 1 0
+4  0 2 3 1
 
 
 DEAL::dim= 3,   spacedim= 3
@@ -79,7 +79,7 @@ DEAL::OFF
 -1.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
 0.00000 0.00000 1.00000
-4  2 3 1 0
+4  0 2 3 1
 3  4 2 0
 3  4 1 3
 3  4 0 1
@@ -100,9 +100,9 @@ DEAL::OFF
 0.00000 1.00000 1.00000
 3  2 0 1
 3  5 4 3
-4  3 4 1 0
-4  4 5 2 1
-4  5 3 0 2
+4  0 3 4 1
+4  1 4 5 2
+4  2 5 3 0
 
 
 DEAL::deal vertices: 8, cgal vertices 8
@@ -119,11 +119,11 @@ DEAL::OFF
 1.00000 0.00000 1.00000
 0.00000 1.00000 1.00000
 1.00000 1.00000 1.00000
-4  4 6 2 0
+4  0 4 6 2
 4  1 3 7 5
-4  1 5 4 0
+4  0 1 5 4
 4  2 6 7 3
-4  2 3 1 0
+4  0 2 3 1
 4  4 5 7 6
 
 
index 98477f1083ee5353e1696f458280a037a1bbd823..305ff71cc9ae28ed46c0a67d18b167ce76acd29e 100644 (file)
@@ -27,6 +27,7 @@
 #include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/tria.h>
 
+#include <CGAL/Polygon_mesh_processing/orientation.h>
 #include <string.h>
 
 #include "../tests.h"
@@ -76,6 +77,16 @@ test()
       Assert(surface_mesh.is_valid(),
              ExcMessage("The CGAL surface mesh is not valid."));
 
+      if (dim == 3)
+        {
+          Assert(CGAL::is_closed(surface_mesh),
+                 ExcMessage("The CGAL mesh is not closed"));
+          Assert(
+            CGAL::Polygon_mesh_processing::is_outward_oriented(surface_mesh),
+            ExcMessage(
+              "The normal vectors of the CGAL mesh are not oriented outwards"));
+        }
+
       // Now back to the original dealii tria.
       cgal_surface_mesh_to_dealii_triangulation(surface_mesh, tria_out);
       std::ofstream out_name(name + std::to_string(spacedim) + ".vtk");
index 397edfe09a3aae9c1e38f5a7b84f4d5c9f81fb6a..247d5d6ac19124c09ece57e1faacc4c3b9d91526 100644 (file)
@@ -57,10 +57,14 @@ test()
   dealii_tria_to_cgal_surface_mesh(tria0, surface_mesh0);
   dealii_tria_to_cgal_surface_mesh(tria1, surface_mesh1);
 
-  // close the surfaces
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh0);
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh1);
-
+  // Ensure the meshes are closed
+  Assert(CGAL::is_closed(surface_mesh0),
+         ExcMessage("The CGAL mesh 0 is not closed"));
+  Assert(CGAL::is_closed(surface_mesh1),
+         ExcMessage("The CGAL mesh 1 is not closed"));
+
+  // Surfaces automatically closed but still need to be triangulated
+  // before using compute_boolean_operation
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0);
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1);
 
@@ -93,10 +97,14 @@ test()
   dealii_tria_to_cgal_surface_mesh(tria0, surface_mesh0);
   dealii_tria_to_cgal_surface_mesh(tria1, surface_mesh1);
 
-  // close the surfaces
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh0);
-  CGAL::Polygon_mesh_processing::stitch_borders(surface_mesh1);
+  // Ensure the meshes are closed
+  Assert(CGAL::is_closed(surface_mesh0),
+         ExcMessage("The CGAL mesh 0 is not closed"));
+  Assert(CGAL::is_closed(surface_mesh1),
+         ExcMessage("The CGAL mesh 1 is not closed"));
 
+  // Surfaces automatically closed but still need to be triangulated
+  // before using compute_boolean_operation
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0);
   CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1);
 

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.