# include <CGAL/Mesh_triangulation_3.h>
# include <CGAL/Polygon_mesh_processing/corefinement.h>
# include <CGAL/Polygon_mesh_processing/measure.h>
+# include <CGAL/Polygon_mesh_processing/remesh.h>
# include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
# include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
# include <CGAL/Simple_cartesian.h>
# include <CGAL/Surface_mesh.h>
# include <CGAL/Triangulation_3.h>
+# include <CGAL/boost/graph/copy_face_graph.h>
+# include <CGAL/boost/graph/selection.h>
# include <CGAL/convex_hull_3.h>
# include <CGAL/make_mesh_3.h>
# include <CGAL/make_surface_mesh.h>
# include <deal.II/cgal/surface_mesh.h>
+//# include <CGAL/tetrahedral_remeshing.h> REQUIRES CGAL_VERSION>=5.1.5
+# include <fstream>
# include <type_traits>
const unsigned int degree,
const Mapping<dim0, spacedim> &mapping0,
const Mapping<dim1, spacedim> &mapping1);
+
+ /**
+ * Remesh a CGAL::Surface_mesh.
+ *
+ * If the domain has 1-dimensional exposed features, the criteria includes a
+ * sizing field to guide the sampling of 1-dimensional features with
+ * protecting balls centers.
+ * - CGAL::parameters::edge_size.
+ *
+ * This utility should be exploited to improve the quality of a mesh coming
+ * from boolean operations. As an example, consider two CGAL::Surface_mesh
+ * objects describing two hyper_spheres that intersect non-trivially. After
+ * computing their corefinement and union with compute_boolean_operation(),
+ * the resulting mesh is the following:
+ *
+ * @image html boolean_union_hyper_spheres.png
+ *
+ * Clearly, elements (triangles) like the ones arising along the intersection
+ * of the two spheres should be avoided as they're quite degenerate and would
+ * decrease the accuracy of the results. A call to the present function will
+ * try to remesh the surface, according to the given named parameters, to
+ * improve its quality. In the present case, the result is the following:
+ *
+ * @image html boolean_union_hyper_spheres_remeshed.png
+ *
+ * @param surface_mesh The input CGAL::Surface_mesh.
+ * @param [in] data AdditionalData object to pass to the CGAL::make_mesh_3 function. See the documentation
+ * of that struct for a description of those parameters.
+ */
+ template <typename CGALPointType>
+ void
+ remesh_surface(CGAL::Surface_mesh<CGALPointType> &surface_mesh,
+ const AdditionalData<3> & data = AdditionalData<3>{});
} // namespace CGALWrappers
# ifndef DOXYGEN
tr.insert(dummy.points().begin(), dummy.points().end());
return compute_quadrature(tr, degree);
}
+
+
+
+ template <typename CGALPointType>
+ void
+ remesh_surface(CGAL::Surface_mesh<CGALPointType> &cgal_mesh,
+ const AdditionalData<3> & data)
+ {
+ using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+ using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
+ // Polyhedron type
+ using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
+ // Triangulation
+ using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
+ using C3t3 =
+ CGAL::Mesh_complex_3_in_triangulation_3<Tr,
+ Mesh_domain::Corner_index,
+ Mesh_domain::Curve_index>;
+ // Criteria
+ using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
+
+ Polyhedron poly;
+ CGAL::copy_face_graph(cgal_mesh, poly);
+ // Create a vector with only one element: the pointer to the
+ // polyhedron.
+ std::vector<Polyhedron *> poly_ptrs_vector(1, &poly);
+ // Create a polyhedral domain, with only one polyhedron,
+ // and no "bounding polyhedron", so the volumetric part of the
+ // domain will be empty.
+ Mesh_domain domain(poly_ptrs_vector.begin(), poly_ptrs_vector.end());
+ // Get sharp features
+ domain.detect_features(); // includes detection of borders
+
+ // Mesh criteria
+ const double default_value_edge_size = std::numeric_limits<double>::max();
+ if (data.edge_size > 0 &&
+ std::abs(data.edge_size - default_value_edge_size) > 1e-12)
+ {
+ Mesh_criteria criteria(CGAL::parameters::edge_size = data.edge_size,
+ CGAL::parameters::facet_size = data.facet_size,
+ CGAL::parameters::facet_angle = data.facet_angle,
+ CGAL::parameters::facet_distance =
+ data.facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio =
+ data.cell_radius_edge_ratio,
+ CGAL::parameters::cell_size = data.cell_size);
+ // Mesh generation
+ C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
+ criteria,
+ CGAL::parameters::no_perturb(),
+ CGAL::parameters::no_exude());
+ cgal_mesh.clear();
+ CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
+ }
+ else if (std::abs(data.edge_size - default_value_edge_size) <= 1e-12)
+ {
+ Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
+ CGAL::parameters::facet_angle = data.facet_angle,
+ CGAL::parameters::facet_distance =
+ data.facet_distance,
+ CGAL::parameters::cell_radius_edge_ratio =
+ data.cell_radius_edge_ratio,
+ CGAL::parameters::cell_size = data.cell_size);
+ // Mesh generation
+ C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
+ criteria,
+ CGAL::parameters::no_perturb(),
+ CGAL::parameters::no_exude());
+ cgal_mesh.clear();
+ CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, cgal_mesh);
+ }
+ else
+ {
+ Assert(false, ExcInternalError());
+ }
+ }
} // namespace CGALWrappers
# endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Remesh the union of two deal.II hyper_spheres in order to avoid bad
+// triangles.
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/cgal/surface_mesh.h>
+#include <deal.II/cgal/triangulation.h>
+#include <string.h>
+
+#include "../tests.h"
+
+using namespace CGALWrappers;
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using CGALPoint = CGAL::Point_3<K>;
+using CGALMesh = CGAL::Surface_mesh<CGALPoint>;
+
+template <int dim, int spacedim>
+void
+test()
+{
+ deallog << "dim= " << dim << ",\t spacedim= " << spacedim << std::endl;
+
+ Triangulation<spacedim> tria0, tria1;
+ Triangulation<2, 3> tria_out;
+ GridOut go;
+ CGALMesh surface_mesh0, surface_mesh1, out_mesh;
+
+ GridGenerator::hyper_ball(tria0, {0., 0., 0.}, 0.4);
+ GridGenerator::hyper_ball(tria1, {0.3, 0.3, 0.3}, 0.4);
+ tria0.refine_global(3);
+ tria1.refine_global(3);
+
+ // Move to CGAL surfaces
+ 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);
+
+ CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh0);
+ CGAL::Polygon_mesh_processing::triangulate_faces(surface_mesh1);
+
+ compute_boolean_operation(surface_mesh0,
+ surface_mesh1,
+ BooleanOperation::compute_union,
+ out_mesh);
+ AdditionalData<3> data;
+ data.edge_size = .02;
+ data.facet_angle = 25;
+ data.facet_size = 0.05;
+ data.facet_distance = 0.05;
+
+ remesh_surface<CGALPoint>(out_mesh, data);
+
+ // Now back to deal.II
+ cgal_surface_mesh_to_dealii_triangulation(out_mesh, tria_out);
+ std::ofstream out_name_spheres("remeshed_union_spheres.vtk");
+ go.write_vtk(tria_out, out_name_spheres);
+ deallog << "OK" << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+ test<3, 3>();
+}
projects / dealii.git / commitdiff
