basic.h 10.7 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 28
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#ifndef __LINEAR_SOLVING_BASIC__
#define __LINEAR_SOLVING_BASIC__

#include "OpenNL/linear_solver.h"
29
#include "Algo/LinearSolving/variablesSetup.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
30 31 32 33 34 35 36 37
#include "Algo/LinearSolving/matrixSetup.h"

namespace CGoGN
{

namespace LinearSolving
{

38 39 40
/*******************************************************************************
 * SOLVER INIT
 *******************************************************************************/
Pierre Kraemer's avatar
Pierre Kraemer committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

template <class SOLVER_TRAITS>
void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool least_squares, bool direct)
{
	// create the solver if it does not exist already
	if(s == NULL)
		s = new LinearSolver<SOLVER_TRAITS>(nb_variables) ;

	// if the number of variables has changed, re-create the solver
	else if(s->nb_variables() != nb_variables)
	{
		delete s ;
		s = new LinearSolver<SOLVER_TRAITS>(nb_variables) ;
	}

	s->set_least_squares(least_squares) ;
	s->set_direct(direct) ;
}

60 61 62 63
/*******************************************************************************
 * VARIABLES SETUP
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
64
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
65 66 67
void setupVariables(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
68 69 70
	const AttributeHandler<unsigned int, VERTEX>& index,
	CellMarker<VERTEX>& lm,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
71
{
72
	FunctorMeshToSolver_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(s, index, lm, attr) ;
73
	m.template foreach_orbit<VERTEX>(fmts) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
74 75 76
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
77 78 79
void setupVariables(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
80 81 82
	const AttributeHandler<unsigned int, VERTEX>& index,
	CellMarker<VERTEX>& lm,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
83
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
84
{
85
	FunctorMeshToSolver_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(s, index, lm, attr, coord) ;
86
	m.template foreach_orbit<VERTEX>(fmts) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
87 88
}

89 90 91 92
/*******************************************************************************
 * START MATRIX DEFINITION
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
93 94 95 96 97 98
template <class SOLVER_TRAITS>
void startMatrix(LinearSolver<SOLVER_TRAITS>* s)
{
	s->begin_system() ;
}

99 100 101 102 103 104 105 106
/*******************************************************************************
 * MATRIX SETUP : EQUALITY
 *******************************************************************************/

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
107 108 109
	const AttributeHandler<unsigned int, VERTEX>& index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
	const AttributeHandler<typename PFP::REAL, VERTEX>& weight)
110
{
Pierre Kraemer's avatar
Pierre Kraemer committed
111
	FunctorEquality_PerVertexWeight_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, weight) ;
112
	m.template foreach_orbit<VERTEX>(ec) ;
113 114 115 116 117 118
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
119 120
	const AttributeHandler<unsigned int, VERTEX>& index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
Pierre Kraemer's avatar
Pierre Kraemer committed
121 122 123
	float weight)
{
	FunctorEquality_UniformWeight_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, weight) ;
124
	m.template foreach_orbit<VERTEX>(ec) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
125 126 127 128 129 130
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
131 132 133
	const AttributeHandler<unsigned int, VERTEX>& index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
	const AttributeHandler<typename PFP::REAL, VERTEX>& weight,
Pierre Kraemer's avatar
Pierre Kraemer committed
134 135 136
	unsigned int coord)
{
	FunctorEquality_PerVertexWeight_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, weight, coord) ;
137
	m.template foreach_orbit<VERTEX>(ec) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
138 139 140 141 142 143
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
144 145
	const AttributeHandler<unsigned int, VERTEX>& index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
Pierre Kraemer's avatar
Pierre Kraemer committed
146
	float weight,
147 148
	unsigned int coord)
{
Pierre Kraemer's avatar
Pierre Kraemer committed
149
	FunctorEquality_UniformWeight_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, weight, coord) ;
150
	m.template foreach_orbit<VERTEX>(ec) ;
151 152 153 154 155 156
}

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

157
template <typename PFP, class SOLVER_TRAITS>
158 159 160
void addRows_Laplacian_Topo(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
161
	const AttributeHandler<unsigned int, VERTEX> index)
162 163
{
	FunctorLaplacianTopo<PFP, SOLVER_TRAITS> flt(m, s, index) ;
164
	m.template foreach_orbit<VERTEX>(flt) ;
165 166 167 168 169 170
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Laplacian_Topo(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
171 172
	const AttributeHandler<unsigned int, VERTEX> index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
173
{
174
	FunctorLaplacianTopoRHS_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr) ;
175
	m.template foreach_orbit<VERTEX>(flt) ;
176 177
}

178 179 180 181
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Laplacian_Topo(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
182 183
	const AttributeHandler<unsigned int, VERTEX> index,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
184 185 186
	unsigned int coord)
{
	FunctorLaplacianTopoRHS_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr, coord) ;
187
	m.template foreach_orbit<VERTEX>(flt) ;
188 189 190 191 192 193
}

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

194
template <typename PFP, class SOLVER_TRAITS>
195 196 197
void addRows_Laplacian_Cotan(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
198 199 200
	const AttributeHandler<unsigned int, VERTEX> index,
	const AttributeHandler<typename PFP::REAL, EDGE>& edgeWeight,
	const AttributeHandler<typename PFP::REAL, VERTEX>& vertexArea)
201
{
202
	FunctorLaplacianCotan<PFP, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea) ;
203
	m.template foreach_orbit<VERTEX>(flc) ;
204
}
Pierre Kraemer's avatar
Pierre Kraemer committed
205 206

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
207 208 209
void addRowsRHS_Laplacian_Cotan(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
210 211 212 213
	const AttributeHandler<unsigned int, VERTEX> index,
	const AttributeHandler<typename PFP::REAL, EDGE>& edgeWeight,
	const AttributeHandler<typename PFP::REAL, VERTEX>& vertexArea,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
214
{
215
	FunctorLaplacianCotanRHS_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea, attr) ;
216
	m.template foreach_orbit<VERTEX>(flc) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
217 218 219
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
220 221 222
void addRowsRHS_Laplacian_Cotan(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
223 224 225 226
	const AttributeHandler<unsigned int, VERTEX> index,
	const AttributeHandler<typename PFP::REAL, EDGE>& edgeWeight,
	const AttributeHandler<typename PFP::REAL, VERTEX>& vertexArea,
	const AttributeHandler<ATTR_TYPE, VERTEX>& attr,
227
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
228
{
229
	FunctorLaplacianCotanRHS_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea, attr, coord) ;
230
	m.template foreach_orbit<VERTEX>(flc) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
231 232
}

233 234 235 236
/*******************************************************************************
 * END MATRIX DEFINITION
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
237 238 239 240 241 242
template <class SOLVER_TRAITS>
void endMatrix(LinearSolver<SOLVER_TRAITS>* s)
{
	s->end_system() ;
}

243 244 245 246
/*******************************************************************************
 * SOLVE SYSTEM
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
247 248 249 250 251 252
template <class SOLVER_TRAITS>
void solve(LinearSolver<SOLVER_TRAITS>* s)
{
	s->solve() ;
}

253 254 255 256
/*******************************************************************************
 * RESET SOLVER
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
257 258 259 260 261 262
template <class SOLVER_TRAITS>
void resetSolver(LinearSolver<SOLVER_TRAITS>* s, bool keepMatrix)
{
	s->reset(keepMatrix) ;
}

263 264 265 266
/*******************************************************************************
 * GET RESULTS
 *******************************************************************************/

Pierre Kraemer's avatar
Pierre Kraemer committed
267
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
268 269 270
void getResult(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
271 272
	const AttributeHandler<unsigned int, VERTEX> index,
	AttributeHandler<ATTR_TYPE, VERTEX>& attr)
Pierre Kraemer's avatar
Pierre Kraemer committed
273
{
274
	FunctorSolverToMesh_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> fstm(s, index, attr) ;
275
	m.template foreach_orbit<VERTEX>(fstm) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
276 277 278
}

template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
279 280 281
void getResult(
	typename PFP::MAP& m,
	LinearSolver<SOLVER_TRAITS>* s,
282 283
	const AttributeHandler<unsigned int, VERTEX> index,
	AttributeHandler<ATTR_TYPE, VERTEX>& attr,
284
	unsigned int coord)
Pierre Kraemer's avatar
Pierre Kraemer committed
285
{
286
	FunctorSolverToMesh_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> fstm(s, index, attr, coord) ;
287
	m.template foreach_orbit<VERTEX>(fstm) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
288 289 290 291 292 293 294
}

} // namespace LinearSolving

} // namespace CGoGN

#endif