matrixSetup.h 14.1 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
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

25
26
#ifndef __LINEAR_SOLVING_MATRIX_SETUP__
#define __LINEAR_SOLVING_MATRIX_SETUP__
Pierre Kraemer's avatar
Pierre Kraemer committed
27
28
29
30
31
32
33

namespace CGoGN
{

namespace LinearSolving
{

34
35
36
37
/*******************************************************************************
 * EQUALITY MATRIX : right-hand-side SCALAR
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
38
template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
Pierre Kraemer's avatar
Pierre Kraemer committed
39
40
41
42
class FunctorEquality_PerVertexWeight_Scalar : public FunctorType
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
43
44
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
45
	const VertexAttribute<typename PFP::REAL>& weightTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
46
47
48
49

public:
	FunctorEquality_PerVertexWeight_Scalar(
		LinearSolver<SOLVER_TRAITS>* s,
50
51
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
52
		const VertexAttribute<typename PFP::REAL>& weight
Pierre Kraemer's avatar
Pierre Kraemer committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
	) :	solver(s), indexTable(index), attrTable(attr), weightTable(weight)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		solver->add_coefficient(indexTable[d], 1) ;
		solver->set_right_hand_side(attrTable[d]) ;
		solver->normalize_row(weightTable[d]) ;
		solver->end_row() ;
		return false ;
	}
} ;

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorEquality_UniformWeight_Scalar : public FunctorType
Pierre Kraemer's avatar
Pierre Kraemer committed
69
70
71
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
72
73
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
74
75
76
	float weight ;

public:
Pierre Kraemer's avatar
Pierre Kraemer committed
77
	FunctorEquality_UniformWeight_Scalar(
78
		LinearSolver<SOLVER_TRAITS>* s,
79
		const VertexAttribute<unsigned int>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
80
		const VertexAttribute<ATTR_TYPE>& attr,
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
		float w
	) :	solver(s), indexTable(index), attrTable(attr), weight(w)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		solver->add_coefficient(indexTable[d], 1) ;
		solver->set_right_hand_side(attrTable[d]) ;
		solver->normalize_row(weight) ;
		solver->end_row() ;
		return false ;
	}
} ;

/*******************************************************************************
 * EQUALITY MATRIX : right-hand-side VECTOR + coordinate
 *******************************************************************************/

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
Pierre Kraemer's avatar
Pierre Kraemer committed
101
102
103
104
class FunctorEquality_PerVertexWeight_Vector : public FunctorType
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
105
106
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
107
	const VertexAttribute<typename PFP::REAL>& weightTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
108
109
110
111
112
	unsigned int coord ;

public:
	FunctorEquality_PerVertexWeight_Vector(
		LinearSolver<SOLVER_TRAITS>* s,
113
114
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
115
		const VertexAttribute<typename PFP::REAL>& weight,
Pierre Kraemer's avatar
Pierre Kraemer committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
		unsigned int c
	) :	solver(s), indexTable(index), attrTable(attr), weightTable(weight), coord(c)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		solver->add_coefficient(indexTable[d], 1) ;
		solver->set_right_hand_side((attrTable[d])[coord]) ;
		solver->normalize_row(weightTable[d]) ;
		solver->end_row() ;
		return false ;
	}
} ;

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorEquality_UniformWeight_Vector : public FunctorType
133
134
135
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
136
137
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
138
	float weight ;
139
140
141
	unsigned int coord ;

public:
Pierre Kraemer's avatar
Pierre Kraemer committed
142
	FunctorEquality_UniformWeight_Vector(
143
		LinearSolver<SOLVER_TRAITS>* s,
144
145
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
		float w,
		unsigned int c
	) :	solver(s), indexTable(index), attrTable(attr), weight(w), coord(c)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		solver->add_coefficient(indexTable[d], 1) ;
		solver->set_right_hand_side((attrTable[d])[coord]) ;
		solver->normalize_row(weight) ;
		solver->end_row() ;
		return false ;
	}
} ;

/*******************************************************************************
 * LAPLACIAN TOPO MATRIX : right-hand-side 0
 *******************************************************************************/

template<typename PFP, class SOLVER_TRAITS>
class FunctorLaplacianTopo : public FunctorMap<typename PFP::MAP>
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
171
	const VertexAttribute<unsigned int>& indexTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
172
173
174

public:
	typedef typename PFP::MAP MAP ;
175
	typedef typename PFP::REAL REAL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
176

177
178
179
	FunctorLaplacianTopo(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
180
		const VertexAttribute<unsigned int>& index
181
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index)
Pierre Kraemer's avatar
Pierre Kraemer committed
182
	{}
183

Pierre Kraemer's avatar
Pierre Kraemer committed
184
185
	bool operator()(Dart d)
	{
186
187
		solver->begin_row() ;
		REAL aii = 0 ;
188
189
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
190
191
192
193
		{
			REAL aij = 1 ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
194
		}
195
196
197
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
		solver->set_right_hand_side(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
198
199
200
201
202
		solver->end_row() ;
		return false ;
	}
} ;

203
204
205
206
/*******************************************************************************
 * LAPLACIAN TOPO MATRIX : right-hand-side SCALAR
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
207
template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
208
class FunctorLaplacianTopoRHS_Scalar : public FunctorMap<typename PFP::MAP>
Pierre Kraemer's avatar
Pierre Kraemer committed
209
210
211
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
212
213
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
214
215
216
217
218
219
220
221

public:
	typedef typename PFP::MAP MAP ;
	typedef typename PFP::REAL REAL ;

	FunctorLaplacianTopoRHS_Scalar(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
222
223
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr
224
225
226
227
228
229
230
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index), attrTable(attr)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		REAL aii = 0 ;
231
232
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
233
234
235
236
		{
			REAL aij = 1 ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
237
		}
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
		solver->set_right_hand_side(attrTable[d]) ;
		solver->end_row() ;
		return false ;
	}
} ;

/*******************************************************************************
 * LAPLACIAN TOPO MATRIX : right-hand-side VECTOR + coordinate
 *******************************************************************************/

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorLaplacianTopoRHS_Vector : public FunctorMap<typename PFP::MAP>
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
255
256
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
257
258
259
260
	unsigned int coord ;

public:
	typedef typename PFP::MAP MAP ;
261
	typedef typename PFP::REAL REAL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
262

263
264
265
	FunctorLaplacianTopoRHS_Vector(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
266
267
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
268
269
		unsigned int c
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index), attrTable(attr), coord(c)
Pierre Kraemer's avatar
Pierre Kraemer committed
270
	{}
271

Pierre Kraemer's avatar
Pierre Kraemer committed
272
273
	bool operator()(Dart d)
	{
274
275
		solver->begin_row() ;
		REAL aii = 0 ;
276
277
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
278
279
280
281
		{
			REAL aij = 1 ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
282
		}
283
284
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
285
286
287
288
289
290
		solver->set_right_hand_side((attrTable[d])[coord]) ;
		solver->end_row() ;
		return false ;
	}
} ;

291
292
293
294
295
296
/*******************************************************************************
 * LAPLACIAN COTAN MATRIX : right-hand-side 0
 *******************************************************************************/

template<typename PFP, class SOLVER_TRAITS>
class FunctorLaplacianCotan : public FunctorMap<typename PFP::MAP>
Pierre Kraemer's avatar
Pierre Kraemer committed
297
{
298
299
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
300
	const VertexAttribute<unsigned int>& indexTable ;
301
302
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
303
304
305
306
307
308
309
310

public:
	typedef typename PFP::MAP MAP ;
	typedef typename PFP::REAL REAL ;

	FunctorLaplacianCotan(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
311
		const VertexAttribute<unsigned int>& index,
312
313
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea
314
315
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index), edgeWeight(eWeight), vertexArea(vArea)
	{}
Pierre Kraemer's avatar
Pierre Kraemer committed
316

317
318
319
320
321
	bool operator()(Dart d)
	{
		solver->begin_row() ;
		REAL vArea = vertexArea[d] ;
		REAL aii = 0 ;
322
323
324
		Dart it = d ;
//		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
//		for(Dart it = t.begin(); it != t.end(); it = t.next())
325
326
327
328
329
		do
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
330
			it = this->m_map.phi2_1(it) ;
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
		} while(it != d) ;
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
		solver->set_right_hand_side(0) ;
		solver->end_row() ;
		return false ;
	}
} ;

/*******************************************************************************
 * LAPLACIAN COTAN MATRIX : right-hand-side SCALAR
 *******************************************************************************/

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorLaplacianCotanRHS_Scalar : public FunctorMap<typename PFP::MAP>
Pierre Kraemer's avatar
Pierre Kraemer committed
346
{
347
348
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
349
	const VertexAttribute<unsigned int>& indexTable ;
350
351
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
352
	const VertexAttribute<ATTR_TYPE>& attrTable ;
353
354
355
356
357
358
359
360

public:
	typedef typename PFP::MAP MAP ;
	typedef typename PFP::REAL REAL ;

	FunctorLaplacianCotanRHS_Scalar(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
361
		const VertexAttribute<unsigned int>& index,
362
363
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea,
364
		const VertexAttribute<ATTR_TYPE>& attr
365
366
367
368
369
370
371
372
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		REAL vArea = vertexArea[d] ;
		REAL aii = 0 ;
373
374
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
375
376
377
378
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
379
		}
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
		solver->set_right_hand_side(attrTable[d]) ;
		solver->end_row() ;
		return false ;
	}
} ;

/*******************************************************************************
 * LAPLACIAN COTAN MATRIX : right-hand-side VECTOR + coordinate
 *******************************************************************************/

template<typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorLaplacianCotanRHS_Vector : public FunctorMap<typename PFP::MAP>
{
protected:
	LinearSolver<SOLVER_TRAITS>* solver ;
397
	const VertexAttribute<unsigned int>& indexTable ;
398
399
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
400
	const VertexAttribute<ATTR_TYPE>& attrTable ;
401
402
403
404
405
406
407
408
409
	unsigned int coord ;

public:
	typedef typename PFP::MAP MAP ;
	typedef typename PFP::REAL REAL ;

	FunctorLaplacianCotanRHS_Vector(
		MAP& m,
		LinearSolver<SOLVER_TRAITS>* s,
410
		const VertexAttribute<unsigned int>& index,
411
412
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea,
413
		const VertexAttribute<ATTR_TYPE>& attr,
414
415
416
417
418
419
420
421
422
		unsigned int c
	) :	FunctorMap<MAP>(m), solver(s), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr), coord(c)
	{}

	bool operator()(Dart d)
	{
		solver->begin_row() ;
		REAL vArea = vertexArea[d] ;
		REAL aii = 0 ;
423
424
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
425
426
427
428
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
			solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
429
		}
430
431
432
433
434
435
436
		solver->add_coefficient(indexTable[d], -aii) ;
		solver->normalize_row() ;
		solver->set_right_hand_side((attrTable[d])[coord]) ;
		solver->end_row() ;
		return false ;
	}
} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
437

438
} // namespace LinearSolving
Pierre Kraemer's avatar
Pierre Kraemer committed
439

440
} // namespace CGoGN
Pierre Kraemer's avatar
Pierre Kraemer committed
441
442

#endif