From 14011ba524705735ec8e572d17a12390b5039fbb Mon Sep 17 00:00:00 2001 From: Sascha Hofstetter Date: Fri, 16 May 2025 10:17:46 +0200 Subject: [PATCH] changes to CGAL functions 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 | 6 ++ source/cgal/surface_mesh.cc | 55 ++++++++++++------- tests/cgal/cgal_remesh_surface.cc | 16 +++++- tests/cgal/cgal_surface_mesh_01.cc | 20 ++++++- tests/cgal/cgal_surface_mesh_01.output | 18 +++--- .../cgal_surface_mesh_01.output.cgal-5.4.1 | 18 +++--- tests/cgal/cgal_surface_mesh_04.cc | 11 ++++ tests/cgal/cgal_surface_mesh_05.cc | 22 +++++--- 8 files changed, 118 insertions(+), 48 deletions(-) create mode 100644 doc/news/changes/minor/20250520Hofstetter diff --git a/doc/news/changes/minor/20250520Hofstetter b/doc/news/changes/minor/20250520Hofstetter new file mode 100644 index 0000000000..8673450439 --- /dev/null +++ b/doc/news/changes/minor/20250520Hofstetter @@ -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. +
+(Sascha Hofstetter, 2025/05/20) diff --git a/source/cgal/surface_mesh.cc b/source/cgal/surface_mesh.cc index aaee70f2c0..973643028d 100644 --- a/source/cgal/surface_mesh.cc +++ b/source/cgal/surface_mesh.cc @@ -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( - 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 diff --git a/tests/cgal/cgal_remesh_surface.cc b/tests/cgal/cgal_remesh_surface.cc index cc014f52e4..e987833390 100644 --- a/tests/cgal/cgal_remesh_surface.cc +++ b/tests/cgal/cgal_remesh_surface.cc @@ -27,6 +27,7 @@ #include #include +#include #include #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); diff --git a/tests/cgal/cgal_surface_mesh_01.cc b/tests/cgal/cgal_surface_mesh_01.cc index 8e7624b6ec..f05e16a0cc 100644 --- a/tests/cgal/cgal_surface_mesh_01.cc +++ b/tests/cgal/cgal_surface_mesh_01.cc @@ -26,6 +26,7 @@ #include #include +#include #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 " diff --git a/tests/cgal/cgal_surface_mesh_01.output b/tests/cgal/cgal_surface_mesh_01.output index 6a1b26a53f..d92dbd368c 100644 --- a/tests/cgal/cgal_surface_mesh_01.output +++ b/tests/cgal/cgal_surface_mesh_01.output @@ -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 diff --git a/tests/cgal/cgal_surface_mesh_01.output.cgal-5.4.1 b/tests/cgal/cgal_surface_mesh_01.output.cgal-5.4.1 index 2bf9a29800..a77fd3f13f 100644 --- a/tests/cgal/cgal_surface_mesh_01.output.cgal-5.4.1 +++ b/tests/cgal/cgal_surface_mesh_01.output.cgal-5.4.1 @@ -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 diff --git a/tests/cgal/cgal_surface_mesh_04.cc b/tests/cgal/cgal_surface_mesh_04.cc index 98477f1083..305ff71cc9 100644 --- a/tests/cgal/cgal_surface_mesh_04.cc +++ b/tests/cgal/cgal_surface_mesh_04.cc @@ -27,6 +27,7 @@ #include #include +#include #include #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"); diff --git a/tests/cgal/cgal_surface_mesh_05.cc b/tests/cgal/cgal_surface_mesh_05.cc index 397edfe09a..247d5d6ac1 100644 --- a/tests/cgal/cgal_surface_mesh_05.cc +++ b/tests/cgal/cgal_surface_mesh_05.cc @@ -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); -- 2.39.5