Mesh_2: Make lloyd_optimize_mesh_2() deterministic (#8768)

## Summary of Changes

Replace a `std::set<Vertex_handle>` by a
`std::set<std::pair<std::size_t,Vertex_handle>>` to make traversing the
set deterministic.
At the same time, do not `erase()` a key, but an `iterator`. 

## Release Management

* Affected package(s): Mesh_2
* Issue(s) solved (if any): fix #8719
* License and copyright ownership: unchanged
This commit is contained in:
Sebastien Loriot 2025-03-28 18:04:34 +01:00 committed by GitHub
commit 3fa6e77002
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
2 changed files with 208 additions and 6 deletions

View File

@ -61,7 +61,14 @@ class Mesh_global_optimizer_2
typedef typename Gt::Vector_2 Vector_2;
typedef typename std::vector<Face_handle> Face_vector;
typedef typename std::set<Vertex_handle> Vertex_set;
typedef std::pair<std::size_t,Vertex_handle> IndexVertexPair;
struct IVP_less {
bool operator()(const IndexVertexPair& ivp1, const IndexVertexPair& ivp2) const
{
return ivp1.first < ivp2.first;
}
};
typedef typename std::set<IndexVertexPair,IVP_less> Vertex_set;
typedef typename std::list<FT> FT_list;
typedef typename std::pair<Vertex_handle,Point_2> Move;
@ -114,13 +121,14 @@ public:
// Fill set containing moving vertices
Vertex_set moving_vertices;
std::size_t ind = 0;
for(typename Tr::Finite_vertices_iterator
vit = cdt_.finite_vertices_begin();
vit != cdt_.finite_vertices_end();
++vit )
{
if(!cdt_.are_there_incident_constraints(vit))
moving_vertices.insert(vit);
moving_vertices.insert(std::make_pair(ind++, vit));
}
double initial_vertices_nb = static_cast<double>(moving_vertices.size());
@ -236,11 +244,12 @@ private:
std::fill(big_moves_.begin(), big_moves_.end(), FT(0));
// Get move for each moving vertex
for ( typename Vertex_set::const_iterator vit = moving_vertices.begin() ;
vit != moving_vertices.end() ; )
for ( typename Vertex_set::iterator vit = moving_vertices.begin() ;
vit != moving_vertices.end() ;)
{
Vertex_handle oldv = *vit;
Vertex_handle oldv = vit->second;
Vector_2 move = compute_move(oldv);
typename Vertex_set::iterator old_vit = vit;
++vit;
if ( CGAL::NULL_VECTOR != move )
@ -249,7 +258,7 @@ private:
moves.push_back(std::make_pair(oldv, new_position));
}
else if(sq_freeze_ratio_ > 0.) //freezing ON
moving_vertices.erase(oldv);
moving_vertices.erase(old_vit);
// Stop if time_limit_ is reached
if ( is_time_limit_reached() )

View File

@ -0,0 +1,193 @@
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/lloyd_optimize_mesh_2.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_mesh_vertex_base_2<K> Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds> CDT;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;
typedef CDT::Vertex_handle Vertex_handle;
typedef CDT::Point Point;
using namespace std;
std::string result;
void CGAL_Remesh(double* vert_xyz_array, size_t vert_count, double criteria_a, double criteria_b, int iteration_number, double*& newVertices, int*& vCount, int*& newFaces, int*& fCount, int fileNum)
{
CDT cdt;
vector<Vertex_handle> cdt_Vh_Boundary;
for (size_t i = 0; i < vert_count; ++i)
{
Vertex_handle vh = cdt.insert(Point(vert_xyz_array[3 * i + 0], vert_xyz_array[3 * i + 1]));
cdt_Vh_Boundary.push_back(vh);
}
// insert Constrain
for (size_t i = 0; i < cdt_Vh_Boundary.size() - 1; ++i)
{
cdt.insert_constraint(cdt_Vh_Boundary[i], cdt_Vh_Boundary[i + 1]);
}
cdt.insert_constraint(cdt_Vh_Boundary[cdt_Vh_Boundary.size() - 1], cdt_Vh_Boundary[0]);
// refine and optimize mesh
Mesher mesher(cdt);
mesher.set_criteria(Criteria(criteria_a, criteria_b));
mesher.refine_mesh();
CGAL::lloyd_optimize_mesh_2(cdt, CGAL::parameters::max_iteration_number = iteration_number);
// make index pair
vector <CDT::Vertex_handle> visitedVertices; // collect visited vertices
map <CDT::Vertex_handle, int> indexList; // create a map to note the index
int i = 0;
for (CDT::Vertex_iterator v_it = cdt.vertices_begin(); v_it != cdt.vertices_end(); ++v_it)
{
CDT::Vertex_handle vh = v_it;
indexList[vh] = i;
visitedVertices.push_back(vh);
i++;
}
// Convert data into double array
int vNum = cdt.number_of_vertices();
newVertices = new double[vNum * 3];
i = 0;
for (CDT::Vertex_iterator vi = cdt.vertices_begin(); vi != cdt.vertices_end(); ++vi)
{
newVertices[i] = vi->point()[0];
i += 1;
newVertices[i] = vi->point()[1];
i += 1;
newVertices[i] = 0;
i += 1;
}
int vertexCount = vNum;
vCount = &vertexCount;
int num_face_in_domain = 0;
for (CDT::Face_iterator f_it = cdt.faces_begin(); f_it != cdt.faces_end(); ++f_it)
{
CDT::Face_handle face = f_it;
if (face->is_in_domain())
{
num_face_in_domain++;
}
}
newFaces = new int[num_face_in_domain * 3];
i = 0;
for (CDT::Face_iterator f_it = cdt.faces_begin(); f_it != cdt.faces_end(); ++f_it)
{
CDT::Face_handle face = f_it;
if (face->is_in_domain())
{
newFaces[i] = int(indexList.find(face->vertex(0))->second);
i += 1;
newFaces[i] = int(indexList.find(face->vertex(1))->second);
i += 1;
newFaces[i] = int(indexList.find(face->vertex(2))->second);
i += 1;
}
}
int faceCount = num_face_in_domain;
fCount = &faceCount;
// print
std::stringstream outputFile;
if (fCount && newFaces) {
outputFile << "\nRemeshed faces (" << *fCount << "):\n";
for (int i = 0; i < (*fCount) * 3; i += 3) {
outputFile << newFaces[i] << " " << newFaces[i + 1] << " " << newFaces[i + 2] << "\n";
}
}
if(fileNum == 1)
{
result = outputFile.str();
}
else
{
assert(result == outputFile.str());
}
}
int main()
{
//
double vert_xyz_array[] = {
3375.4981, 1935.35224056, 0.0,
3350.77259333, 2066.38188194, 0.0,
3210.83712383, 2054.53004190, 0.0,
3068.88, 2060.98161842, 0.0,
3034.24939361, 2066.55369658, 0.0,
3025.6776, 2008.40156297, 0.0,
3013.23241519, 1927.9864, 0.0,
3033.36312437, 1924.87291062, 0.0,
3078.68871988, 1917.86131994, 0.0,
3124.01437021, 1910.85008365, 0.0,
3167.66255113, 1908.86629111, 0.0,
3169.83178711, 1908.76770020, 0.0,
3215.64920401, 1906.68531674, 0.0,
3260.96808047, 1913.74020504, 0.0,
3306.03738982, 1922.24487242, 0.0,
3351.10669914, 1930.74953992, 0.0
};
size_t vert_count = 16;
double criteria_a = 0.125;
double criteria_b = 36.691771392;
int iteration_number = 20;
double* newVertices = nullptr;
int* vCount = nullptr;
int* newFaces = nullptr;
int* fCount = nullptr;
for (int i = 1; i <= 2; ++i)
{
CGAL_Remesh(vert_xyz_array, vert_count, criteria_a, criteria_b, iteration_number, newVertices, vCount, newFaces, fCount, i);
}
std::cout << "\nDone";
return 0;
}