1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include "GmshConfig.h"
7
8 #if defined(HAVE_REVOROPT)
9
10 #include "CVTRemesh.h"
11
12 #include "Revoropt/Mesh/builder_def.hpp"
13 #include "Revoropt/Mesh/sampling_def.hpp"
14 #include "Revoropt/Mesh/normals_def.hpp"
15
16 #include "Revoropt/RVD/rvd.hpp"
17 #include "Revoropt/RVD/rdt.hpp"
18
19 #include "Revoropt/CVT/minimizer.hpp"
20
21 #include "Revoropt/Solver/alglib_lbfgs.hpp"
22
23 #include "GModel.h"
24 #include "MTriangle.h"
25 #include "discreteFace.h"
26
27 #include <vector>
28 #include <iostream>
29
30 StringXNumber CVTRemeshOptions_Number[] = {
31 {GMSH_FULLRC, "Sites", nullptr, 20000.},
32 {GMSH_FULLRC, "Iterations", nullptr, 20.},
33 {GMSH_FULLRC, "Anisotropy", nullptr, 0.03},
34 {GMSH_FULLRC, "Variable density", nullptr, 0.3},
35 {GMSH_FULLRC, "Feature sensitivity", nullptr, 5.},
36 {GMSH_FULLRC, "Normal computation radius", nullptr, 0.005}};
37
38 extern "C" {
GMSH_RegisterCVTRemeshPlugin()39 GMSH_Plugin *GMSH_RegisterCVTRemeshPlugin()
40 {
41 return new GMSH_CVTRemeshPlugin();
42 }
43 }
44
getHelp() const45 std::string GMSH_CVTRemeshPlugin::getHelp() const
46 {
47 return "Plugin(CVTRemesh) triangulates an input geometry using"
48 "centroïdal Voronoï tesselations. The STL mesh of the geometry"
49 "is generated and randomly sampled. An objective function derived"
50 "from centroïdal Voronoï tesselations is then defined on the"
51 "generated sampling, and optimized through LBFGS to obtain a"
52 "regular sampling of the surface. The triangulation is extracted"
53 "as the restricted Delaunay triangulation of the samples and the"
54 "STL mesh.\n\n"
55 "If `View' < 0, the plugin is run on the current view.\n\n"
56 "Plugin(CVTRemesh) creates one new view.";
57 }
58
getNbOptions() const59 int GMSH_CVTRemeshPlugin::getNbOptions() const
60 {
61 return sizeof(CVTRemeshOptions_Number) / sizeof(StringXNumber);
62 }
63
getOption(int iopt)64 StringXNumber *GMSH_CVTRemeshPlugin::getOption(int iopt)
65 {
66 return &CVTRemeshOptions_Number[iopt];
67 }
68
69 // solver callback
70 class SolverCallback {
71 public:
operator ()(Data * data)72 template <typename Data> void operator()(Data *data)
73 {
74 // Output current iteration data
75 Msg::Info("[CVTRemesh] : iteration %d, objective function value is %f",
76 data->k, data->fx);
77 }
78 };
79
execute(PView * v)80 PView *GMSH_CVTRemeshPlugin::execute(PView *v)
81 {
82 // TODO normalization
83
84 GModel *m = GModel::current();
85
86 std::vector<double> vertices;
87 std::vector<unsigned int> faces;
88
89 unsigned int offset = 0;
90 for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
91 (*it)->buildSTLTriangulation();
92 for(std::size_t i = 0; i < (*it)->stl_vertices_uv.size(); ++i) {
93 GPoint p = (*it)->point((*it)->stl_vertices_uv[i]);
94 vertices.push_back(p.x());
95 vertices.push_back(p.y());
96 vertices.push_back(p.z());
97 }
98 for(std::size_t i = 0; i < (*it)->stl_triangles.size(); ++i) {
99 faces.push_back((*it)->stl_triangles[i] + offset);
100 }
101 offset += (*it)->stl_vertices_uv.size();
102 }
103
104 Revoropt::MeshBuilder<3> mesh;
105 mesh.swap_vertices(vertices);
106 mesh.swap_faces(faces);
107
108 double mesh_center[3];
109 double mesh_scale;
110 Revoropt::normalize_mesh(&mesh, mesh_center, &mesh_scale);
111
112 double nradius = (double)CVTRemeshOptions_Number[5].def;
113
114 // normals
115 std::vector<double> normals(3 * mesh.vertices_size());
116 Revoropt::full_robust_vertex_normals(&mesh, nradius, normals.data());
117
118 // lifted vertices
119 std::vector<double> lifted_vertices(6 * mesh.vertices_size(), 0);
120 for(std::size_t vertex = 0; vertex < mesh.vertices_size(); ++vertex) {
121 std::copy(mesh.vertex(vertex), mesh.vertex(vertex) + 3,
122 lifted_vertices.data() + 6 * vertex);
123 std::copy(normals.data() + 3 * vertex, normals.data() + 3 * vertex + 3,
124 lifted_vertices.data() + 6 * vertex + 3);
125 }
126
127 // setup lifted mesh
128 Revoropt::ROMeshWrapper<3, 6> lifted_mesh(lifted_vertices.data(),
129 lifted_vertices.size() / 6, &mesh);
130
131 // triangle weight factor
132 double twfactor = (double)CVTRemeshOptions_Number[3].def;
133
134 // face ratios
135 std::vector<double> triangle_weights(lifted_mesh.faces_size());
136 if(twfactor > 0) {
137 for(std::size_t f = 0; f < lifted_mesh.faces_size(); ++f) {
138 // vertices of the initial triangle
139 const unsigned int *fverts = mesh.face(f);
140
141 // positions
142 const double *x[3];
143 for(int i = 0; i < 3; ++i) { x[i] = lifted_mesh.vertex(fverts[i]); }
144
145 // ratio
146 double ratio = 1;
147
148 // vectors
149 typedef Eigen::Matrix<double, 3, 1> Vector3;
150
151 Eigen::Map<const Vector3> v0(x[0]);
152 Eigen::Map<const Vector3> v1(x[1]);
153 Eigen::Map<const Vector3> v2(x[2]);
154
155 // triangle frame
156 Vector3 U = (v1 - v0);
157 const double U_len = U.norm();
158 if(U_len > 0) {
159 U /= U_len;
160 Vector3 H = (v2 - v0);
161 H = H - H.dot(U) * U;
162 const double H_len = H.norm();
163 if(H_len > 0) {
164 // we know that the triangle is not flat
165 H /= H_len;
166
167 // gradient of the barycentric weights in the triangle
168 Eigen::Matrix<double, 3, 2> bar_grads;
169 bar_grads(2, 0) = 0;
170 bar_grads(2, 1) = 1 / H_len;
171
172 // gradient norms of every normal component
173 for(int i = 0; i < 2; ++i) {
174 // reference frame for the vertex
175 Eigen::Map<const Vector3> vi0(x[(i + 1) % 3]);
176 Eigen::Map<const Vector3> vi1(x[(i + 2) % 3]);
177 Eigen::Map<const Vector3> vi2(x[i]);
178
179 Vector3 Ui = (vi1 - vi0);
180 Ui /= Ui.norm();
181 Vector3 Hi = (vi2 - vi0);
182 Hi = Hi - Hi.dot(Ui) * Ui;
183 const double Hi_invlen = 1 / Hi.norm();
184 Hi *= Hi_invlen;
185 bar_grads(i, 0) = Hi.dot(U) * Hi_invlen;
186 bar_grads(i, 1) = Hi.dot(H) * Hi_invlen;
187 }
188
189 // gradient of each component of the normal
190 Eigen::Map<const Vector3> n0(x[0] + 3);
191 Eigen::Map<const Vector3> n1(x[1] + 3);
192 Eigen::Map<const Vector3> n2(x[2] + 3);
193
194 Eigen::Matrix<double, 3, 2> n_grads =
195 Eigen::Matrix<double, 3, 2>::Zero();
196
197 n_grads = n0 * bar_grads.row(0);
198 n_grads += n1 * bar_grads.row(1);
199 n_grads += n2 * bar_grads.row(2);
200
201 // maximal gradient norm
202 double g_max = n_grads.row(0).dot(n_grads.row(0));
203 double g_other = n_grads.row(1).dot(n_grads.row(1));
204 g_max = g_max > g_other ? g_max : g_other;
205 g_other = n_grads.row(2).dot(n_grads.row(2));
206 g_max = g_max > g_other ? g_max : g_other;
207
208 if(g_max == g_max) { // prevent nan
209 ratio += g_max;
210 }
211 }
212 }
213 triangle_weights[f] = pow(ratio, twfactor);
214 }
215 }
216
217 // normal factor
218 double nfactor = (double)CVTRemeshOptions_Number[2].def;
219 ;
220
221 // weight the normal component by the provided factor
222 for(std::size_t i = 0; i < lifted_mesh.vertices_size(); ++i) {
223 double *v = lifted_vertices.data() + 6 * i;
224 v[3] *= nfactor;
225 v[4] *= nfactor;
226 v[5] *= nfactor;
227 }
228
229 // number of sites
230 unsigned int nsites = (unsigned int)CVTRemeshOptions_Number[0].def;
231
232 // lifted sites
233 std::vector<double> lifted_sites(6 * nsites);
234 if(twfactor > 0) {
235 Revoropt::generate_random_sites<Revoropt::ROMesh<3, 6> >(
236 &lifted_mesh, nsites, lifted_sites.data(), triangle_weights.data());
237 }
238 else {
239 Revoropt::generate_random_sites<Revoropt::ROMesh<3, 6> >(
240 &lifted_mesh, nsites, lifted_sites.data());
241 }
242
243 // setup the cvt minimizer
244 Revoropt::CVT::DirectMinimizer<Revoropt::ROMesh<3, 6> > cvt;
245 cvt.set_sites(lifted_sites.data(), nsites);
246 cvt.set_mesh(&lifted_mesh);
247 if(twfactor > 0) { cvt.set_triangle_weights(triangle_weights.data()); }
248
249 // setup the callback
250 SolverCallback callback;
251
252 // number of iterations
253 unsigned int niter = (unsigned int)CVTRemeshOptions_Number[1].def;
254 ;
255 unsigned int aniso_niter = std::min<unsigned int>(10, niter);
256
257 // solver status
258 int status = 0;
259
260 // isotropic iterations
261 if(niter > 10) {
262 aniso_niter = std::max(aniso_niter, niter * 10 / 100);
263 cvt.set_anisotropy(1);
264 status =
265 cvt.minimize<Revoropt::Solver::AlgLBFGS>(niter - aniso_niter, &callback);
266 }
267
268 // anisotropic iterations
269 if(niter > 0) {
270 // tangent space anisotropy
271 double tanisotropy = (double)CVTRemeshOptions_Number[4].def;
272 ;
273
274 // anisotropic iterations
275 cvt.set_anisotropy(tanisotropy);
276 status = cvt.minimize<Revoropt::Solver::AlgLBFGS>(aniso_niter, &callback);
277 }
278
279 // rdt
280 std::vector<unsigned int> rdt_triangles;
281 Revoropt::RDTBuilder<Revoropt::ROMesh<3, 6> > build_rdt(rdt_triangles);
282 Revoropt::RVD<Revoropt::ROMesh<3, 6> > rvd;
283 rvd.set_sites(lifted_sites.data(), nsites);
284 rvd.set_mesh(&lifted_mesh);
285 rvd.compute(build_rdt);
286
287 GFace *res_face = new discreteFace(m, m->getMaxElementaryNumber(2) + 1);
288 m->add(res_face);
289
290 // scale back and transfer to gmsh
291 std::vector<MVertex *> m_verts(nsites);
292 for(std::size_t i = 0; i < nsites; ++i) {
293 m_verts[i] =
294 new MVertex(lifted_sites[6 * i] * mesh_scale + mesh_center[0],
295 lifted_sites[6 * i + 1] * mesh_scale + mesh_center[1],
296 lifted_sites[6 * i + 2] * mesh_scale + mesh_center[2]);
297 res_face->addMeshVertex(m_verts[i]);
298 }
299 for(std::size_t i = 0; i < rdt_triangles.size() / 3; ++i) {
300 res_face->addTriangle(new MTriangle(m_verts[rdt_triangles[3 * i]],
301 m_verts[rdt_triangles[3 * i + 1]],
302 m_verts[rdt_triangles[3 * i + 2]]));
303 }
304
305 res_face->setAllElementsVisible(true);
306
307 return v;
308 }
309
310 #endif
311