matrixSetup.h 13.9 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 38 39 40 41 42
template <typename CoeffType>
struct Coeff
{
	unsigned int index;
	CoeffType value;
	Coeff(unsigned int i, CoeffType v) : index(i), value(v)
	{}
} ;

43 44 45 46
/*******************************************************************************
 * EQUALITY MATRIX : right-hand-side SCALAR
 *******************************************************************************/

47
template<typename PFP, typename ATTR_TYPE>
Pierre Kraemer's avatar
Pierre Kraemer committed
48 49 50
class FunctorEquality_PerVertexWeight_Scalar : public FunctorType
{
protected:
51 52
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
53
	const VertexAttribute<typename PFP::REAL>& weightTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
54 55 56

public:
	FunctorEquality_PerVertexWeight_Scalar(
57 58
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
59
		const VertexAttribute<typename PFP::REAL>& weight
60
	) :	indexTable(index), attrTable(attr), weightTable(weight)
Pierre Kraemer's avatar
Pierre Kraemer committed
61 62 63 64
	{}

	bool operator()(Dart d)
	{
65 66 67 68 69
		nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[d]) ;
		nlRowParameterd(NL_ROW_SCALING, weightTable[d]) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(indexTable[d], 1) ;
		nlEnd(NL_ROW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
70 71 72 73
		return false ;
	}
} ;

74
template<typename PFP, typename ATTR_TYPE>
Pierre Kraemer's avatar
Pierre Kraemer committed
75
class FunctorEquality_UniformWeight_Scalar : public FunctorType
Pierre Kraemer's avatar
Pierre Kraemer committed
76 77
{
protected:
78 79
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
80 81 82
	float weight ;

public:
Pierre Kraemer's avatar
Pierre Kraemer committed
83
	FunctorEquality_UniformWeight_Scalar(
84
		const VertexAttribute<unsigned int>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
85
		const VertexAttribute<ATTR_TYPE>& attr,
86
		float w
87
	) :	indexTable(index), attrTable(attr), weight(w)
88 89 90 91
	{}

	bool operator()(Dart d)
	{
92 93 94 95 96
		nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[d]) ;
		nlRowParameterd(NL_ROW_SCALING, weight) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(indexTable[d], 1) ;
		nlEnd(NL_ROW) ;
97 98 99 100 101 102 103 104
		return false ;
	}
} ;

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

105
template<typename PFP, typename ATTR_TYPE>
Pierre Kraemer's avatar
Pierre Kraemer committed
106 107 108
class FunctorEquality_PerVertexWeight_Vector : public FunctorType
{
protected:
109 110
	const VertexAttribute<unsigned int>& indexTable ;
	const VertexAttribute<ATTR_TYPE>& attrTable ;
111
	const VertexAttribute<typename PFP::REAL>& weightTable ;
Pierre Kraemer's avatar
Pierre Kraemer committed
112 113 114 115
	unsigned int coord ;

public:
	FunctorEquality_PerVertexWeight_Vector(
116 117
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
118
		const VertexAttribute<typename PFP::REAL>& weight,
Pierre Kraemer's avatar
Pierre Kraemer committed
119
		unsigned int c
120
	) :	indexTable(index), attrTable(attr), weightTable(weight), coord(c)
Pierre Kraemer's avatar
Pierre Kraemer committed
121 122 123 124
	{}

	bool operator()(Dart d)
	{
125 126 127 128 129
		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[d])[coord]) ;
		nlRowParameterd(NL_ROW_SCALING, weightTable[d]) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(indexTable[d], 1) ;
		nlEnd(NL_ROW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
130 131 132 133
		return false ;
	}
} ;

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

public:
Pierre Kraemer's avatar
Pierre Kraemer committed
144
	FunctorEquality_UniformWeight_Vector(
145 146
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
147 148
		float w,
		unsigned int c
149
	) :	indexTable(index), attrTable(attr), weight(w), coord(c)
150 151 152 153
	{}

	bool operator()(Dart d)
	{
154 155 156 157 158
		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[d])[coord]) ;
		nlRowParameterd(NL_ROW_SCALING, weight) ;
		nlBegin(NL_ROW) ;
		nlCoefficient(indexTable[d], 1) ;
		nlEnd(NL_ROW) ;
159 160 161 162 163 164 165 166
		return false ;
	}
} ;

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

167
template<typename PFP>
168 169 170
class FunctorLaplacianTopo : public FunctorMap<typename PFP::MAP>
{
protected:
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
	FunctorLaplacianTopo(
		MAP& m,
179
		const VertexAttribute<unsigned int>& index
180
	) :	FunctorMap<MAP>(m), indexTable(index)
Pierre Kraemer's avatar
Pierre Kraemer committed
181
	{}
182

Pierre Kraemer's avatar
Pierre Kraemer committed
183 184
	bool operator()(Dart d)
	{
185 186
		nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ;
		nlBegin(NL_ROW);
187
		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
		{
			REAL aij = 1 ;
			aii += aij ;
193
			nlCoefficient(indexTable[this->m_map.phi1(it)], aij) ;
194
		}
195 196
		nlCoefficient(indexTable[d], -aii) ;
		nlEnd(NL_ROW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
197 198 199 200
		return false ;
	}
} ;

201 202 203 204
/*******************************************************************************
 * LAPLACIAN TOPO MATRIX : right-hand-side SCALAR
 *******************************************************************************/

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

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

	FunctorLaplacianTopoRHS_Scalar(
		MAP& m,
218 219
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr
220
	) :	FunctorMap<MAP>(m), indexTable(index), attrTable(attr)
221 222 223 224
	{}

	bool operator()(Dart d)
	{
225 226 227 228
		std::vector<Coeff<REAL> > coeffs ;
		coeffs.reserve(12) ;

		REAL norm2 = 0 ;
229
		REAL aii = 0 ;
230 231
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
232 233 234
		{
			REAL aij = 1 ;
			aii += aij ;
235 236
			coeffs.push_back(Coeff<REAL>(indexTable[this->m_map.phi1(it)], aij)) ;
			norm2 += aij * aij ;
237
		}
238 239 240 241 242 243 244 245 246
		coeffs.push_back(Coeff<REAL>(indexTable[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[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) ;

247 248 249 250 251 252 253 254
		return false ;
	}
} ;

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

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

public:
	typedef typename PFP::MAP MAP ;
265
	typedef typename PFP::REAL REAL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
266

267 268
	FunctorLaplacianTopoRHS_Vector(
		MAP& m,
269 270
		const VertexAttribute<unsigned int>& index,
		const VertexAttribute<ATTR_TYPE>& attr,
271
		unsigned int c
272
	) :	FunctorMap<MAP>(m), indexTable(index), attrTable(attr), coord(c)
Pierre Kraemer's avatar
Pierre Kraemer committed
273
	{}
274

Pierre Kraemer's avatar
Pierre Kraemer committed
275 276
	bool operator()(Dart d)
	{
277 278 279 280
		std::vector<Coeff<REAL> > coeffs ;
		coeffs.reserve(12) ;

		REAL norm2 = 0 ;
281
		REAL aii = 0 ;
282 283
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
284 285 286
		{
			REAL aij = 1 ;
			aii += aij ;
287 288
			coeffs.push_back(Coeff<REAL>(indexTable[this->m_map.phi1(it)], aij)) ;
			norm2 += aij * aij ;
289
		}
290 291 292 293 294 295 296 297 298
		coeffs.push_back(Coeff<REAL>(indexTable[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[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) ;

Pierre Kraemer's avatar
Pierre Kraemer committed
299 300 301 302
		return false ;
	}
} ;

303 304 305 306
/*******************************************************************************
 * LAPLACIAN COTAN MATRIX : right-hand-side 0
 *******************************************************************************/

307
template<typename PFP>
308
class FunctorLaplacianCotan : public FunctorMap<typename PFP::MAP>
Pierre Kraemer's avatar
Pierre Kraemer committed
309
{
310
protected:
311
	const VertexAttribute<unsigned int>& indexTable ;
312 313
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
314 315 316 317 318 319 320

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

	FunctorLaplacianCotan(
		MAP& m,
321
		const VertexAttribute<unsigned int>& index,
322 323
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea
324
	) :	FunctorMap<MAP>(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea)
325
	{}
Pierre Kraemer's avatar
Pierre Kraemer committed
326

327 328
	bool operator()(Dart d)
	{
329 330
		nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ;
		nlBegin(NL_ROW);
331 332
		REAL vArea = vertexArea[d] ;
		REAL aii = 0 ;
333 334
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
335 336 337
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
338 339 340 341 342
			nlCoefficient(indexTable[this->m_map.phi1(it)], aij) ;
		}
		nlCoefficient(indexTable[d], -aii) ;
		nlEnd(NL_ROW) ;

343 344 345 346 347 348 349 350
		return false ;
	}
} ;

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

351
template<typename PFP, typename ATTR_TYPE>
352
class FunctorLaplacianCotanRHS_Scalar : public FunctorMap<typename PFP::MAP>
Pierre Kraemer's avatar
Pierre Kraemer committed
353
{
354
protected:
355
	const VertexAttribute<unsigned int>& indexTable ;
356 357
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
358
	const VertexAttribute<ATTR_TYPE>& attrTable ;
359 360 361 362 363 364 365

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

	FunctorLaplacianCotanRHS_Scalar(
		MAP& m,
366
		const VertexAttribute<unsigned int>& index,
367 368
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea,
369
		const VertexAttribute<ATTR_TYPE>& attr
370
	) :	FunctorMap<MAP>(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr)
371 372 373 374
	{}

	bool operator()(Dart d)
	{
375 376 377
		std::vector<Coeff<REAL> > coeffs ;
		coeffs.reserve(12) ;

378
		REAL vArea = vertexArea[d] ;
379
		REAL norm2 = 0 ;
380
		REAL aii = 0 ;
381 382
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
383 384 385
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
386 387
			coeffs.push_back(Coeff<REAL>(indexTable[this->m_map.phi1(it)], aij)) ;
			norm2 += aij * aij ;
388
		}
389 390 391 392 393 394 395 396 397
		coeffs.push_back(Coeff<REAL>(indexTable[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[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) ;

398 399 400 401 402 403 404 405
		return false ;
	}
} ;

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

406
template<typename PFP, typename ATTR_TYPE>
407 408 409
class FunctorLaplacianCotanRHS_Vector : public FunctorMap<typename PFP::MAP>
{
protected:
410
	const VertexAttribute<unsigned int>& indexTable ;
411 412
	const EdgeAttribute<typename PFP::REAL>& edgeWeight ;
	const VertexAttribute<typename PFP::REAL>& vertexArea ;
413
	const VertexAttribute<ATTR_TYPE>& attrTable ;
414 415 416 417 418 419 420 421
	unsigned int coord ;

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

	FunctorLaplacianCotanRHS_Vector(
		MAP& m,
422
		const VertexAttribute<unsigned int>& index,
423 424
		const EdgeAttribute<typename PFP::REAL>& eWeight,
		const VertexAttribute<typename PFP::REAL>& vArea,
425
		const VertexAttribute<ATTR_TYPE>& attr,
426
		unsigned int c
427
	) :	FunctorMap<MAP>(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr), coord(c)
428 429 430 431
	{}

	bool operator()(Dart d)
	{
432 433 434
		std::vector<Coeff<REAL> > coeffs ;
		coeffs.reserve(12) ;

435
		REAL vArea = vertexArea[d] ;
436
		REAL norm2 = 0 ;
437
		REAL aii = 0 ;
438 439
		Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
440 441 442
		{
			REAL aij = edgeWeight[it] / vArea ;
			aii += aij ;
443 444
			coeffs.push_back(Coeff<REAL>(indexTable[this->m_map.phi1(it)], aij)) ;
			norm2 += aij * aij ;
445
		}
446 447 448 449 450 451 452 453 454
		coeffs.push_back(Coeff<REAL>(indexTable[d], -aii)) ;
		norm2 += aii * aii ;

		nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[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) ;

455 456 457
		return false ;
	}
} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
458

459
} // namespace LinearSolving
Pierre Kraemer's avatar
Pierre Kraemer committed
460

461
} // namespace CGoGN
Pierre Kraemer's avatar
Pierre Kraemer committed
462 463

#endif