--- /dev/null
+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)
{
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, "
{
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
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
#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"
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);
#include <deal.II/grid/tria.h>
#include <CGAL/IO/io.h>
+#include <CGAL/Polygon_mesh_processing/orientation.h>
#include "../tests.h"
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 "
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
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
-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
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
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
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
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
-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
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
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
#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"
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");
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);
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);