basic.h 12.6 KB
Newer Older
Pierre Kraemer's avatar
Pierre Kraemer committed
1
2
3
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
* version 0.1                                                                  *
4
* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
Pierre Kraemer's avatar
Pierre Kraemer committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
*                                                                              *
* This library is free software; you can redistribute it and/or modify it      *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 of the License, or (at your     *
* option) any later version.                                                   *
*                                                                              *
* This library is distributed in the hope that it will be useful, but WITHOUT  *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
* for more details.                                                            *
*                                                                              *
* You should have received a copy of the GNU Lesser General Public License     *
* along with this library; if not, write to the Free Software Foundation,      *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
*                                                                              *
20
* Web site: http://cgogn.unistra.fr/                                           *
Pierre Kraemer's avatar
Pierre Kraemer committed
21
22
23
24
25
26
27
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#ifndef __LINEAR_SOLVING_BASIC__
#define __LINEAR_SOLVING_BASIC__

28
#include "NL/nl.h"
29
#include "Topology/generic/traversor/traversorCell.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
30
31
32
33
34
35
36

namespace CGoGN
{

namespace LinearSolving
{

37
38
39
40
41
42
43
44
45
template <typename CoeffType>
struct Coeff
{
	unsigned int index;
	CoeffType value;
	Coeff(unsigned int i, CoeffType v) : index(i), value(v)
	{}
} ;

46
47
48
49
/*******************************************************************************
 * VARIABLES SETUP
 *******************************************************************************/

50
template <typename PFP, typename ATTR_TYPE>
51
52
void setupVariables(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
53
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
54
	const CellMarker<typename PFP::MAP, VERTEX>& freeMarker,
Sylvain Thery's avatar
Sylvain Thery committed
55
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
56
{
57
58
59
60
61
62
	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlSetVariable(index[d], attr[d]);
		if(!freeMarker.isMarked(d))
			nlLockVariable(index[d]);
	});
Pierre Kraemer's avatar
Pierre Kraemer committed
63
64
}

65
template <typename PFP, typename ATTR_TYPE>
66
67
void setupVariables(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
68
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
69
	const CellMarker<typename PFP::MAP, VERTEX>& freeMarker,
Sylvain Thery's avatar
Sylvain Thery committed
70
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
71
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
72
{
73
74
75
76
77
78
	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlSetVariable(index[d], (attr[d])[coord]);
		if(!freeMarker.isMarked(d))
			nlLockVariable(index[d]);
	});
Pierre Kraemer's avatar
Pierre Kraemer committed
79
80
}

81
82
83
84
/*******************************************************************************
 * MATRIX SETUP : EQUALITY
 *******************************************************************************/

85
template <typename PFP, typename ATTR_TYPE>
86
87
void addRowsRHS_Equality(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
88
89
90
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP>& weight)
91
{
92
	nlEnable(NL_NORMALIZE_ROWS) ;
93
94
95
96
97
98
99
100
101
102

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, attr[d]) ;
		nlRowParameterd(NL_ROW_SCALING, weight[d]) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(index[d], 1) ;
		nlEnd(NL_ROW) ;
	});

103
	nlDisable(NL_NORMALIZE_ROWS) ;
104
105
}

106
template <typename PFP, typename ATTR_TYPE>
107
108
void addRowsRHS_Equality(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
109
110
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
Pierre Kraemer's avatar
Pierre Kraemer committed
111
112
	float weight)
{
113
	nlEnable(NL_NORMALIZE_ROWS) ;
114
115
116
117
118
119
120
121
122
123

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, attr[d]) ;
		nlRowParameterd(NL_ROW_SCALING, weight) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(index[d], 1) ;
		nlEnd(NL_ROW) ;
	});

124
	nlDisable(NL_NORMALIZE_ROWS) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
125
126
}

127
template <typename PFP, typename ATTR_TYPE>
Pierre Kraemer's avatar
Pierre Kraemer committed
128
129
void addRowsRHS_Equality(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
130
131
132
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP>& weight,
Pierre Kraemer's avatar
Pierre Kraemer committed
133
134
	unsigned int coord)
{
135
	nlEnable(NL_NORMALIZE_ROWS) ;
136
137
138
139
140
141
142
143
144
145

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attr[d])[coord]) ;
		nlRowParameterd(NL_ROW_SCALING, weight[d]) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(index[d], 1) ;
		nlEnd(NL_ROW) ;
	});

146
	nlDisable(NL_NORMALIZE_ROWS) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
147
148
}

149
template <typename PFP, typename ATTR_TYPE>
Pierre Kraemer's avatar
Pierre Kraemer committed
150
151
void addRowsRHS_Equality(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
152
153
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
Pierre Kraemer's avatar
Pierre Kraemer committed
154
	float weight,
155
156
	unsigned int coord)
{
157
	nlEnable(NL_NORMALIZE_ROWS) ;
158
159
160
161
162
163
164
165
166
167

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attr[d])[coord]) ;
		nlRowParameterd(NL_ROW_SCALING, weight) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(index[d], 1) ;
		nlEnd(NL_ROW) ;
	});

168
	nlDisable(NL_NORMALIZE_ROWS) ;
169
170
171
172
173
174
}

/*******************************************************************************
 * MATRIX SETUP : LAPLACIAN TOPO
 *******************************************************************************/

175
template <typename PFP>
176
177
void addRows_Laplacian_Topo(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
178
	const VertexAttribute<unsigned int, typename PFP::MAP>& index)
179
{
180
	nlEnable(NL_NORMALIZE_ROWS) ;
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ;
		nlBegin(NL_ROW);
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = 1 ;
			aii += aij ;
			nlCoefficient(index[m.phi1(it)], aij) ;
		}
		nlCoefficient(index[d], -aii) ;
		nlEnd(NL_ROW) ;
	});

198
	nlDisable(NL_NORMALIZE_ROWS) ;
199
200
}

201
template <typename PFP, typename ATTR_TYPE>
202
203
void addRowsRHS_Laplacian_Topo(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
204
205
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
206
{
207
	nlEnable(NL_NORMALIZE_ROWS) ;
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		std::vector<Coeff<typename PFP::REAL> > coeffs ;
		coeffs.reserve(12) ;

		typename PFP::REAL norm2 = 0 ;
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = 1 ;
			aii += aij ;
			coeffs.push_back(Coeff<typename PFP::REAL>(index[m.phi1(it)], aij)) ;
			norm2 += aij * aij ;
		}
		coeffs.push_back(Coeff<typename PFP::REAL>(index[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, attr[d] * sqrt(norm2)) ;
		nlBegin(NL_ROW);
		for(unsigned int i = 0; i < coeffs.size(); ++i)
			nlCoefficient(coeffs[i].index, coeffs[i].value) ;
		nlEnd(NL_ROW) ;
	});

234
	nlDisable(NL_NORMALIZE_ROWS) ;
235
236
}

237
template <typename PFP, typename ATTR_TYPE>
238
239
void addRowsRHS_Laplacian_Topo(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
240
241
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
242
243
	unsigned int coord)
{
244
	nlEnable(NL_NORMALIZE_ROWS) ;
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		std::vector<Coeff<typename PFP::REAL> > coeffs ;
		coeffs.reserve(12) ;

		typename PFP::REAL norm2 = 0 ;
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = 1 ;
			aii += aij ;
			coeffs.push_back(Coeff<typename PFP::REAL>(index[m.phi1(it)], aij)) ;
			norm2 += aij * aij ;
		}
		coeffs.push_back(Coeff<typename PFP::REAL>(index[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attr[d])[coord] * sqrt(norm2)) ;
		nlBegin(NL_ROW);
		for(unsigned int i = 0; i < coeffs.size(); ++i)
			nlCoefficient(coeffs[i].index, coeffs[i].value) ;
		nlEnd(NL_ROW) ;
	});

271
	nlDisable(NL_NORMALIZE_ROWS) ;
272
273
274
275
276
277
}

/*******************************************************************************
 * MATRIX SETUP : LAPLACIAN COTAN
 *******************************************************************************/

278
template <typename PFP>
279
280
void addRows_Laplacian_Cotan(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
281
282
283
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP>& vertexArea)
284
{
285
	nlEnable(NL_NORMALIZE_ROWS) ;
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ;
		nlBegin(NL_ROW);
		typename PFP::REAL vArea = vertexArea[d] ;
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			nlCoefficient(index[m.phi1(it)], aij) ;
		}
		nlCoefficient(index[d], -aii) ;
		nlEnd(NL_ROW) ;
	});

304
	nlDisable(NL_NORMALIZE_ROWS) ;
305
}
Pierre Kraemer's avatar
Pierre Kraemer committed
306

307
template <typename PFP, typename ATTR_TYPE>
308
309
void addRowsRHS_Laplacian_Cotan(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
310
311
312
313
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP>& vertexArea,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
314
{
315
	nlEnable(NL_NORMALIZE_ROWS) ;
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342

	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		std::vector<Coeff<typename PFP::REAL> > coeffs ;
		coeffs.reserve(12) ;

		typename PFP::REAL vArea = vertexArea[d] ;
		typename PFP::REAL norm2 = 0 ;
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			coeffs.push_back(Coeff<typename PFP::REAL>(index[m.phi1(it)], aij)) ;
			norm2 += aij * aij ;
		}
		coeffs.push_back(Coeff<typename PFP::REAL>(index[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, attr[d] * sqrt(norm2)) ;
		nlBegin(NL_ROW);
		for(unsigned int i = 0; i < coeffs.size(); ++i)
			nlCoefficient(coeffs[i].index, coeffs[i].value) ;
		nlEnd(NL_ROW) ;
	});

343
	nlDisable(NL_NORMALIZE_ROWS) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
344
345
}

346
template <typename PFP, typename ATTR_TYPE>
347
348
void addRowsRHS_Laplacian_Cotan(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
349
350
351
352
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP>& vertexArea,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
353
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
354
{
355
	nlEnable(NL_NORMALIZE_ROWS) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
356

357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		std::vector<Coeff<typename PFP::REAL> > coeffs ;
		coeffs.reserve(12) ;

		typename PFP::REAL vArea = vertexArea[d] ;
		typename PFP::REAL norm2 = 0 ;
		typename PFP::REAL aii = 0 ;
		Traversor2VE<typename PFP::MAP> t(m, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
		{
			typename PFP::REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			coeffs.push_back(Coeff<typename PFP::REAL>(index[m.phi1(it)], aij)) ;
			norm2 += aij * aij ;
		}
		coeffs.push_back(Coeff<typename PFP::REAL>(index[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attr[d])[coord] * sqrt(norm2)) ;
		nlBegin(NL_ROW);
		for(unsigned int i = 0; i < coeffs.size(); ++i)
			nlCoefficient(coeffs[i].index, coeffs[i].value) ;
		nlEnd(NL_ROW) ;
	});

383
	nlDisable(NL_NORMALIZE_ROWS) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
384
385
}

386
387
388
389
/*******************************************************************************
 * GET RESULTS
 *******************************************************************************/

390
template <typename PFP, typename ATTR_TYPE>
391
392
void getResult(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
393
394
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
395
{
396
397
398
399
	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		attr[d] = nlGetVariable(index[d]) ;
	});
Pierre Kraemer's avatar
Pierre Kraemer committed
400
401
}

402
template <typename PFP, typename ATTR_TYPE>
403
404
void getResult(
	typename PFP::MAP& m,
Sylvain Thery's avatar
Sylvain Thery committed
405
406
	const VertexAttribute<unsigned int, typename PFP::MAP>& index,
	VertexAttribute<ATTR_TYPE, typename PFP::MAP>& attr,
407
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
408
{
409
410
411
412
	foreach_cell<VERTEX>(m, [&] (Dart d)
	{
		(attr[d])[coord] = nlGetVariable(index[d]) ;
	});
Pierre Kraemer's avatar
Pierre Kraemer committed
413
414
415
416
417
418
419
}

} // namespace LinearSolving

} // namespace CGoGN

#endif