basic.h 12.8 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
53
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
54
	const CellMarker<typename PFP::MAP, VERTEX>& freeMarker,
Pierre Kraemer's avatar
Pierre Kraemer committed
55
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
68
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
69
	const CellMarker<typename PFP::MAP, VERTEX>& freeMarker,
Pierre Kraemer's avatar
Pierre Kraemer committed
70
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
88 89 90
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& attr,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
109 110
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
130 131 132
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& attr,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& 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,
Pierre Kraemer's avatar
Pierre Kraemer committed
152 153
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
178
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& 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,
204
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
205
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
240
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
241
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
281
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
282 283
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& 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,
310
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
311 312 313
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& vertexArea,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
349
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
350 351 352
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeWeight,
	const VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& vertexArea,
	const VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
393
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
394
	VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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,
405
	const VertexAttribute<unsigned int, typename PFP::MAP::IMPL>& index,
Pierre Kraemer's avatar
Pierre Kraemer committed
406
	VertexAttribute<ATTR_TYPE, typename PFP::MAP::IMPL>& 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