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
41			#ifndef H_HEXDOM_ALGO_MIXED_CONSTRAINED_SOLVER_H
42			#define H_HEXDOM_ALGO_MIXED_CONSTRAINED_SOLVER_H
43
44
45			#include <exploragram/basic/common.h>
46			#include <exploragram/hexdom/basic.h>
47			#include <exploragram/hexdom/id_map.h>
48			#include <geogram/NL/nl.h>
49			#include <algorithm>
50
51
52
53			//***************************************************
54			// HOW TO USE:
55			//***************************************************
56			// MatrixMixedConstrainedSolver intsolver(nb_vars);
57			// // tell that variable var_id is a multiple of m
58			// intsolver.set_multiplicity(var_id, m);
59			// // tell 4 times the linear constraints
60			// FOR(pass, 4) {
61			// // construct each constraint
62			// intsolver.begin_constraint();
63			// intsolver.add_constraint_coeff(var_id, coeff);
64			// ....
65			// intsolver.end_constraint();
66			// // tell the solver to start the next pass
67			// intsolver.end_pass(pass);
68			//}
69			// // solve in two pass
70			// FOR(iter, 2) {
71			// if (iter == 0) intsolver.start_full_real_iter();
72			// else intsolver.start_mixed_iter();
73			// intsolver.begin_energy();
74			// intsolver.add_energy_coeff(var_id, coeff);
75			// ....
76			// intsolver.add_energy_rhs(rhs);
77			// intsolver.end_energy();
78			// }
79			// intsolver.end_iter();
80			//}
81			// // get the result
82			//double var_value = intsolver.value(var_id);
83
84
85			namespace GEO {
86
87			struct Coeff {
88		✗	Coeff(index_t p_index, double p_a) : index(p_index), a(p_a) {
89			}
90
91		✗	Coeff() : index(index_t(-1)), a(-1.0) {
92			}
93
94			index_t index;
95			double a;
96			};
97
98			inline std::ostream& operator<<(
99			std::ostream& os, const GEO::vector<Coeff>& obj
100			) {
101			// write obj to stream
102			FOR(i, obj.size())
103			os <<(i?" + ":"")<< obj[i].a << ".X_"<<obj[i].index;
104			return os;
105			}
106
107			//inline double nint(double x) {
108			// // return floor(x + .5);
109			// double r1 = floor(x);
110			// double r2 = floor(x + 1.0);
111			// return (x - r1) < (r2 - x) ? r1 : r2;
112			//}
113
114		✗	inline bool coeffCmp(const Coeff& A, const Coeff& B) {
115		✗	return (A.index < B.index);
116			}
117
118		✗	inline void order_AND_make_unique(vector<Coeff>& c) {
119		✗	std::sort(c.begin(), c.end(), coeffCmp);
120			double sum = 0;
121			index_t ires = 0;
122		✗	FOR(i, c.size()) {
123		✗	sum += c[i].a;
124		✗	if (i + 1 != c.size() && c[i].index == c[i + 1].index) {
125		✗	continue;
126			}
127		✗	if (std::fabs(sum)>1e-10) {
128		✗	c[ires] = Coeff(c[i].index, sum);
129		✗	ires++;
130			}
131			sum = 0;
132			}
133		✗	c.resize(ires);
134		✗	}
135
136			/**
137			* \brief Adds a coefficient in a sparse vector.
138			* \details If the coefficient already exists, then \p val
139			* is added to the existing coefficient, else a new coefficient
140			* is created and inserted in such a way that the vector remains
141			* sorted.
142			* \param[in,out] L a reference to the sparse vector,
143			* represented by a vector<Coeff>.
144			* \param[in] index the index of the coefficient to add.
145			* \param[in] val the value of the coefficient to add.
146			*/
147		✗	inline void add_to_sorted_sparse_vector(
148			vector<Coeff>& L, index_t index, double val
149			) {
150			// [BL] what follows is an optimized version of:
151			// L.push_back(Coeff(index,val));
152			// order_AND_make_unique(L);
153		✗	for (index_t jj = 0; jj < L.size(); ++jj) {
154		✗	if (L[jj].index == index) {
155		✗	L[jj].a += val;
156		✗	return;
157			}
158		✗	if (L[jj].index > index) {
159		✗	L.insert(L.begin() + long(jj), Coeff(index, val));
160		✗	return;
161			}
162			}
163		✗	L.push_back(Coeff(index, val));
164			}
165
166		✗	struct NaiveSparseMatrix {
167
168		✗	void init_zero(index_t p_nb_lines) {
169		✗	M.resize(p_nb_lines);
170			// [BL] the following line was missing (resize does
171			// not reset the existing items in a vector).
172		✗	for (index_t i = 0; i < p_nb_lines; ++i) {
173		✗	M[i].resize(0);
174			}
175		✗	}
176
177		✗	void init_identity(index_t size) {
178		✗	init_zero(size);
179		✗	FOR(i, size) {
180		✗	M[i].push_back(Coeff(i, 1));
181			}
182		✗	}
183
184			index_t nb_coeffs_in_line(index_t l) const {
185			return M[l].size();
186			}
187
188			const Coeff& get(index_t line, index_t i_th_non_null_coeff) const {
189			return M[line][i_th_non_null_coeff];
190			}
191
192			void add(index_t line, index_t col, double coeff) {
193		✗	add_to_sorted_sparse_vector(M[line], col, coeff);
194			// [BL] the code above replaces:
195			// L.push_back(Coeff(col,coeff));
196			// order_AND_make_unique(L);
197		✗	}
198
199			index_t nb_lines() const {
200			return index_t(M.size());
201			}
202
203			vector<vector<Coeff> > M;
204			};
205
206
207
208			// a RemoveColMatrix B is a matrix that zeros a column of M by M <-- MB
209		✗	struct RemoveColMatrix {
210
211		✗	RemoveColMatrix(vector<Coeff> C, index_t ind_to_remove = 0) {
212		✗	geo_assert(ind_to_remove < C.size());
213			/* if (ind_to_remove >= 0) [BL] TODO: Check with Nico !!!*/
214			std::swap(C.back(), C[ind_to_remove]);
215		✗	zero_col = C.back().index;
216		✗	double multi = -1. / C.back().a;
217			C.pop_back();
218		✗	FOR(c, C.size()) {
219		✗	line.push_back(Coeff(C[c].index, C[c].a*multi));
220			}
221		✗	}
222
223			index_t nb_coeffs_in_line(index_t l) const {
224			if (zero_col == l) {
225			return line.size();
226			}
227			return 1;
228			}
229
230			Coeff get(index_t l, int i_th_non_null_coeff) {
231			if (zero_col == l) {
232			return line[i_th_non_null_coeff];
233			}
234			return Coeff(l, 1.);
235			}
236
237			// [BL] can we remove this function ?
238			/*
239			index_t nb_lines() const {
240			geo_assert_not_reached;
241			return NOT_AN_ID;
242			}
243			*/
244
245			index_t zero_col;
246			vector<Coeff > line;
247			};
248
249
250			// [BL] passed result by reference instead of making it returned.
251		✗	inline void mult(const vector<Coeff>& c, NaiveSparseMatrix& B, vector<Coeff>& result) {
252		✗	result.resize(0);
253		✗	FOR(co_c, c.size()) {
254		✗	index_t lB = c[co_c].index;
255		✗	FOR(co_b, B.nb_coeffs_in_line(lB)) {
256		✗	result.push_back(Coeff(B.get(lB, co_b).index, B.get(lB, co_b).a * c[co_c].a));
257			}
258			}
259		✗	order_AND_make_unique(result);
260		✗	}
261
262			// MatrixM is filled in 3 passes
263			// pass 1: creates M0 that sets some variables to zero
264			// pass 2: creates M1 that creates an indirection
265			// pass 3: creates M2 that manage all remaining problems
266
267			struct MatrixM {
268
269		✗	MatrixM(index_t size) {
270		✗	multiplicity.resize(size, 0);
271		✗	DEBUG_init_multiplicity.resize(size, 0);
272		✗	M0.resize(size);
273		✗	FOR(i, size) {
274		✗	M0[i] = i;
275			}
276		✗	pass = 0;
277		✗	}
278
279		✗	void finalize_M0() {
280		✗	geo_assert(M1.size() == 0);
281		✗	geo_assert(M2.nb_lines() == 0);
282			index_t M0_nb_col = 0;
283		✗	FOR(i, M0.size()) {
284		✗	if (M0[i] != NOT_AN_ID) {
285		✗	multiplicity[M0_nb_col] = multiplicity[i];
286		✗	M0[i] = M0_nb_col;
287		✗	M0_nb_col++;
288			}
289			}
290		✗	GEO::Logger::out("HexDom") << "M0 final size is " << M0.size() << " x " << M0_nb_col << std::endl;
291		✗	multiplicity.resize(M0_nb_col);
292		✗	M1.resize(M0_nb_col);
293		✗	FOR(i, M0_nb_col) {
294		✗	M1[i] = Coeff(i, 1);
295			}
296		✗	}
297
298		✗	void finalize_M1() {
299		✗	geo_assert(M2.nb_lines() == 0);
300
301			// find M1^n (n-> infty)
302		✗	FOR(i, M1.size()) {
303		✗	while (M1[i].index != M1[M1[i].index].index) {
304		✗	M1[i].a *= M1[M1[i].index].a;
305		✗	M1[i].index = M1[M1[i].index].index;
306			}
307			}
308
309			// compress
310			index_t M1_nb_col = 0;
311			vector<index_t> to_grp(M1.size(), NOT_AN_ID);
312		✗	FOR(i, M1.size()) if (M1[i].index == i) {
313		✗	to_grp[i] = M1_nb_col;
314		✗	M1_nb_col++;
315			}
316		✗	vector<int> new_multiplicity(M1_nb_col, 0);
317		✗	FOR(i, M1.size()) {
318		✗	index_t g = to_grp[M1[i].index];
319		✗	new_multiplicity[g] = std::max(new_multiplicity[g], multiplicity[M1[i].index]);
320		✗	M1[i].index = g;
321			}
322			multiplicity = new_multiplicity;
323		✗	GEO::Logger::out("HexDom") << "M1 final size is " << M1.size() << " x " << M1_nb_col << std::endl;
324		✗	M2.init_identity(M1_nb_col);
325		✗	M2t.init_identity(M1_nb_col);
326
327		✗	}
328
329		✗	void finalize_M2() {
330			// WARNING !
331			// we do not compress here (there will be unused variables)
332			// some linear solver may not appreciate it...
333			index_t M2_nb_col = 0;
334			M2_nb_col = M2.M.size();
335		✗	kernel_size = M2_nb_col;
336		✗	multiplicity.resize(M2_nb_col);
337		✗	GEO::Logger::out("HexDom") << "M2 final size is " << M2.nb_lines() << " x " << M2_nb_col << std::endl;
338		✗	}
339
340		✗	void end_pass(index_t p_pass) {
341		✗	geo_assert(pass == p_pass);
342		✗	if (pass == 0) {
343		✗	finalize_M0();
344			}
345		✗	else if (pass == 1) {
346		✗	finalize_M1();
347			}
348		✗	else if (pass == 2) {
349		✗	finalize_M2();
350			} else {
351		✗	check_multiplicity();
352			}
353		✗	pass++;
354		✗	}
355
356
357		✗	vector<Coeff> mult_by_M0(vector<Coeff>& c) {
358			vector<Coeff> cM0;
359		✗	FOR(i, c.size()) if (M0[c[i].index] != NOT_AN_ID) cM0.push_back(Coeff(M0[c[i].index], c[i].a));
360		✗	return cM0;
361			}
362
363		✗	vector<Coeff> mult_by_M1(vector<Coeff>& cM0) {
364			vector<Coeff> cM0M1;
365		✗	FOR(i, cM0.size()) cM0M1.push_back(Coeff(M1[cM0[i].index].index, cM0[i].a*M1[cM0[i].index].a));
366		✗	order_AND_make_unique(cM0M1);
367		✗	return cM0M1;
368			}
369
370		✗	vector<Coeff> mult_by_M2(vector<Coeff>& cM0M1) {
371			vector<Coeff> cM0M1M2;
372		✗	mult(cM0M1, M2, cM0M1M2);
373			// order_AND_make_unique(cM0M1M2); [BL]: already done in mult()
374		✗	return cM0M1M2;
375			}
376
377		✗	void inplace_mult_M2_by_B(RemoveColMatrix& B) {
378			vector<Coeff> diffline = B.line;
379		✗	add_to_sorted_sparse_vector(diffline, B.zero_col, -1.0);
380			// [BL] replaced with "insert_in_sorted_vector"
381			// diffline.push_back(Coeff(B.zero_col, -1.));
382			// order_AND_make_unique(diffline);
383
384		✗	vector<Coeff> M2col = M2t.M[B.zero_col];
385		✗	FOR(nc, diffline.size()) {
386		✗	FOR(nl, M2col.size()) {
387		✗	index_t c = diffline[nc].index;
388		✗	index_t l = M2col[nl].index;
389		✗	double coeff = diffline[nc].a*M2col[nl].a;
390			M2t.add(c, l, coeff);
391			M2.add(l, c, coeff);
392			}
393			}
394		✗	}
395
396		✗	bool add_constraint(vector<Coeff>& c) {
397
398		✗	if (pass == 0) {
399		✗	if (c.size() == 1) {
400		✗	geo_assert(c[0].a != 0);
401		✗	M0[c[0].index] = NOT_AN_ID;
402			}
403		✗	return true;
404			}
405
406		✗	vector<Coeff> cM0 = mult_by_M0(c);
407
408		✗	vector<Coeff> cM0M1 = mult_by_M1(cM0);
409		✗	if (cM0M1.size() == 0) return true;
410
411		✗	if (pass == 1) {
412		✗	if (cM0M1.size() == 2) {
413		✗	FOR(i, 2) if ((std::abs(cM0M1[i].a) != 1)) return false;
414		✗	RemoveColMatrix B(cM0M1);
415		✗	Coeff co[2] = { Coeff(B.zero_col, 1.0),B.line[0] };
416			// find root
417		✗	FOR(i, 2) {
418		✗	while (co[i].index != M1[co[i].index].index) {
419		✗	co[i].a = co[i].a*M1[co[i].index].a;
420		✗	co[i].index = M1[co[i].index].index;
421			}
422			}
423			// link
424		✗	if (multiplicity[co[0].index] < multiplicity[co[1].index])
425		✗	M1[co[0].index] = Coeff(co[1].index, co[1].a*co[0].a);
426		✗	else M1[co[1].index] = Coeff(co[0].index, co[1].a*co[0].a);
427
428			}
429		✗	return true;
430			}
431
432		✗	vector<Coeff> cM0M1M2 = mult_by_M2(cM0M1);
433		✗	if (cM0M1M2.size() == 0) {
434			return true;
435			}
436
437		✗	if (pass == 2) {
438			index_t ind_to_remove = 0;
439		✗	FOR(i, cM0M1M2.size()) {
440		✗	if (cM0M1M2[i].a!=0 &&
441		✗	double(multiplicity[cM0M1M2[i].index])* std::abs(cM0M1M2[i].a) <
442		✗	double(multiplicity[cM0M1M2[ind_to_remove].index]) * std::abs(cM0M1M2[ind_to_remove].a)) {
443			ind_to_remove = i;
444			}
445			}
446		✗	RemoveColMatrix B(cM0M1M2, ind_to_remove);
447		✗	inplace_mult_M2_by_B(B);
448			}
449			// debug
450		✗	geo_assert(pass != 3 || cM0M1M2.empty());// "Kernel must be complete ";
451			return true;
452
453			}
454
455		✗	void check_multiplicity() {
456		✗	std::cerr << " check_that multiplicity are respected\n";
457		✗	FOR(x, M0.size()) if (debug_var_multiplicity(x) < DEBUG_init_multiplicity[x]) {
458		✗	plop(x);
459		✗	plop(debug_var_multiplicity(x));
460		✗	plop(DEBUG_init_multiplicity[x]);
461		✗	show_line(x);
462		✗	geo_assert_not_reached;
463			}
464		✗	}
465		✗	void show_line(index_t x) {
466		✗	vector<Coeff> res = get_line(x);
467		✗	plop(DEBUG_init_multiplicity[x]);
468		✗	FOR(i, res.size()) {
469		✗	std::cerr<< res[i].a<<" . "<< res[i].index<< " ( mul= " << multiplicity[res[i].index] <<" )\n";
470			}
471		✗	}
472
473		✗	vector<Coeff> get_line(index_t l) {
474		✗	vector<Coeff> c(1, Coeff(l, 1.));
475		✗	vector<Coeff> cM0 = mult_by_M0(c);
476		✗	if (M1.empty()) return cM0;
477		✗	vector<Coeff> cM0M1 = mult_by_M1(cM0);
478		✗	if (M2.M.empty()) return cM0M1;
479		✗	return mult_by_M2(cM0M1);
480			}
481
482		✗	double debug_var_multiplicity(index_t x) {
483		✗	double min_mult = 10000;
484		✗	vector<Coeff> res = get_line(x);
485		✗	FOR(i, res.size()) {
486		✗	min_mult = std::min(
487			min_mult,
488		✗	double(multiplicity[res[i].index]) * std::abs(res[i].a)
489			);
490			}
491		✗	return min_mult;
492			}
493
494			vector<index_t> M0; // indirection map that removes null variables (set when growing a sphere)
495			vector<Coeff> M1; // indirection map + coefficient that produces groups of variables (for boundary & cuts)
496			NaiveSparseMatrix M2; // removes extra DOF when joining boundaries & cuts
497			NaiveSparseMatrix M2t; // transposed of M2: speed up computations
498			vector<int> multiplicity; // the multiplicity constraint for each variable (solver variable, not user variable)
499			index_t kernel_size; // number of final DOF
500			index_t pass; // iteration in filling the constraint matrix
501			vector<int> DEBUG_init_multiplicity; // the multiplicity constraint for each variable (user variable)
502			};
503
504
505
506			struct EXPLORAGRAM_API MatrixMixedConstrainedSolver {
507
508		✗	MatrixMixedConstrainedSolver(index_t p_nb_vars) :
509		✗	M(p_nb_vars),
510		✗	size_(p_nb_vars) {
511		✗	snap_size = 0;
512		✗	}
513
514			/**
515			* \brief Gets the number of variables.
516			* \return the number of variables.
517			*/
518			index_t size() const {
519			return size_;
520			}
521
522			// set multiple (including integer) constraints
523			void set_multiplicity(index_t id, int period) {
524			geo_debug_assert(id < size());
525		✗	M.multiplicity[id] = period;
526		✗	M.DEBUG_init_multiplicity[id] = period;
527		✗	}
528
529			// constraints are included by increasing nuber of variables (from 1 to 3). The last pass (==3) is just for checking.
530			void end_pass(index_t p_pass) {
531		✗	M.end_pass(p_pass);
532		✗	}
533
534			// add a new constraint
535			void begin_constraint() {
536			new_constraint.clear();
537			}
538
539			void add_constraint_coeff(index_t id, double coeff) {
540			geo_debug_assert(id < size());
541		✗	if (coeff != 0.0) {
542		✗	new_constraint.push_back(Coeff(id, coeff));
543			}
544			}
545
546			bool end_constraint() {
547		✗	return M.add_constraint(new_constraint);
548			}
549			vector<Coeff> new_constraint;
550
551			// start / end solving iterations
552		✗	void start_new_iter() {
553		✗	std::cerr << "New LS iteration \n";
554			bool first_iter = V.empty();
555		✗	if (first_iter) {
556		✗	fixed.resize(M.kernel_size, false);
557		✗	V.resize(M.kernel_size);
558			}
559		✗	nlNewContext();
560		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
561		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(M.kernel_size));
562		✗	nlBegin(NL_SYSTEM);
563
564		✗	if (!first_iter) {
565			static const double auto_snap_threshold = 1.;// .005;// 25;
566			double snap_threshold = 1e20;
567			int nb_fixed_var = 0;
568			bool snap_size_has_changed = false;
569		✗	do {
570		✗	std::cerr << "Try enforce mutiplicity equality with snap size " << snap_size << "\n";
571		✗	FOR(pass, 2) {
572			//FOR(i, M.kernel_size) if (M.multiplicity[i]>0)std::cerr << i << " " << M.multiplicity[i] << "\n";
573		✗	FOR(i, M.kernel_size) {
574		✗	nlSetVariable(i, V[i]);
575
576		✗	if (M.multiplicity[i] < snap_size) continue;
577		✗	if (M.multiplicity[i] ==0 ) continue;
578		✗	if (M.M2t.nb_coeffs_in_line(i) == 0) continue;
579		✗	if (pass == 0 && fixed[i]) continue;
580		✗	double val = V[i];
581		✗	double snapped_val = double(M.multiplicity[i]) * nint(V[i] / double(M.multiplicity[i]));
582		✗	double dist = std::abs(snapped_val - val) / double(M.multiplicity[i]);
583		✗	if (pass == 0) {
584		✗	if (dist < snap_threshold
585		✗	&& snap_threshold>auto_snap_threshold) {
586		✗	snap_threshold = std::max(dist + .001, auto_snap_threshold);
587			}
588		✗	continue;
589			}
590
591			//if (!fixed[i]) {
592			// std::cerr << std::setprecision(2);
593			// std::cerr << snap_threshold << " " << dist << " " << val << "(" << M.multiplicity[i] << ")\n";
594			// if (dist < snap_threshold)
595			// std::cerr << snap_threshold << " < " << dist << " < " << 1.- snap_threshold << "\n";
596			//}
597
598		✗	if (dist < snap_threshold && M.multiplicity[i] != 1000) {
599		✗	if (!fixed[i]) nb_fixed_var++;
600			fixed[i] = true;
601		✗	GEO::Logger::out("HexDom") << " i " << i << " snapped_val" << snapped_val << " (V[i] " << V[i] << std::endl;
602		✗	nlSetVariable(i, snapped_val);
603		✗	nlLockVariable(i);
604			}
605			}
606			}
607
608			snap_size_has_changed = false;
609		✗	if (nb_fixed_var == 0) {
610		✗	if (snap_size >= 1) {
611			snap_size_has_changed = true;
612		✗	snap_size /= 2;
613			}
614			}
615		✗	plop(nb_fixed_var);
616			} while (snap_size_has_changed);
617			}
618		✗	nlBegin(NL_MATRIX);
619		✗	}
620
621		✗	bool converged() {
622
623			if (fixed.empty()) return false;
624		✗	FOR(i, M.kernel_size)
625		✗	if (M.multiplicity[i] > 0
626		✗	&& !fixed[i]
627		✗	&& M.M2t.nb_coeffs_in_line(i) > 0) return false;
628			return true;
629			}
630
631		✗	void end_iter() {
632		✗	nlEnd(NL_MATRIX);
633		✗	nlEnd(NL_SYSTEM);
634		✗	nlEnable(NL_VERBOSE);
635		✗	nlSolve();
636		✗	FOR(i, M.kernel_size) {
637		✗	V[i] = nlGetVariable(i);
638			}
639		✗	nlDeleteContext(nlGetCurrent());
640		✗	}
641
642			// set up energy
643		✗	void begin_energy() { nlBegin(NL_ROW); }
644		✗	void add_energy_coeff(index_t id, double coeff) {
645			geo_debug_assert(id < size());
646		✗	if (coeff == 0.0) {
647		✗	return;
648			}
649		✗	vector<Coeff> line = M.get_line(id);
650		✗	FOR(m, line.size()) {
651		✗	nlCoefficient(line[m].index, line[m].a*coeff);
652			}
653			}
654		✗	void add_energy_rhs(double rhs) { nlRightHandSide(rhs); }
655
656		✗	void end_energy() { nlEnd(NL_ROW); }
657
658			// access to the result ;)
659		✗	double value(index_t i) {
660			geo_debug_assert(i < size());
661			double res = 0;
662		✗	vector<Coeff> line = M.get_line(i);
663		✗	FOR(m, line.size()) {
664		✗	res += V[line[m].index] * line[m].a;
665			}
666		✗	return res;
667			}
668
669			double show_var(index_t i) {
670			geo_debug_assert(i < size());
671			double res = 0;
672			vector<Coeff> line = M.get_line(i);
673			std::cerr << " var " << i << " = ";
674			FOR(m, line.size()) {
675			std::cerr << line[m].a << " X " << V[line[m].index] << " + ";
676			res += V[line[m].index] * line[m].a;
677			}
678			GEO::Logger::out("HexDom") << std::endl;
679			return res;
680			}
681
682
683			bool solved;
684			MatrixM M; // Matrix + multiplicity of variables
685			index_t size_;
686			std::vector<double> V; // inner variables
687			std::vector<bool> fixed; // is it fixed already ?
688			int snap_size;
689			};
690
691			}
692
693
694			#endif
695