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