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/parameterization/mesh_PGP_2d.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/NL/nl.h>
44			#include <geogram/bibliography/bibliography.h>
45
46			#include <stack>
47
48			namespace {
49			using namespace GEO;
50
51
52			/**
53			* \brief computes gradients in a triangle in 2D.
54			*/
55			class ParamTrglGradient {
56			public:
57		✗	ParamTrglGradient(
58			const vec2& p1, const vec2& p2, const vec2& p3
59		✗	) {
60		✗	vertex_[0] = p1 ;
61		✗	vertex_[1] = p2 ;
62		✗	vertex_[2] = p3 ;
63
64		✗	double x1 = p1.x ;
65		✗	double y1 = p1.y ;
66		✗	double x2 = p2.x ;
67		✗	double y2 = p2.y ;
68		✗	double x3 = p3.x ;
69		✗	double y3 = p3.y ;
70
71		✗	double d = x2*y3 - y2*x3 + x3*y1 - y3*x1 + x1*y2 - y1*x2 ;
72
73		✗	if(fabs(d) < 1e-10) {
74			d = 1.0 ;
75		✗	is_flat_ = true ;
76			} else {
77		✗	is_flat_ = false ;
78			}
79
80		✗	TX_[0] = (y2 - y3)/d ;
81		✗	TX_[1] = (y3 - y1)/d ;
82		✗	TX_[2] = (y1 - y2)/d ;
83
84		✗	TY_[0] = -(x2 - x3)/d ;
85		✗	TY_[1] = -(x3 - x1)/d ;
86		✗	TY_[2] = -(x1 - x2)/d ;
87		✗	}
88
89			double TX(int i) const {
90			geo_debug_assert(i<3);
91		✗	return TX_[i];
92			}
93
94			double TY(int i) const {
95			geo_debug_assert(i<3);
96		✗	return TY_[i];
97			}
98
99			bool is_flat() const {
100		✗	return is_flat_ ;
101			}
102
103			private:
104			double TX_[3] ;
105			double TY_[3] ;
106			vec2 vertex_[3] ;
107			bool is_flat_ ;
108			} ;
109
110
111			/**
112			* \brief Retrieves a coordinate from an angle computed by PGP.
113			* \param[in] alpha the input variable.
114			* \param[in] ref the reference variable.
115			* \return a number congruent to \p alpha modulo 2 pi in the inverval
116			* [ \p ref - M_PI, \p ref + M_PI]
117			*/
118		✗	double normalize_periodic_variable(
119			double alpha, double ref
120			) {
121			int count = 0;
122		✗	if(Numeric::is_nan(alpha)) {
123			return 0.0 ;
124			}
125		✗	if(Numeric::is_nan(ref)) {
126			return 0.0 ;
127			}
128			double result = alpha ;
129			count = 0 ;
130		✗	while(ref - result > M_PI) {
131		✗	result += 2.0 * M_PI ;
132		✗	count ++ ;
133		✗	if(count > 100) {
134			return 0.0 ;
135			}
136			}
137			count = 0 ;
138		✗	while(result - ref > M_PI) {
139		✗	result -= 2.0 * M_PI ;
140		✗	count ++ ;
141		✗	if(count > 100) {
142			return 0.0 ;
143			}
144			}
145			return result ;
146			}
147			}
148
149			namespace GEO {
150
151			namespace GlobalParam2d {
152
153		✗	void PGP(
154			Mesh* mesh,
155			Attribute<vec3>& B, Attribute<vec2>& U,
156			double scaling, bool constrain_hard_edges,
157			bool use_direct_solver,
158			double maximum_scaling_correction
159			) {
160		✗	geo_cite("DBLP:journals/tog/RayLLSA06");
161
162			bool do_brush = true;
163
164			// We use f = c/3
165			// (could be fixed in the future, using a c2f array...).
166		✗	geo_assert(mesh->facets.are_simplices());
167
168			// Step 0: Preparation
169
170		✗	scaling *= 2.0;
171		✗	scaling *= surface_average_edge_length(*mesh);
172
173			// Step 0.1: Preparation / Brushing
174			if(do_brush) {
175		✗	Internal::brush(mesh,B);
176			}
177
178			// Step 0.2: Preparation / Compute relative rotation of B between
179			// pairs of adjacent facets
180		✗	Attribute<index_t> R_ff(mesh->facet_corners.attributes(),"R");
181		✗	Internal::compute_R_ff(mesh,B,R_ff);
182		✗	Attribute<index_t> R_fv(mesh->facet_corners.attributes(),"R_fv");
183		✗	Internal::compute_R_fv(mesh,R_ff,R_fv);
184
185			Attribute<double> CC;
186		✗	if(maximum_scaling_correction != 1.0) {
187		✗	Logger::out("PGP") << "Computing scaling correction"
188			<< std::endl;
189		✗	CC.bind(mesh->vertices.attributes(), "CC");
190		✗	Attribute<vec3> Bv(mesh->vertices.attributes(), "B");
191		✗	Internal::transfer_B_to_vertices(mesh, B, Bv, R_fv);
192		✗	curl_correction(
193			mesh, Bv, R_fv, CC,
194			use_direct_solver, maximum_scaling_correction
195			);
196			}
197
198
199		✗	Logger::out("PGP") << "Solving for PGP" << std::endl;
200
201
202			// Step 1: Determine the structure of the problem:
203			// - There are four variables per vertex (cu, su, cv, sv)
204			// - Each vertex has a reference corner (v2c[v])
205			// - Each corner c knows the number of rotations Rc[c] required
206			// to make the B of its triangle match the B of the triangle
207			// of the reference corner attached to its vertex (clear
208			// enough ?)
209
210			// Rot[ R_ff[c] ][i][j] corresponds to the transform to
211			// be applied to the variables associated with a
212			// vertex to express their coordinates in the frame
213			// of corner c (it corresponds to the same Rot[][][]
214			// matrices as QuadCover, with the exception that each
215			// coefficient 1 is replaced with the 2x2 identity matrix
216			// and each coefficient -1 with Mat2x2(1,0,0,-1) (there is
217			// only one "-1" coefficient because cos(-x) = cos(x) !!).
218			// Note: in the PGP paper, the wrong coefficient is negated.
219
220			static const double Rot[4][4][4] = {
221			{{ 1, 0, 0, 0},
222			{ 0, 1, 0, 0},
223			{ 0, 0, 1, 0},
224			{ 0, 0, 0, 1}},
225
226			{{ 0, 0, 1, 0},
227			{ 0, 0, 0, 1},
228			{ 1, 0, 0, 0},
229			{ 0,-1, 0, 0}},
230
231			{{ 1, 0, 0, 0},
232			{ 0,-1, 0, 0},
233			{ 0, 0, 1, 0},
234			{ 0, 0, 0,-1}},
235
236			{{ 0, 0, 1, 0},
237			{ 0, 0, 0,-1},
238			{ 1, 0, 0, 0},
239			{ 0, 1, 0, 0}}
240			};
241
242			static double Rotc2[4][4];
243
244			// Step 2: setup and solve linear system
245			{
246
247		✗	nlNewContext();
248
249		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
250		✗	nlSolverParameteri(
251			NL_NB_VARIABLES, NLint(mesh->vertices.nb()*4)
252			);
253
254		✗	if(use_direct_solver) {
255		✗	if(nlInitExtension("CHOLMOD")) {
256		✗	nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
257		✗	} else if(nlInitExtension("SUPERLU")) {
258		✗	nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT);
259			} else {
260		✗	Logger::warn("PGP")
261			<< "Could not initialize direct sovlver"
262			<< std::endl;
263		✗	Logger::warn("PGP")
264			<< "Falling back to Jacobi pre-CG"
265			<< std::endl;
266			use_direct_solver = false;
267			}
268			}
269
270			if(!use_direct_solver) {
271			// With the iterative solver, a very small threshold is
272			// needed, because of the very high scaling on the
273			// solution due to the varying modulus of the complex
274			// numbers far away from the constrained point.
275		✗	nlSolverParameterd(NL_THRESHOLD, 1e-20);
276			}
277
278		✗	nlBegin(NL_SYSTEM);
279
280		✗	if(constrain_hard_edges) {
281		✗	for(index_t c: mesh->facet_corners) {
282			index_t cnstr =
283		✗	Internal::get_edge_constraints(mesh, c, B);
284			index_t v = mesh->facet_corners.vertex(c);
285			// Inverse R, because when the axis turns
286			// clockwise, coordinates turn anticlockwise.
287		✗	index_t Rcc = Internal::inverse_R(R_fv[c]);
288		✗	if(cnstr & Internal::CNSTR_U) {
289		✗	if(Rcc == 0 || Rcc == 2) {
290		✗	nlLockVariable(4*v);
291			nlSetVariable(4*v,1e4);
292		✗	nlLockVariable(4*v+1);
293			nlSetVariable(4*v+1,0.0);
294			} else {
295		✗	nlLockVariable(4*v+2);
296			nlSetVariable(4*v+2,1e4);
297		✗	nlLockVariable(4*v+3);
298			nlSetVariable(4*v+3,0.0);
299			}
300			}
301		✗	if(cnstr & Internal::CNSTR_V) {
302		✗	if(Rcc == 0 || Rcc == 2) {
303		✗	nlLockVariable(4*v+2);
304			nlSetVariable(4*v+2,1e4);
305		✗	nlLockVariable(4*v+3);
306			nlSetVariable(4*v+3,0.0);
307			} else {
308		✗	nlLockVariable(4*v);
309			nlSetVariable(4*v,1e4);
310		✗	nlLockVariable(4*v+1);
311			nlSetVariable(4*v+1,0.0);
312			}
313			}
314			}
315			}
316
317			// Lock at least one variable per connected component.
318			{
319		✗	std::vector<bool> f_visited(mesh->facets.nb(),false);
320		✗	for(index_t f: mesh->facets) {
321		✗	if(!f_visited[f]) {
322			index_t nb_locked=0;
323			index_t first_v = mesh->facets.vertex(f,0);
324			std::stack<index_t> S;
325			f_visited[f] = true;
326			S.push(f);
327		✗	while(!S.empty()) {
328		✗	index_t f_top = S.top();
329			S.pop();
330		✗	FOR(le,mesh->facets.nb_vertices(f_top)) {
331			index_t v = mesh->facets.vertex(f_top,le);
332			if(
333		✗	nlVariableIsLocked(4*v) ||
334		✗	nlVariableIsLocked(4*v+2)) {
335		✗	++nb_locked;
336			}
337			index_t f_neigh =
338		✗	mesh->facets.adjacent(f_top,le);
339
340		✗	if(f_neigh != NO_FACET &&
341			!f_visited[f_neigh]
342			) {
343			f_visited[f_neigh] = true;
344			S.push(f_neigh);
345			}
346			}
347			}
348
349			// Lock one of the points in each
350			// connected component to make sure that
351			// the minimum is well defined.
352		✗	if(nb_locked == 0) {
353		✗	nlLockVariable(4*first_v);
354			nlSetVariable(4*first_v,1e4);
355
356		✗	nlLockVariable(4*first_v+1);
357			nlSetVariable(4*first_v+1,0.0);
358
359		✗	nlLockVariable(4*first_v+2);
360			nlSetVariable(4*first_v+2,1e4);
361
362		✗	nlLockVariable(4*first_v+3);
363			nlSetVariable(4*first_v+3,0.0);
364			}
365			}
366			}
367			}
368
369		✗	nlBegin(NL_MATRIX);
370
371			// This one will be replaced in-place
372			// with the rotation that encodes the
373			// delta u and delta v along each edge.
374
375			double RotDelta[4][4];
376		✗	FOR(i,4) {
377		✗	FOR(j,4) {
378		✗	RotDelta[i][j] = ((i==j) ? 1.0 : 0.0);
379			}
380			}
381
382		✗	for(index_t f: mesh->facets) {
383
384			vec3 Bf, BTf;
385
386			for(index_t c1=mesh->facets.corners_begin(f);
387		✗	c1 < mesh->facets.corners_end(f); ++c1) {
388
389		✗	Internal::get_B_on_edge(mesh, B, R_ff, f, c1, Bf, BTf);
390
391			index_t c2 =
392			mesh->facets.next_corner_around_facet(f,c1);
393			index_t v1 = mesh->facet_corners.vertex(c1);
394			index_t v2 = mesh->facet_corners.vertex(c2);
395
396			vec3 E =
397			mesh->vertices.point(v2) -
398			mesh->vertices.point(v1) ;
399		✗	double delta_u = 2.0 * M_PI * dot(E,Bf)/scaling;
400		✗	double delta_v = 2.0 * M_PI * dot(E,BTf)/scaling;
401
402			if(CC.is_bound()) {
403		✗	double s = (0.5*(CC[v1] + CC[v2]));
404		✗	delta_u /= s;
405		✗	delta_v /= s;
406			}
407
408		✗	double sdu = sin(delta_u);
409		✗	double cdu = cos(delta_u);
410		✗	double sdv = sin(delta_v);
411		✗	double cdv = cos(delta_v);
412
413		✗	RotDelta[0][0] = cdu;
414		✗	RotDelta[0][1] = -sdu;
415		✗	RotDelta[1][0] = sdu;
416		✗	RotDelta[1][1] = cdu;
417
418		✗	RotDelta[2][2] = cdv;
419		✗	RotDelta[2][3] = -sdv;
420		✗	RotDelta[3][2] = sdv;
421		✗	RotDelta[3][3] = cdv;
422
423			// Inverse R, because when the axis turns
424			// clockwise, coordinates turn anticlockwise.
425		✗	index_t Rc1 = Internal::inverse_R(R_fv[c1]);
426		✗	index_t Rc2 = Internal::inverse_R(R_fv[c2]);
427
428			// Compute the product of the "delta u, delta v"
429			// rotation with the Rc2 "90 degrees rotation" matrix,
430			// exactly like in the PGP article.
431
432		✗	for(index_t i=0; i<4; ++i) {
433		✗	for(index_t j=0; j<4; ++j) {
434		✗	Rotc2[i][j] = 0.0;
435		✗	for(index_t k=0; k<4; ++k) {
436		✗	Rotc2[i][j] +=
437		✗	RotDelta[i][k] * Rot[Rc2][k][j];
438			}
439			}
440			}
441
442		✗	for(index_t i=0; i<4; ++i) {
443		✗	nlBegin(NL_ROW);
444		✗	for(index_t j=0; j<4; ++j) {
445		✗	double a1 = Rot[Rc1][i][j];
446		✗	if(a1 != 0.0) {
447		✗	nlCoefficient(v1*4+j, a1);
448			}
449			}
450		✗	for(index_t j=0; j<4; ++j) {
451		✗	double a2 = Rotc2[i][j];
452		✗	if(a2 != 0.0) {
453		✗	nlCoefficient(v2*4+j, -a2);
454			}
455			}
456		✗	nlEnd(NL_ROW);
457			}
458			}
459			}
460
461		✗	nlEnd(NL_MATRIX);
462		✗	nlEnd(NL_SYSTEM);
463
464		✗	nlSolve();
465
466			Attribute<double> PGP;
467		✗	PGP.bind_if_is_defined(mesh->facet_corners.attributes(),"PGP");
468			if(!PGP.is_bound()) {
469		✗	PGP.create_vector_attribute(
470		✗	mesh->facet_corners.attributes(), "PGP", 4
471			);
472			}
473
474		✗	for(index_t c: mesh->facet_corners) {
475			index_t v = mesh->facet_corners.vertex(c);
476			// Inverse R, because when the axis turns
477			// clockwise, coordinates turn anticlockwise.
478		✗	index_t Rcc = Internal::inverse_R(R_fv[c]);
479			double vars[4];
480		✗	for(index_t i=0; i<4; ++i) {
481		✗	vars[i] = 0.0;
482		✗	for(index_t j=0; j<4; ++j) {
483		✗	vars[i] += Rot[Rcc][i][j] * nlGetVariable(4*v+j);
484			}
485			}
486		✗	double s = sqrt(vars[0]*vars[0]+vars[1]*vars[1]);
487		✗	vars[0] /= s;
488		✗	vars[1] /= s;
489
490		✗	s = sqrt(vars[2]*vars[2]+vars[3]*vars[3]);
491		✗	vars[2] /= s;
492		✗	vars[3] /= s;
493
494		✗	U[c].x = atan2(vars[1], vars[0]);
495		✗	U[c].y = atan2(vars[3], vars[2]);
496
497		✗	FOR(i,4) {
498		✗	PGP[4*c+i] = vars[i];
499			}
500			}
501
502			Attribute<index_t> singular(
503		✗	mesh->facets.attributes(),"is_singular"
504			);
505
506		✗	for(index_t f: mesh->facets) {
507		✗	singular[f] = false;
508			for(index_t c1=mesh->facets.corners_begin(f);
509		✗	c1<mesh->facets.corners_end(f); ++c1) {
510			index_t c2 =
511			mesh->facets.next_corner_around_facet(f,c1);
512
513			index_t v1 = mesh->facet_corners.vertex(c1);
514			index_t v2 = mesh->facet_corners.vertex(c2);
515			vec3 E =
516			mesh->vertices.point(v2) -
517			mesh->vertices.point(v1);
518
519			vec3 Bf, BTf;
520		✗	Internal::get_B_on_edge(mesh, B, R_ff, f, c1, Bf, BTf);
521
522		✗	double delta_u = 2.0 * M_PI * dot(E,Bf)/scaling;
523		✗	double delta_v = 2.0 * M_PI * dot(E,BTf)/scaling;
524
525			if(CC.is_bound()) {
526		✗	double s = (0.5*(CC[v1] + CC[v2]));
527		✗	delta_u /= s;
528		✗	delta_v /= s;
529			}
530
531		✗	double expected_u = U[c1].x-delta_u;
532		✗	double expected_v = U[c1].y-delta_v;
533			vec2 uv(
534		✗	normalize_periodic_variable(U[c2].x,expected_u),
535		✗	normalize_periodic_variable(U[c2].y,expected_v)
536			);
537		✗	if(c2 == mesh->facets.corners_begin(f)) {
538		✗	singular[f] =
539		✗	(length(uv - U[c2]) > 1e-20) ? 1 : 0;
540			} else {
541		✗	U[c2] = uv;
542			}
543
544			}
545			}
546
547		✗	for(index_t c: mesh->facet_corners) {
548		✗	U[c].x /= M_PI;
549		✗	U[c].y /= M_PI;
550		✗	Internal::snap_tex_coord(U[c].x);
551		✗	Internal::snap_tex_coord(U[c].y);
552			}
553
554		✗	nlDeleteContext(nlGetCurrent());
555			}
556
557
558		✗	}
559
560
561		✗	void curl_correction(
562			Mesh* mesh, Attribute<vec3>& Bv,
563			Attribute<index_t>& R_fv, Attribute<double>& CC,
564			bool use_direct_solver, double max_scaling_correction
565			) {
566		✗	geo_assert(Bv.manager() == &mesh->vertices.attributes());
567		✗	nlNewContext();
568
569		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
570		✗	nlSolverParameteri(
571			NL_NB_VARIABLES, NLint(mesh->vertices.nb()*4)
572			);
573
574		✗	if(use_direct_solver) {
575		✗	if(nlInitExtension("CHOLMOD")) {
576		✗	nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
577		✗	} else if(nlInitExtension("SUPERLU")) {
578		✗	nlSolverParameteri(NL_SOLVER, NL_PERM_SUPERLU_EXT);
579			} else {
580		✗	Logger::warn("PGP")
581			<< "Could not initialize direct solver"
582			<< std::endl;
583		✗	Logger::warn("PGP")
584			<< "Falling back to Jacobi pre-CG"
585			<< std::endl;
586			use_direct_solver = false;
587			}
588			}
589
590			if(!use_direct_solver) {
591			// With the iterative solver, a very small threshold is
592			// needed, because of the very high scaling on the
593			// solution due to the varying modulus of the complex
594			// numbers far away from the constrained point.
595		✗	nlSolverParameterd(NL_THRESHOLD, 1e-20);
596			}
597
598			const double locked_value = 1.0;
599			const double solver_scale = 1e3;
600
601		✗	nlBegin(NL_SYSTEM);
602		✗	nlLockVariable(0u);
603		✗	nlSetVariable(0u,locked_value);
604		✗	nlBegin(NL_MATRIX);
605		✗	for(index_t f: mesh->facets) {
606
607			index_t va = mesh->facets.vertex(f,0);
608			index_t vb = mesh->facets.vertex(f,1);
609			index_t vc = mesh->facets.vertex(f,2);
610
611		✗	vec3 N = normalize(Geom::mesh_facet_normal(*mesh,f));
612
613		✗	vec3 fieldA3d = normalize(Bv[va]);
614		✗	vec3 fieldB3d = normalize(Bv[vb]);
615		✗	vec3 fieldC3d = normalize(Bv[vc]);
616
617		✗	FOR(i, Internal::inverse_R(R_fv[f*3])) {
618			fieldA3d = cross(N,fieldA3d);
619			}
620
621		✗	FOR(i, Internal::inverse_R(R_fv[f*3+1])) {
622			fieldB3d = cross(N,fieldB3d);
623			}
624
625		✗	FOR(i, Internal::inverse_R(R_fv[f*3+2])) {
626			fieldC3d = cross(N,fieldC3d);
627			}
628
629		✗	vec3 U = normalize(cross(cross(N, fieldA3d),N));
630		✗	vec3 V = normalize(cross(N, U));
631
632
633			// fieldx is the 2d field at point x. The fields are
634			// rotated in a coherent way to take the modulus into
635			// account notice that it will be better if the field
636			// is rotated to the facet instead of being projected
637			vec2 fieldA(1.0,0.0);
638			vec2 fieldB(dot(fieldB3d, U), dot(fieldB3d, V)) ;
639		✗	fieldB = normalize(fieldB);
640
641			vec2 fieldC (dot(fieldC3d, U), dot(fieldC3d, V)) ;
642		✗	fieldC = normalize(fieldC);
643
644		✗	vec3 A3d = mesh->vertices.point(va);
645		✗	vec3 B3d = mesh->vertices.point(vb);
646		✗	vec3 C3d = mesh->vertices.point(vc);
647
648			vec3 AB3d = B3d-A3d;
649			vec3 AC3d = C3d-A3d;
650
651			vec2 A(0.0,0.0);
652			vec2 B(dot(AB3d,U), dot(AB3d,V));
653			vec2 C(dot(AC3d,U), dot(AC3d,V));
654
655			// hummm... flat triangles ?
656			double a=0;
657			double b=0;
658		✗	ParamTrglGradient trg(A,B,C);
659		✗	if (trg.is_flat() ||
660		✗	length(cross(AB3d,AC3d))<1e-10 ||
661		✗	::fabs(fieldB.x) < 1e-20 ||
662		✗	::fabs(fieldC.x) < 1e-20
663			) {
664			std::cerr<<"bad triangle ..." << std::endl ;
665
666		✗	nlBegin(NL_ROW);
667			nlCoefficient(va,1e-2) ;
668			nlCoefficient(vb,-1e-2) ;
669		✗	nlEnd(NL_ROW);
670
671		✗	nlBegin(NL_ROW);
672			nlCoefficient(va,1e-2) ;
673			nlCoefficient(vc,-1e-2) ;
674		✗	nlEnd(NL_ROW);
675
676		✗	nlBegin(NL_ROW);
677			nlCoefficient(vc,1e-2) ;
678			nlCoefficient(vb,-1e-2) ;
679		✗	nlEnd(NL_ROW);
680		✗	continue ;
681			} else {
682
683			// the direction field is represented by angles
684			// double alpha_A = 0;
685			// double alpha_B = atan(fieldB.y/fieldB.x);
686			// double alpha_C = atan(fieldC.y/fieldC.x);
687
688			// Order 1 Taylor Expansion ....
689			double alpha_A = 0;
690		✗	double alpha_B = fieldB.y/fieldB.x;
691		✗	double alpha_C = fieldC.y/fieldC.x;
692
693			vec2 grad_alpha (
694		✗	trg.TX(0)*alpha_A +
695		✗	trg.TX(1)*alpha_B +
696		✗	trg.TX(2)*alpha_C ,
697		✗	trg.TY(0)*alpha_A +
698		✗	trg.TY(1)*alpha_B +
699		✗	trg.TY(2)*alpha_C
700			);
701			a = grad_alpha.x;
702			b = grad_alpha.y;
703			}
704
705		✗	double sqrt_area = ::sqrt(Geom::mesh_facet_area(*mesh,f)) ;
706
707			if(
708		✗	Numeric::is_nan(sqrt_area) ||
709		✗	Numeric::is_nan(trg.TX(0)) ||
710		✗	Numeric::is_nan(trg.TX(1)) ||
711		✗	Numeric::is_nan(trg.TX(2)) ||
712		✗	Numeric::is_nan(trg.TY(0)) ||
713		✗	Numeric::is_nan(trg.TY(1)) ||
714		✗	Numeric::is_nan(trg.TY(2)) ||
715		✗	Numeric::is_nan(a) ||
716		✗	Numeric::is_nan(b)
717			) {
718			std::cerr << "Found NAN !!" << std::endl ;
719			}
720
721
722		✗	nlRowScaling(sqrt_area) ;
723		✗	nlBegin(NL_ROW);
724			nlCoefficient(va, trg.TX(0)) ;
725			nlCoefficient(vb, trg.TX(1)) ;
726			nlCoefficient(vc, trg.TX(2)) ;
727		✗	nlRightHandSide(solver_scale*b);
728		✗	nlEnd(NL_ROW);
729
730		✗	nlRowScaling(sqrt_area) ;
731		✗	nlBegin(NL_ROW);
732			nlCoefficient(va, trg.TY(0)) ;
733			nlCoefficient(vb, trg.TY(1)) ;
734			nlCoefficient(vc, trg.TY(2)) ;
735		✗	nlRightHandSide(solver_scale*-a);
736		✗	nlEnd(NL_ROW);
737			}
738		✗	nlEnd(NL_MATRIX);
739		✗	nlEnd(NL_SYSTEM);
740
741		✗	nlSolve();
742		✗	double min_val = Numeric::max_float64();
743		✗	double max_val = Numeric::min_float64();
744		✗	for(index_t v: mesh->vertices) {
745		✗	CC[v] = ::exp(
746		✗	(nlGetVariable(v)-locked_value)/solver_scale)
747			;
748		✗	min_val = std::min(min_val,CC[v]);
749		✗	max_val = std::max(max_val,CC[v]);
750			}
751		✗	double avg_val = 0.5 * (min_val + max_val);
752		✗	double min_limit = 1.0 / ::sqrt(max_scaling_correction);
753		✗	double max_limit = ::sqrt(max_scaling_correction);
754		✗	for(index_t v: mesh->vertices) {
755		✗	CC[v] = CC[v] / avg_val;
756			geo_clamp(CC[v], min_limit, max_limit);
757			}
758		✗	nlDeleteContext(nlGetCurrent());
759		✗	}
760
761			} // namespace GlobalParam2d
762			} // namespace OGF
763