GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <geogram/mesh/mesh_manifold_harmonics.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/basic/file_system.h>
44			#include <geogram/basic/geometry_nd.h>
45			#include <geogram/bibliography/bibliography.h>
46			#include <geogram/NL/nl.h>
47
48			namespace {
49			using namespace GEO;
50
51			/**
52			* \brief Computes a coefficient of the P1 Laplacian.
53			* \param[in] M a reference to a surface mesh
54			* \param[in] f facet index
55			* \param[in] v1 , v2 global indices of two vertices of \p f
56			* \return the cotangent of the angle at the corner of \p f opposite to
57			* \p v1 and \p v2
58			*/
59		✗	inline double P1_FEM_coefficient(
60			const Mesh& M, index_t f, index_t v1, index_t v2
61			) {
62		✗	index_t v3 = NO_VERTEX;
63		✗	for(index_t lv=0; lv<M.facets.nb_vertices(f); ++lv) {
64		✗	index_t v = M.facets.vertex(f,lv);
65		✗	if(v != v1 && v != v2) {
66		✗	v3 = v;
67		✗	break;
68			}
69			}
70		✗	geo_assert(v3 != NO_VERTEX);
71
72			// cotan weights, in arbitrary dimension
73		✗	const double* p1 = M.vertices.point_ptr(v1);
74		✗	const double* p2 = M.vertices.point_ptr(v2);
75		✗	const double* p3 = M.vertices.point_ptr(v3);
76
77		✗	double Lu = 0.0;
78		✗	double Lv = 0.0;
79		✗	double cosangle = 0.0;
80		✗	for(index_t d=0; d<M.vertices.dimension(); ++d) {
81		✗	double u = p2[d] - p1[d];
82		✗	double v = p3[d] - p1[d];
83		✗	Lu += u*u;
84		✗	Lv += v*v;
85		✗	cosangle += u*v;
86			}
87		✗	double Luv = ::sqrt(Lu*Lv);
88		✗	if(Luv < 1e-50) {
89		✗	cosangle = 1.0;
90			} else {
91		✗	cosangle /= Luv;
92			}
93		✗	geo_clamp(cosangle, -1.0, 1.0);
94		✗	return 1.0 / ::tan(::acos(cosangle));
95			}
96
97
98			/**
99			* \brief Assemble the stiffness and mass matrices
100			* of the Laplacian in OpenNL.
101			* \details This function is supposed to be called
102			* between nlBegin(NL_SYSTEM) and nlEnd().
103			* \param[in] M a const reference to a surface mesh.
104			* \param[in] discretization the discretization of the Laplace-Beltrami
105			* operator, one of:
106			* - COMBINATORIAL: 1.0 everywhere
107			* - UNIFORM: combinatorial divided by node degree
108			* - FEM_P1: linear finite elements
109			* - FEM_P1_LUMPED: linear finite elements with lumped mass matrix
110			*/
111		✗	void assemble_Laplacian_matrices(
112			const Mesh& M,
113			LaplaceBeltramiDiscretization discretization
114			) {
115
116			// Step 1: compute vertices degrees (used by
117			// uniform weights).
118			// **************************************************
119
120		✗	vector<index_t> v_degree;
121		✗	if(discretization == UNIFORM) {
122		✗	v_degree.assign(M.vertices.nb(), 0);
123		✗	for(index_t c: M.facet_corners) {
124		✗	index_t v = M.facet_corners.vertex(c);
125		✗	++v_degree[v];
126			}
127			}
128
129			// Sum of row coefficient associated with each vertex
130		✗	vector<double> v_row_sum(M.vertices.nb(), 0.0);
131
132			// Step 2: compute stiffness matrix
133			// **************************************************
134
135		✗	nlMatrixMode(NL_STIFFNESS_MATRIX);
136		✗	nlBegin(NL_MATRIX);
137
138		✗	for(index_t f: M.facets) {
139		✗	index_t fnv = M.facets.nb_vertices(f);
140		✗	for(index_t lv=0; lv<fnv; ++lv) {
141		✗	index_t v1 = M.facets.vertex(f,lv);
142		✗	index_t v2 = M.facets.vertex(f,(lv+1)%fnv);
143		✗	switch(discretization) {
144		✗	case COMBINATORIAL: {
145		✗	double w = 1.0;
146		✗	nlAddIJCoefficient(v1,v2,w);
147		✗	v_row_sum[v1] += w;
148		✗	} break;
149		✗	case UNIFORM: {
150		✗	double w = 1.0 / double(v_degree[v1]);
151		✗	nlAddIJCoefficient(v1,v2,w);
152		✗	v_row_sum[v1] += w;
153		✗	} break;
154		✗	case FEM_P1:
155			case FEM_P1_LUMPED: {
156		✗	double w = 0.5 * P1_FEM_coefficient(M,f,v1,v2);
157		✗	nlAddIJCoefficient(v1,v2,w);
158		✗	nlAddIJCoefficient(v2,v1,w);
159		✗	v_row_sum[v1] += w;
160		✗	v_row_sum[v2] += w;
161		✗	} break;
162			}
163			}
164			}
165		✗	for(index_t v: M.vertices) {
166			// Diagonal term is minus row sum
167			// plus small number to make M non-singular
168		✗	nlAddIJCoefficient(v,v,-v_row_sum[v] + 1e-6);
169			}
170		✗	nlEnd(NL_MATRIX);
171
172			// Step 3: compute mass matrix
173			// **************************************************
174
175		✗	if(discretization == FEM_P1 || discretization == FEM_P1_LUMPED) {
176		✗	nlMatrixMode(NL_MASS_MATRIX);
177		✗	nlBegin(NL_MATRIX);
178		✗	for(index_t f: M.facets) {
179		✗	index_t v1 = M.facets.vertex(f,0);
180		✗	index_t v2 = M.facets.vertex(f,1);
181		✗	index_t v3 = M.facets.vertex(f,2);
182		✗	const double* p1 = M.vertices.point_ptr(v1);
183		✗	const double* p2 = M.vertices.point_ptr(v2);
184		✗	const double* p3 = M.vertices.point_ptr(v3);
185		✗	double A = Geom::triangle_area(
186		✗	p1,p2,p3, coord_index_t(M.vertices.dimension())
187			);
188
189		✗	if(discretization == FEM_P1_LUMPED) {
190
191		✗	nlAddIJCoefficient(v1,v1,A/3.0);
192		✗	nlAddIJCoefficient(v2,v2,A/3.0);
193		✗	nlAddIJCoefficient(v3,v3,A/3.0);
194
195		✗	} else if(discretization == FEM_P1) {
196
197		✗	nlAddIJCoefficient(v1,v2,A/12.0);
198		✗	nlAddIJCoefficient(v1,v3,A/12.0);
199		✗	nlAddIJCoefficient(v2,v3,A/12.0);
200		✗	nlAddIJCoefficient(v2,v1,A/12.0);
201		✗	nlAddIJCoefficient(v3,v1,A/12.0);
202		✗	nlAddIJCoefficient(v3,v2,A/12.0);
203
204		✗	nlAddIJCoefficient(v1,v1,A/6.0);
205		✗	nlAddIJCoefficient(v2,v2,A/6.0);
206		✗	nlAddIJCoefficient(v3,v3,A/6.0);
207			}
208			}
209		✗	nlEnd(NL_MATRIX);
210			}
211		✗	}
212			}
213
214
215			namespace GEO {
216
217
218
219		✗	void mesh_compute_manifold_harmonics(
220			Mesh& M, index_t nb_eigens,
221			LaplaceBeltramiDiscretization discretization,
222			const std::string& attribute_name,
223			double shift,
224			bool print_spectrum
225			) {
226
227		✗	geo_cite("DBLP:conf/smi/Levy06");
228		✗	geo_cite("DBLP:journals/cgf/ValletL08");
229
230		✗	if(M.vertices.attributes().is_defined(attribute_name)) {
231		✗	M.vertices.attributes().delete_attribute_store(attribute_name);
232			}
233
234			// Step 1: configure eigen solver
235			// **************************************************
236
237
238		✗	if(!nlInitExtension("ARPACK")) {
239		✗	Logger::err("MH")
240		✗	<< "Could not initialize OpenNL ARPACK extension"
241		✗	<< std::endl;
242		✗	return;
243			}
244
245		✗	nlNewContext();
246
247		✗	nlEigenSolverParameteri(NL_EIGEN_SOLVER, NL_ARPACK_EXT);
248		✗	nlEigenSolverParameteri(NL_NB_VARIABLES, NLint(M.vertices.nb()));
249		✗	nlEigenSolverParameteri(NL_NB_EIGENS, (NLint)nb_eigens);
250		✗	nlEigenSolverParameterd(NL_EIGEN_SHIFT, shift);
251
252		✗	if(discretization == COMBINATORIAL) {
253		✗	nlEigenSolverParameteri(NL_SYMMETRIC, NL_TRUE);
254			}
255
256		✗	nlEnable(NL_VARIABLES_BUFFER);
257
258		✗	nlBegin(NL_SYSTEM);
259
260		✗	Attribute<double> eigen_vector;
261		✗	eigen_vector.create_vector_attribute(
262			M.vertices.attributes(), attribute_name, nb_eigens
263			);
264
265		✗	for(index_t eigen=0; eigen<nb_eigens; ++eigen) {
266			// Bind directly the variables buffer to the attribute in
267			// the mesh, to avoid copying data.
268		✗	nlBindBuffer(
269			NL_VARIABLES_BUFFER,
270			NLuint(eigen),
271		✗	&eigen_vector[0] + eigen, // base address for eigenvector
272			NLuint(sizeof(double)*nb_eigens) // number of bytes between two
273			// consecutive components in current eigenvector
274			);
275			}
276
277			// Step 2: assemble matrices
278			// *************************
279
280		✗	assemble_Laplacian_matrices(M, discretization);
281
282		✗	nlEnd(NL_SYSTEM);
283
284			// Step 3: solve and cleanup
285			// *************************
286
287		✗	nlEigenSolve();
288
289		✗	if(print_spectrum) {
290		✗	for(index_t i=0; i<nb_eigens; ++i) {
291		✗	Logger::out("MH") << i << ":" << nlGetEigenValue(i)
292		✗	<< std::endl;
293			}
294			}
295
296		✗	nlDeleteContext(nlGetCurrent());
297		✗	}
298
299
300		✗	void mesh_compute_manifold_harmonics_by_bands(
301			Mesh& M, index_t nb_eigens,
302			LaplaceBeltramiDiscretization discretization,
303			ManifoldHarmonicsCallback callback,
304			index_t nb_eigens_per_band,
305			double initial_shift,
306			void* client_data
307			) {
308
309			// Step 1: configure eigen solver and assemble matrices
310			// ****************************************************
311
312		✗	if(!nlInitExtension("ARPACK")) {
313		✗	Logger::err("MH")
314		✗	<< "Could not initialize OpenNL ARPACK extension"
315		✗	<< std::endl;
316		✗	return;
317			}
318
319		✗	nlNewContext();
320
321		✗	nlEigenSolverParameteri(NL_EIGEN_SOLVER, NL_ARPACK_EXT);
322		✗	nlEigenSolverParameteri(NL_NB_VARIABLES, NLint(M.vertices.nb()));
323		✗	nlEigenSolverParameteri(NL_NB_EIGENS, (NLint)nb_eigens_per_band);
324
325		✗	if(discretization == COMBINATORIAL) {
326		✗	nlEigenSolverParameteri(NL_SYMMETRIC, NL_TRUE);
327			}
328
329		✗	nlBegin(NL_SYSTEM);
330		✗	assemble_Laplacian_matrices(M, discretization);
331		✗	nlEnd(NL_SYSTEM);
332
333
334			// Step 2: main loop
335			// *****************
336
337		✗	double shift = initial_shift;
338		✗	index_t current_eigen = 0;
339		✗	index_t current_band = 0;
340		✗	double latest_eigen = 0.0;
341
342		✗	vector<double> eigen_vector(M.vertices.nb());
343
344			for(;;) {
345
346		✗	bool compute_band = true;
347		✗	while(compute_band) {
348		✗	Logger::out("MH")
349		✗	<< "Compute band, shift=" << shift << std::endl;
350		✗	nlEigenSolverParameterd(NL_EIGEN_SHIFT, shift);
351		✗	nlEigenSolve();
352		✗	compute_band = false;
353
354			// Test whether the current band overlaps the previous one.
355			// If this is not the case, go back (move shift towards zero)
356			// a little bit.
357
358		✗	if(current_band != 0) {
359		✗	if(::fabs(nlGetEigenValue(0)) > ::fabs(latest_eigen)) {
360		✗	Logger::out("MH")
361		✗	<< "Bands do no overlap (going back a little bit)"
362		✗	<< std::endl;
363		✗	shift -= 0.2*(
364		✗	nlGetEigenValue(nb_eigens_per_band - 1) -
365		✗	nlGetEigenValue(0)
366			);
367		✗	compute_band = true;
368			}
369			}
370
371			}
372
373		✗	for(index_t i=0; i<nb_eigens_per_band; ++i) {
374			// Output all the eigenpairs with an eigenvalue that was
375			// not previously seen (ignore the part of the current band
376			// that overlaps the previous band).
377		✗	if(
378		✗	current_eigen == 0 ||
379		✗	::fabs(nlGetEigenValue(i)) > ::fabs(latest_eigen)
380			) {
381		✗	latest_eigen = nlGetEigenValue(i);
382		✗	for(index_t j: M.vertices) {
383		✗	eigen_vector[j] = nlMultiGetVariable(j,i);
384			}
385		✗	callback(
386			current_eigen,
387		✗	nlGetEigenValue(i), eigen_vector.data(), client_data
388			);
389		✗	++current_eigen;
390		✗	if(current_eigen >= nb_eigens) {
391		✗	nlDeleteContext(nlGetCurrent());
392		✗	return;
393			}
394			}
395			}
396
397			// Move to next band / next eigen shift.
398		✗	++current_band;
399		✗	shift += 0.8 * (
400		✗	nlGetEigenValue(nb_eigens_per_band - 1) - nlGetEigenValue(0)
401			);
402		✗	}
403		✗	}
404
405
406			}
407