Merge branch 'Tet_remeshing-flips_on_surface-jtournois' into Tet_remeshing-wip-jtournois

# Conflicts:
#	Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp
#	Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h
#	Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h
This commit is contained in:
Jane Tournois 2024-02-02 10:38:06 +01:00
commit 72fb087bd5
5 changed files with 57 additions and 91 deletions

View File

@ -12,27 +12,27 @@
#include <CGAL/IO/File_medit.h>
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Mesh_polyhedron_3<K>::type Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K> Mesh_domain;
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron = CGAL::Mesh_polyhedron_3<K>::type;
using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3<K>;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
using Concurrency_tag = CGAL::Parallel_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
using Concurrency_tag = CGAL::Sequential_tag;
#endif
// Triangulation for Meshing
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;
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
typedef CGAL::Triangulation_3<typename Tr::Geom_traits,
typename Tr::Triangulation_data_structure> Triangulation_3;
using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
using Vertex_handle = Triangulation_3::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
@ -68,6 +68,7 @@ int main(int argc, char* argv[])
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Property map of constraints
Constraints_set constraints;
Constraints_pmap constraints_pmap(constraints);

View File

@ -1,6 +1,3 @@
#define CGAL_MESH_3_VERBOSE 1
#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
@ -13,52 +10,42 @@
#include <CGAL/make_mesh_3.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/IO/File_binary_mesh_3.h>
#include <boost/container/flat_set.hpp>
int nb_surface_flip_candidates = 0;
int nb_surface_flip_done = 0;
int nb_surface_44_configs = 0;
int nb_surface_nm_configs = 0;
int nb_surface_44_flips_done = 0;
int nb_surface_nm_flips_done = 0;
std::ostringstream oss_flip("flipped_surface_edges.polylines.txt");
// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Polyhedron;
typedef CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron> Mesh_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>;
#ifdef CGAL_CONCURRENT_MESH_3
typedef CGAL::Parallel_tag Concurrency_tag;
using Concurrency_tag = CGAL::Parallel_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
using Concurrency_tag = CGAL::Sequential_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;
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain, CGAL::Default, Concurrency_tag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<
Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>;
// Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
// Triangulation for Remeshing
using Triangulation_3 = CGAL::Triangulation_3<Tr::Geom_traits,
Tr::Triangulation_data_structure>;
using Vertex_handle = Triangulation_3::Vertex_handle;
using Vertex_pair = std::pair<Vertex_handle, Vertex_handle>;
using Constraints_set = std::unordered_set<Vertex_pair, boost::hash<Vertex_pair>>;
using Constraints_pmap = CGAL::Boolean_property_map<Constraints_set>;
// To avoid verbose function and named parameters call
using namespace CGAL::parameters;
int main(int argc, char* argv[])
{
//"data/Shape1/shape1.off" 40. "data/Shape1/mesh_40.binary.cgal" "data/Shape1/remesh_40.binary.cgal"
const char* fname = (argc > 1) ? argv[1] : "data/tensileASCII.off";
const double target = (argc > 2) ? atof(argv[2]) : 200.;
const int nb_iter = (argc > 3) ? atoi(argv[3]) : 10;
const char* fname_mesh3 = (argc > 4) ? argv[4] : "data/Tensile/mesh_200.binary.cgal";
const char* fname_remesh = (argc > 5) ? argv[5] : "data/Tensile/remesh_200.binary.cgal";
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/anchor.off");
const int nb_iter = (argc > 2) ? atoi(argv[2]) : 5;
std::ifstream input(fname);
Polyhedron polyhedron;
@ -83,10 +70,7 @@ int main(int argc, char* argv[])
domain.detect_features();
// Mesh criteria
double CellSize = 31.3082;
int mesh_factor = 50;
double size = (CellSize * mesh_factor) / 100;
//const double size = 40.;
const double size = 0.072;
Mesh_criteria criteria(edge_size = size,
facet_angle = 25,
facet_size = size,
@ -94,37 +78,21 @@ int main(int argc, char* argv[])
cell_size = size);
// Mesh generation
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);// , no_perturb(), no_exude());
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria, no_perturb(), no_exude());
// Output
// std::ofstream os(fname_mesh3, std::ios_base::out | std::ios_base::binary);
// CGAL::Mesh_3::save_binary_file(os, c3t3);
CGAL::dump_c3t3(c3t3, "out_after_meshing");
Constraints_set constraints;
Constraints_pmap constraints_pmap(constraints);
const double min_dihedral_angle = CGAL::Tetrahedral_remeshing::min_dihedral_angle(c3t3);
std::cout << "Min dihedral angle after meshing = " << min_dihedral_angle << std::endl;
Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3),
CGAL::parameters::edge_is_constrained_map(constraints_pmap));
// Remeshing
CGAL::tetrahedral_isotropic_remeshing(c3t3, size,
CGAL::tetrahedral_isotropic_remeshing(tr, size,
CGAL::parameters::number_of_iterations(nb_iter));
// .smooth_constrained_edges(true));
// Output
// std::ofstream osr(fname_remesh, std::ios_base::out | std::ios_base::binary);
// CGAL::Mesh_3::save_binary_file(osr, c3t3);
CGAL::dump_c3t3(c3t3, "out_after_remeshing");
std::cout << "Meshing and remeshing done." << std::endl;
std::cout << "Surface flip attempts : " << nb_surface_flip_candidates << std::endl;
std::cout << "Surface flip done : " << nb_surface_flip_done << std::endl;
std::cout << "Surface 4/4 flips candidates proportion : "
<< 100.*nb_surface_44_configs / (double)(nb_surface_nm_configs + nb_surface_44_configs)
<< " percent" << std::endl;
std::cout << "nb surface flips 4/4 done = " << nb_surface_44_flips_done << std::endl;
std::cout << "nb surface flips n/m done = " << nb_surface_nm_flips_done << std::endl;
std::ofstream out("out_remeshed.mesh");
CGAL::IO::write_MEDIT(out, tr);
out.close();
return EXIT_SUCCESS;
}

View File

@ -1169,11 +1169,9 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge,
if (are_edge_lengths_valid(edge, c3t3, collapse_type, new_pos, sizing, cell_selector/*, adaptive = false*/)
&& collapse_preserves_surface_star(edge, c3t3, new_pos, cell_selector))
{
CGAL_assertion_code(typename Tr::Cell_handle dc);
CGAL_assertion_code(int di);
CGAL_assertion_code(int dj);
CGAL_assertion(c3t3.triangulation().is_edge(edge.first->vertex(edge.second),
edge.first->vertex(edge.third), dc, di, dj));
CGAL_assertion(c3t3.triangulation().tds().is_edge(
edge.first->vertex(edge.second),
edge.first->vertex(edge.third)));
Vertex_handle v0_init = edge.first->vertex(edge.second);
Vertex_handle v1_init = edge.first->vertex(edge.third);

View File

@ -424,9 +424,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3,
indices(facet_circulator->second, j));
if (curr != vh0 && curr != vh1)
{
Cell_handle ch;
int i0, i1;
if (tr.is_edge(curr, vh, ch, i0, i1))
if (tr.tds().is_edge(curr, vh))
is_edge = true;
}
}
@ -1864,7 +1862,6 @@ std::size_t flipBoundaryEdges(
Tr& tr = c3t3.triangulation();
std::vector<Edge_vv> candidate_edges_for_flip;
for (const Edge& e : boundary_edges)
{
if (!c3t3.is_in_complex(e))
@ -1919,14 +1916,7 @@ std::size_t flipBoundaryEdges(
const Vertex_handle vh2 = third_vertex(f0, vh0, vh1, tr);
const Vertex_handle vh3 = third_vertex(f1, vh0, vh1, tr);
CGAL_assertion_code(int li);
CGAL_assertion_code(int lj);
CGAL_assertion_code(int lk);
CGAL_assertion_code(Cell_handle cc);
CGAL_assertion(tr.is_facet(vh0, vh1, vh2, cc, li, lj, lk));
CGAL_assertion(c3t3.is_in_complex(cc, (6 - li - lj - lk)));
CGAL_assertion(tr.is_facet(vh0, vh1, vh3, cc, li, lj, lk));
CGAL_assertion(c3t3.is_in_complex(cc, (6 - li - lj - lk)));
CGAL_assertion_code(debug::check_facets(vh0, vh1, vh2, vh3, c3t3));
if (!tr.tds().is_edge(vh2, vh3)) // most-likely to happen early exit
{

View File

@ -97,7 +97,16 @@ third_vertex(const typename Tr::Facet& f,
typename Tr::Vertex_handle();
}
// returns angle in degrees
{
for(auto v : tr.vertices(f))
if(v != v0 && v != v1)
return v;
CGAL_assertion(false);
return
typename Tr::Vertex_handle();
}
template<typename Gt, typename Point>
typename Gt::FT dihedral_angle(const Point& p,
const Point& q,