Merge pull request #7307 from janetournois/Mesh_3-sizing_field_with_aabb_tree-GF

Mesh_3 - Document Sizing_field_with_aabb_tree
This commit is contained in:
Laurent Rineau 2023-05-04 17:22:59 +02:00
commit 628d8ae43d
12 changed files with 262 additions and 119 deletions

View File

@ -53,13 +53,6 @@ CGAL tetrahedral Delaunay refinement algorithm.
- This new package wraps all the existing code that deals with a `MeshComplex_3InTriangulation_3` to describe 3D simplicial meshes, and makes the data structure independent from the tetrahedral mesh generation package.
### [Tetrahedral Mesh Generation](https://doc.cgal.org/5.6/Manual/packages.html#PkgMesh3)
- Added two new named parameters to the named constructor `CGAL::create_labeled_image_mesh_domain()`
for automatic detection and protection
of 1D-curves that lie at the intersection of three or more subdomains,
extracted from labeled images.
### [Shape Detection](https://doc.cgal.org/5.6/Manual/packages.html#PkgShapeDetection) (breaking change, major changes)
- **Breaking change**: The region growing part of the package have been reworked to fix design issues introduced with the handling `FaceGraph` models.
@ -104,6 +97,11 @@ extracted from labeled images.
### [3D Mesh Generation](https://doc.cgal.org/5.6/Manual/packages.html#PkgMesh3)
- Deprecated usage of boost parameters in favor of function named parameters.
- Added two new named parameters to the named constructor `CGAL::create_labeled_image_mesh_domain()`
for automatic detection and protection
of 1D-curves that lie at the intersection of three or more subdomains,
extracted from labeled images.
- Added `CGAL::Sizing_field_with_aabb_tree`, a geometry-aware sizing field for feature edges in polyhedral domains.
### [3D Periodic Mesh Generation](https://doc.cgal.org/5.6/Manual/packages.html#PkgPeriodic3Mesh3)
- Deprecated usage of boost parameters in favor of function named parameters.

View File

@ -1,73 +0,0 @@
namespace CGAL {
/*!
\ingroup PkgMesh3Domains
The class `Mesh_constant_domain_field_3` is a model of concept `MeshDomainField_3`. It provides
a constant field accessible using queries on 3D-points.
The class `Mesh_constant_domain_field_3` can also be customized through `set_size()` operations to become
a piecewise constant field, i.e.\ a sizing field with a constant size on each subpart
of the domain.
\tparam Gt is the geometric traits class. It must match the type `Triangulation::Geom_traits`,
where `Triangulation` is the nested type of the model of `MeshComplex_3InTriangulation_3` used
in the meshing process.
\tparam Index is the type of index of the vertices of the triangulation.
It must match the type `%Index` of the model of `MeshDomain_3` used in the meshing process.
\cgalModels `MeshDomainField_3`
\sa `MeshDomainField_3`
*/
template< typename Gt, typename Index >
class Mesh_constant_domain_field_3 {
public:
/// \name Types
/// @{
/*!
Numerical type.
*/
typedef Gt::FT FT;
/*!
Point type.
*/
typedef Gt::Point_3 Point_3;
/*!
Type of index of the vertices of the triangulation.
*/
typedef Index Index;
/// @}
/// \name Creation
/// @{
/*!
Builds a constant domain field with size `size`.
*/
Mesh_constant_domain_field_3(FT size);
/// @}
/// \name Operations
/// @{
/*!
Sets the size such as `operator()` will return size `size`
at any query point of dimension `dimension` and index `index`.
*/
void set_size(FT size, int dimension, const Index& index);
/// @}
}; /* end Mesh_constant_domain_field_3 */
} /* end namespace CGAL */

View File

@ -46,10 +46,9 @@ on the input feature including the query point.
/// @{
/*!
returns the value of the sizing field at the point `p`,
assumed to be included in the input complex feature with dimension `dimension`
and mesh vertex index `index`.
and mesh subcomplex index `index`.
*/
FT operator()(const Point_3& p, int dimension, const Index& index) const;

View File

@ -7,6 +7,8 @@ INPUT += \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Polyhedral_complex_mesh_domain_3.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_domain_with_polyline_features_3.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_3/generate_label_weights.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Sizing_field_with_aabb_tree.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Mesh_constant_domain_field_3.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/exude_mesh_3.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/odt_optimize_mesh_3.h \
${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/lloyd_optimize_mesh_3.h \

View File

@ -1013,7 +1013,24 @@ Right: the mesh with added 1D-features along sharp edges and edges shared
by at least two of the input polyhedral surfaces forming the complex.
\cgalFigureCaptionEnd
\subsubsection Mesh_3EdgesSizingField Sizing Field for Feature Edges
The following example shows how to generate a mesh from a polyhedral complex
or a polyhedral surface with polyline features.
Polyline features are covered by a set of so-called "protecting balls" which
sizes are highly related to the edge length, and
are driven by the size component of `Mesh_edge_criteria_3`
(see Section \ref Mesh_3Protectionof0and1dimensionalExposed).
The ideal size can be computed using
`Sizing_field_with_aabb_tree` that helps start the feature protection
and one-dimensional meshing process with a good initial guess.
To fit the protecting balls requirements,
no protecting ball can have its radius
larger than half of the distance from the corresponding vertex (its center),
to surface patches the
current polyline feature does not belong to.
\cgalExample{Mesh_3/mesh_polyhedral_domain_with_features_sizing.cpp}
\subsection Mesh_3TuningMeshOptimization Tuning Mesh Optimization

View File

@ -20,6 +20,10 @@
/// \ingroup PkgMesh3Ref
/// The functors in this group perform polyline features detection in input domains.
/// \defgroup PkgMesh3DomainFields Mesh Domain Fields
/// \ingroup PkgMesh3Ref
/// The classes in this group are models of `MeshDomainField_3`
/// \defgroup PkgMesh3Functions Mesh Generation Functions
/// \ingroup PkgMesh3Ref
/// The two main functions to generate a mesh are `make_mesh_3()` and `refine_mesh_3()`.

View File

@ -23,6 +23,7 @@
\example Mesh_3/mesh_polyhedral_complex.cpp
\example Mesh_3/remesh_polyhedral_surface.cpp
\example Mesh_3/mesh_polyhedral_domain_with_features.cpp
\example Mesh_3/mesh_polyhedral_domain_with_features_sizing.cpp
\example Mesh_3/mesh_polyhedral_domain_with_surface_inside.cpp
\example Mesh_3/mesh_polyhedral_domain_with_lipschitz_sizing.cpp
\example Mesh_3/mesh_two_implicit_spheres_with_balls.cpp

View File

@ -96,6 +96,10 @@ create_single_source_cgal_program("mesh_polyhedral_domain_with_features.cpp")
target_link_libraries(mesh_polyhedral_domain_with_features
PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program("mesh_polyhedral_domain_with_features_sizing.cpp")
target_link_libraries(mesh_polyhedral_domain_with_features_sizing
PUBLIC CGAL::Eigen3_support)
create_single_source_cgal_program("mesh_polyhedral_domain_with_features_sm.cpp")
target_link_libraries(mesh_polyhedral_domain_with_features_sm
PUBLIC CGAL::Eigen3_support)
@ -207,6 +211,7 @@ if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support)
mesh_polyhedral_domain
mesh_polyhedral_domain_sm
mesh_polyhedral_domain_with_features
mesh_polyhedral_domain_with_features_sizing
mesh_polyhedral_domain_with_features_sm
mesh_polyhedral_domain_with_lipschitz_sizing
mesh_two_implicit_spheres_with_balls)

View File

@ -0,0 +1,66 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Sizing_field_with_aabb_tree.h>
// Domain
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Surface_mesh<K::Point_3>;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron>;
using Features_sizing_field = CGAL::Sizing_field_with_aabb_tree<K, Mesh_domain>;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain,CGAL::Default,Concurrency_tag>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<
Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_index> C3t3;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
namespace params = CGAL::parameters;
int main(int argc, char*argv[])
{
const std::string fname = (argc>1)?argv[1]:CGAL::data_file_path("meshes/fandisk.off");
std::ifstream input(fname);
Polyhedron polyhedron;
input >> polyhedron;
if(input.fail()){
std::cerr << "Error: Cannot read file " << fname << std::endl;
return EXIT_FAILURE;
}
// Create domain
Mesh_domain domain(polyhedron);
// Get sharp features
domain.detect_features();
// Mesh criteria
Features_sizing_field edges_sizing_field(0.07, domain);
Mesh_criteria criteria(params::edge_size(edges_sizing_field).
facet_distance(0.0072).
cell_radius_edge_ratio(3));
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, params::no_exude().no_perturb());
// Output
CGAL::dump_c3t3(c3t3, "out_sizing_field_with_aabb_tree");
return EXIT_SUCCESS;
}

View File

@ -28,27 +28,77 @@
namespace CGAL {
/*!
* @ingroup PkgMesh3DomainFields
* The class `Mesh_constant_domain_field_3` is a model of concept `MeshDomainField_3`. It provides
* a constant field accessible using queries on 3D-points.
*
* The class `Mesh_constant_domain_field_3` can also be customized through `set_size()` operations to become
* a piecewise constant field, i.e. a sizing field with a constant size on each subpart
* of the domain.
*
* @tparam Gt is the geometric traits class. It must match the type `Triangulation::Geom_traits`,
* where `Triangulation` is the nested type of the model of `MeshComplex_3InTriangulation_3` used
* in the meshing process.
*
* @tparam Index_ is the type of index of the vertices of the triangulation.
* It must match the type `%Index` of the model of `MeshDomain_3` used in the meshing process.
*
* @cgalModels `MeshDomainField_3`
*/
template <typename Gt, typename Index_>
class Mesh_constant_domain_field_3
{
public:
/// \name Types
/// @{
/*!
* Numerical type.
*/
typedef typename Gt::FT FT;
typedef typename boost::mpl::eval_if_c<
/*!
* Point type.
*/
#ifdef DOXYGEN_RUNNING
typedef Gt::Point_3 Point_3;
#else
typedef typename boost::mpl::eval_if_c<
internal::Has_nested_type_Bare_point<Gt>::value,
typename internal::Bare_point_type<Gt>,
boost::mpl::identity<typename Gt::Point_3>
>::type Point_3;
#endif
/*!
* Type of index of the vertices of the triangulation.
*/
typedef Index_ Index;
/// @}
private:
// Map to store field values
typedef std::map<std::pair<int,Index>,FT> Values;
public:
/// Constructor
Mesh_constant_domain_field_3(const FT& d) : d_(d) {}
/// \name Creation
/// @{
/// Returns size
/*!
* Builds a constant domain field with size `size`.
*/
Mesh_constant_domain_field_3(const FT& size) : d_(size) {}
/// @}
/// \name Operations
/// @{
/*!
* Returns the size of query points of dimension `dim` and index `index`.
*/
FT operator()(const Point_3&, const int dim, const Index& index) const
{
typename Values::const_iterator it = values_.find(std::make_pair(dim,index));
@ -57,11 +107,15 @@ public:
return d_;
}
/// Sets size at any point of dimension `dim` and index `index`.
/*!
* Sets the size such that `operator()` for any query point
* of dimension `dim` and index `index` returns `size`.
*/
void set_size(const FT& size, const int dim, const Index& index)
{
values_.insert(std::make_pair(std::make_pair(dim,index),size));
}
/// @}
private:
// default value

View File

@ -17,34 +17,65 @@
#include <CGAL/Profile_counter.h>
#include <CGAL/Delaunay_triangulation_3.h>
#include "Facet_patch_id_map.h"
#include "Get_curve_index.h"
#include <CGAL/Mesh_3/Protect_edges_sizing_field.h> // for weight_modifier
#include <cstddef>
#include <memory>
#include <limits>
#include <CGAL/Mesh_3/experimental/Facet_patch_id_map.h>
#include <CGAL/Mesh_3/experimental/Get_curve_index.h>
#include <CGAL/Mesh_3/experimental/AABB_filtered_projection_traits.h>
#include <boost/container/flat_set.hpp>
#if defined(CGAL_MESH_3_PROTECTION_HIGH_VERBOSITY) || defined(CGAL_NO_ASSERTIONS) == 0
# include <sstream>
# include <boost/format.hpp>
#endif // CGAL_MESH_3_PROTECTION_HIGH_VERBOSITY || (! CGAL_NO_ASSERTIONS)
#include "AABB_filtered_projection_traits.h"
namespace CGAL {
/*!
* @ingroup PkgMesh3DomainFields
*
* The class `Sizing_field_with_aabb_tree` is a model of concept `MeshDomainField_3`.
* It provides a sizing field to be used in the meshing process of a polyhedral surface with features.
*
* At each query point `p`, the field value is the minimum of the distances to
* - the surface patches `p` does not belong to, and
* - for each curve `C` it does belong to, the intersection points of `C` with
* the plane orthogonal to `C` and passing through `p`.
*
* This sizing field is designed to be used for the edges of the mesh complex,
* in `Mesh_edge_criteria_3` or as `edge_size` in `Mesh_criteria_3`.
* Using this sizing field for complex edges provides a high-quality hint
* to the protecting balls placement algorithm, since it ensures that no
* protecting ball will intersect a surface patch to which the
* corresponding vertex does not belong.
*
* @tparam GeomTraits is the geometric traits class. It must match the type `Triangulation::Geom_traits`,
* where `Triangulation` is the nested type of the model of `MeshComplex_3InTriangulation_3` used
* in the meshing process.
* @tparam MeshDomain is the type of the domain. It must be a model of `MeshDomainWithFeatures_3`,
* and derive from `CGAL::Polyhedral_mesh_domain_with_features_3`
template <typename GeomTraits, typename MeshDomain,
typename Input_facets_AABB_tree = typename MeshDomain::AABB_tree,
typename Get_curve_index_ = CGAL::Default,
typename Facet_patch_id_map_ = CGAL::Default
* @cgalModels `MeshDomainField_3`
*/
template <typename GeomTraits,
typename MeshDomain
#ifndef DOXYGEN_RUNNING
, typename Input_facets_AABB_tree_ = CGAL::Default
, typename Get_curve_index_ = CGAL::Default
, typename Facet_patch_id_map_ = CGAL::Default
#endif
>
struct Sizing_field_with_aabb_tree
{
typedef GeomTraits Kernel_;
typedef typename Kernel_::FT FT;
typedef typename Kernel_::Point_3 Point_3;
using FT = typename GeomTraits::FT;
using Point_3 = typename GeomTraits::Point_3;
using Index = typename MeshDomain::Index;
typedef typename MeshDomain::Index Index;
private:
typedef typename MeshDomain::Corner_index Corner_index;
typedef typename MeshDomain::Curve_index Curve_index;
typedef typename MeshDomain::Surface_patch_index Patch_index;
@ -58,28 +89,52 @@ struct Sizing_field_with_aabb_tree
typedef std::vector<Patches_ids> Corners_incident_patches;
typedef std::vector<Patches_ids> Curves_incident_patches;
typedef GeomTraits Kernel_;
typedef CGAL::Delaunay_triangulation_3<Kernel_> Dt;
typedef Input_facets_AABB_tree Input_facets_AABB_tree_;
typedef typename Input_facets_AABB_tree_::AABB_traits AABB_traits;
using Point_and_primitive_id = typename AABB_traits::Point_and_primitive_id;
typedef typename Input_facets_AABB_tree_::Primitive Input_facets_AABB_tree_primitive_;
typedef typename MeshDomain::Curves_AABB_tree Input_curves_AABB_tree_;
using Input_facets_AABB_tree = typename CGAL::Default::Get<
Input_facets_AABB_tree_,
typename MeshDomain::AABB_tree
>::type;
using AABB_traits = typename Input_facets_AABB_tree::AABB_traits;
using Point_and_primitive_id = typename AABB_traits::Point_and_primitive_id;
typedef typename Input_facets_AABB_tree::Primitive Input_facets_AABB_tree_primitive_;
typedef typename MeshDomain::Curves_AABB_tree Input_curves_AABB_tree_;
typedef typename Input_curves_AABB_tree_::Primitive Input_curves_AABB_tree_primitive_;
typedef typename CGAL::Default::Get<
using Get_curve_index = typename CGAL::Default::Get<
Get_curve_index_,
CGAL::Mesh_3::Get_curve_index<typename MeshDomain::Curves_AABB_tree::Primitive>
>::type Get_curve_index;
typedef typename CGAL::Default::Get<
>::type;
using Facet_patch_id_map = typename CGAL::Default::Get<
Facet_patch_id_map_,
CGAL::Mesh_3::Facet_patch_id_map<MeshDomain,
typename Input_facets_AABB_tree::Primitive>
>::type Facet_patch_id_map;
>::type;
public:
/// \name Creation
/// @{
/*!
* Constructor for the sizing field.
* @param d the maximal value returned by `operator()`, corresponding
* to an upper bound on feature edge lengths in the output mesh.
* @param domain the mesh domain to be meshed
*/
Sizing_field_with_aabb_tree(const FT& d, const MeshDomain& domain)
: Sizing_field_with_aabb_tree(d,
domain,
domain.aabb_tree(),
Get_curve_index(),
Facet_patch_id_map())
{}
/// @}
struct Face_ids_traversal_traits {
using Limits = std::numeric_limits<Patch_index>;
Patch_index min{Limits::max()};
Patch_index min{(Limits::max)()};
Patch_index max{Limits::lowest()};
Facet_patch_id_map index_map{};
@ -93,16 +148,17 @@ struct Sizing_field_with_aabb_tree
}
};
#ifndef DOXYGEN_RUNNING
Sizing_field_with_aabb_tree
(typename Kernel_::FT d,
const Input_facets_AABB_tree_& aabb_tree,
(const typename Kernel_::FT d,
const MeshDomain& domain,
const Input_facets_AABB_tree& aabb_tree,
Get_curve_index get_curve_index = Get_curve_index(),
Facet_patch_id_map facet_patch_id_map = Facet_patch_id_map()
)
: d_ptr{std::make_shared<Private_data>(d,
aabb_tree,
domain,
aabb_tree,
get_curve_index,
facet_patch_id_map)}
{
@ -258,7 +314,17 @@ struct Sizing_field_with_aabb_tree
result = projection_traits.closest_point_and_primitive();
return result;
}
#endif // DOXYGEN_RUNNING
public:
/// \name Operations
/// @{
/*!
* returns the value of the sizing field at the point `p`,
* assumed to be included in the input complex feature with dimension `dimension`
* and mesh subcomplex index `id`.
*/
double operator()(const Point_3& p,
const int dim,
const Index& id) const
@ -597,23 +663,27 @@ struct Sizing_field_with_aabb_tree
#endif // CGAL_MESH_3_PROTECTION_HIGH_VERBOSITY
return result;
}
/// @}
private:
using Kd_tree = CGAL::AABB_search_tree<AABB_traits>;
struct Private_data {
using FT = typename Kernel_::FT;
Private_data(FT d, const Input_facets_AABB_tree_& aabb_tree,
Private_data(FT d,
const MeshDomain& domain,
const Input_facets_AABB_tree& aabb_tree,
Get_curve_index get_curve_index,
Facet_patch_id_map facet_patch_id_map)
: d_(d)
, aabb_tree(aabb_tree)
, domain(domain)
, aabb_tree(aabb_tree)
, get_curve_index(get_curve_index)
, facet_patch_id_map(facet_patch_id_map)
{}
FT d_;
const Input_facets_AABB_tree_& aabb_tree;
const MeshDomain& domain;
const Input_facets_AABB_tree& aabb_tree;
Get_curve_index get_curve_index;
Facet_patch_id_map facet_patch_id_map;
Dt dt{};
@ -628,4 +698,6 @@ private:
std::shared_ptr<Private_data> d_ptr;
};
}//end namespace CGAL
#endif // CGAL_MESH_3_SIZING_FIELD_WITH_AABB_TREE_H

View File

@ -244,7 +244,7 @@ edge_criteria(double b, Mesh_fnt::Domain_tag)
return Edge_criteria(b);
}
#include <CGAL/Mesh_3/experimental/Sizing_field_with_aabb_tree.h>
#include <CGAL/Sizing_field_with_aabb_tree.h>
#include <CGAL/Mesh_3/experimental/Facet_topological_criterion_with_adjacency.h>
template < typename D_, typename Tag >
@ -254,7 +254,7 @@ edge_criteria(double edge_size, Mesh_fnt::Polyhedral_domain_tag)
{
if(p_.use_sizing_field_with_aabb_tree) {
typedef typename Domain::Surface_patch_index_set Set_of_patch_ids;
typedef Sizing_field_with_aabb_tree<Kernel, Domain> Mesh_sizing_field; // type of sizing field for 0D and 1D features
typedef CGAL::Sizing_field_with_aabb_tree<Kernel, Domain> Mesh_sizing_field; // type of sizing field for 0D and 1D features
typedef std::vector<Set_of_patch_ids> Patches_ids_vector;
typedef typename Domain::Curve_index Curve_index;
const Curve_index max_index = domain_->maximal_curve_index();
@ -265,9 +265,7 @@ edge_criteria(double edge_size, Mesh_fnt::Polyhedral_domain_tag)
(*patches_ids_vector_p)[curve_id] = domain_->get_incidences(curve_id);
}
Mesh_sizing_field* sizing_field_ptr =
new Mesh_sizing_field(edge_size,
domain_->aabb_tree(),
*domain_);
new Mesh_sizing_field(edge_size, *domain_, domain_->aabb_tree());
// The sizing field object, as well as the `patch_ids_vector` are
// allocated on the heap, and the following `boost::any` object,
// containing two shared pointers, is used to make the allocated