diff --git a/SCHNApps/CMakeLists.txt b/SCHNApps/CMakeLists.txt index 5d839bec66d6355ac5e78ae43f3940253f423836..d9c1d83760719890053d17cc8da03fae2de8ea74 100644 --- a/SCHNApps/CMakeLists.txt +++ b/SCHNApps/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 2.8) project(SCHNApps) -SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -fPIC") +SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -fPIC -std=c++11") SET(CGoGN_ROOT_DIR ${CMAKE_SOURCE_DIR}/.. CACHE STRING "CGoGN root dir") diff --git a/include/Algo/LinearSolving/basic.h b/include/Algo/LinearSolving/basic.h index 209796c455c9098f31c12c71bd52f9a1cd0634a7..472b7e4870f69be31a7915cd1407779226c4c38d 100644 --- a/include/Algo/LinearSolving/basic.h +++ b/include/Algo/LinearSolving/basic.h @@ -26,9 +26,7 @@ #define __LINEAR_SOLVING_BASIC__ #include "NL/nl.h" -#include "Algo/LinearSolving/variablesSetup.h" -#include "Algo/LinearSolving/matrixSetup.h" -#include "Algo/Topo/basic.h" +#include "Topology/generic/traversor/traversorCell.h" namespace CGoGN { @@ -36,6 +34,15 @@ namespace CGoGN namespace LinearSolving { +template +struct Coeff +{ + unsigned int index; + CoeffType value; + Coeff(unsigned int i, CoeffType v) : index(i), value(v) + {} +} ; + /******************************************************************************* * VARIABLES SETUP *******************************************************************************/ @@ -44,29 +51,31 @@ template void setupVariables( typename PFP::MAP& m, const VertexAttribute& index, - const CellMarker& fm, + const CellMarker& freeMarker, const VertexAttribute& attr) { -// TraversorV t(m); -// for (Dart d = t.begin(); d != t.end(); d = t.next()) -// { - -// } - - FunctorMeshToSolver_Scalar fmts(index, fm, attr) ; - Algo::Topo::foreach_orbit(m, fmts) ; + foreach_cell(m, [&] (Dart d) + { + nlSetVariable(index[d], attr[d]); + if(!freeMarker.isMarked(d)) + nlLockVariable(index[d]); + }); } template void setupVariables( typename PFP::MAP& m, const VertexAttribute& index, - const CellMarker& fm, + const CellMarker& freeMarker, const VertexAttribute& attr, unsigned int coord) { - FunctorMeshToSolver_Vector fmts(index, fm, attr, coord) ; - Algo::Topo::foreach_orbit(m, fmts) ; + foreach_cell(m, [&] (Dart d) + { + nlSetVariable(index[d], (attr[d])[coord]); + if(!freeMarker.isMarked(d)) + nlLockVariable(index[d]); + }); } /******************************************************************************* @@ -81,8 +90,16 @@ void addRowsRHS_Equality( const VertexAttribute& weight) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorEquality_PerVertexWeight_Scalar feq(index, attr, weight) ; - Algo::Topo::foreach_orbit(m, feq) ; + + foreach_cell(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -94,8 +111,16 @@ void addRowsRHS_Equality( float weight) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorEquality_UniformWeight_Scalar feq(index, attr, weight) ; - Algo::Topo::foreach_orbit(m, feq) ; + + foreach_cell(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -108,8 +133,16 @@ void addRowsRHS_Equality( unsigned int coord) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorEquality_PerVertexWeight_Vector feq(index, attr, weight, coord) ; - Algo::Topo::foreach_orbit(m, feq) ; + + foreach_cell(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -122,8 +155,16 @@ void addRowsRHS_Equality( unsigned int coord) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorEquality_UniformWeight_Vector feq(index, attr, weight, coord) ; - Algo::Topo::foreach_orbit(m, feq) ; + + foreach_cell(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -134,36 +175,99 @@ void addRowsRHS_Equality( template void addRows_Laplacian_Topo( typename PFP::MAP& m, - const VertexAttribute index) + const VertexAttribute& index) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianTopo flt(m, index) ; - Algo::Topo::foreach_orbit(m, flt) ; + + foreach_cell(m, [&] (Dart d) + { + nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; + nlBegin(NL_ROW); + typename PFP::REAL aii = 0 ; + Traversor2VE 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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } template void addRowsRHS_Laplacian_Topo( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, const VertexAttribute& attr) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianTopoRHS_Scalar flt(m, index, attr) ; - Algo::Topo::foreach_orbit(m, flt) ; + + foreach_cell(m, [&] (Dart d) + { + std::vector > coeffs ; + coeffs.reserve(12) ; + + typename PFP::REAL norm2 = 0 ; + typename PFP::REAL aii = 0 ; + Traversor2VE 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(index[m.phi1(it)], aij)) ; + norm2 += aij * aij ; + } + coeffs.push_back(Coeff(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } template void addRowsRHS_Laplacian_Topo( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, const VertexAttribute& attr, unsigned int coord) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianTopoRHS_Vector flt(m, index, attr, coord) ; - Algo::Topo::foreach_orbit(m, flt) ; + + foreach_cell(m, [&] (Dart d) + { + std::vector > coeffs ; + coeffs.reserve(12) ; + + typename PFP::REAL norm2 = 0 ; + typename PFP::REAL aii = 0 ; + Traversor2VE 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(index[m.phi1(it)], aij)) ; + norm2 += aij * aij ; + } + coeffs.push_back(Coeff(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -174,57 +278,108 @@ void addRowsRHS_Laplacian_Topo( template void addRows_Laplacian_Cotan( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, const EdgeAttribute& edgeWeight, const VertexAttribute& vertexArea) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianCotan flc(m, index, edgeWeight, vertexArea) ; - Algo::Topo::foreach_orbit(m, flc) ; + + foreach_cell(m, [&] (Dart d) + { + nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; + nlBegin(NL_ROW); + typename PFP::REAL vArea = vertexArea[d] ; + typename PFP::REAL aii = 0 ; + Traversor2VE 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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } template void addRowsRHS_Laplacian_Cotan( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, const EdgeAttribute& edgeWeight, const VertexAttribute& vertexArea, const VertexAttribute& attr) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianCotanRHS_Scalar flc(m, index, edgeWeight, vertexArea, attr) ; - Algo::Topo::foreach_orbit(m, flc) ; + + foreach_cell(m, [&] (Dart d) + { + std::vector > coeffs ; + coeffs.reserve(12) ; + + typename PFP::REAL vArea = vertexArea[d] ; + typename PFP::REAL norm2 = 0 ; + typename PFP::REAL aii = 0 ; + Traversor2VE 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(index[m.phi1(it)], aij)) ; + norm2 += aij * aij ; + } + coeffs.push_back(Coeff(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } template void addRowsRHS_Laplacian_Cotan( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, const EdgeAttribute& edgeWeight, const VertexAttribute& vertexArea, const VertexAttribute& attr, unsigned int coord) { nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianCotanRHS_Vector flc(m, index, edgeWeight, vertexArea, attr, coord) ; - Algo::Topo::foreach_orbit(m, flc) ; - nlDisable(NL_NORMALIZE_ROWS) ; -} -template -void addRowsRHS_Laplacian_Cotan_NL( - typename PFP::MAP& m, - const VertexAttribute index, - const EdgeAttribute& edgeWeight, - const VertexAttribute& vertexArea, - const VertexAttribute& attr, - unsigned int coord) -{ - nlEnable(NL_NORMALIZE_ROWS) ; - FunctorLaplacianCotanRHS_Vector flc(m, index, edgeWeight, vertexArea, attr, coord) ; - Algo::Topo::foreach_orbit(m, flc) ; + foreach_cell(m, [&] (Dart d) + { + std::vector > coeffs ; + coeffs.reserve(12) ; + + typename PFP::REAL vArea = vertexArea[d] ; + typename PFP::REAL norm2 = 0 ; + typename PFP::REAL aii = 0 ; + Traversor2VE 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(index[m.phi1(it)], aij)) ; + norm2 += aij * aij ; + } + coeffs.push_back(Coeff(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) ; + }); + nlDisable(NL_NORMALIZE_ROWS) ; } @@ -235,22 +390,26 @@ void addRowsRHS_Laplacian_Cotan_NL( template void getResult( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, VertexAttribute& attr) { - FunctorSolverToMesh_Scalar fstm(index, attr) ; - Algo::Topo::foreach_orbit(m, fstm) ; + foreach_cell(m, [&] (Dart d) + { + attr[d] = nlGetVariable(index[d]) ; + }); } template void getResult( typename PFP::MAP& m, - const VertexAttribute index, + const VertexAttribute& index, VertexAttribute& attr, unsigned int coord) { - FunctorSolverToMesh_Vector fstm(index, attr, coord) ; - Algo::Topo::foreach_orbit(m, fstm) ; + foreach_cell(m, [&] (Dart d) + { + (attr[d])[coord] = nlGetVariable(index[d]) ; + }); } } // namespace LinearSolving diff --git a/include/Algo/LinearSolving/matrixSetup.h b/include/Algo/LinearSolving/matrixSetup.h deleted file mode 100644 index 6fd6e1715576ad9894c3eed27ae61845e5397026..0000000000000000000000000000000000000000 --- a/include/Algo/LinearSolving/matrixSetup.h +++ /dev/null @@ -1,477 +0,0 @@ -/******************************************************************************* -* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps * -* version 0.1 * -* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg * -* * -* 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. * -* * -* Web site: http://cgogn.unistra.fr/ * -* Contact information: cgogn@unistra.fr * -* * -*******************************************************************************/ - -#ifndef __LINEAR_SOLVING_MATRIX_SETUP__ -#define __LINEAR_SOLVING_MATRIX_SETUP__ - -namespace CGoGN -{ - -namespace LinearSolving -{ - -template -struct Coeff -{ - unsigned int index; - CoeffType value; - Coeff(unsigned int i, CoeffType v) : index(i), value(v) - {} -} ; - -/******************************************************************************* - * EQUALITY MATRIX : right-hand-side SCALAR - *******************************************************************************/ - -template -class FunctorEquality_PerVertexWeight_Scalar : public FunctorType -{ - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - const VertexAttribute& weightTable ; - -public: - FunctorEquality_PerVertexWeight_Scalar( - const VertexAttribute& index, - const VertexAttribute& attr, - const VertexAttribute& weight - ) : indexTable(index), attrTable(attr), weightTable(weight) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[d]) ; - nlRowParameterd(NL_ROW_SCALING, weightTable[d]) ; - nlBegin(NL_ROW) ; - nlCoefficient(indexTable[d], 1) ; - nlEnd(NL_ROW) ; - return false ; - } -} ; - -template -class FunctorEquality_UniformWeight_Scalar : public FunctorType -{ - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - float weight ; - -public: - FunctorEquality_UniformWeight_Scalar( - const VertexAttribute& index, - const VertexAttribute& attr, - float w - ) : indexTable(index), attrTable(attr), weight(w) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, attrTable[d]) ; - nlRowParameterd(NL_ROW_SCALING, weight) ; - nlBegin(NL_ROW) ; - nlCoefficient(indexTable[d], 1) ; - nlEnd(NL_ROW) ; - return false ; - } -} ; - -/******************************************************************************* - * EQUALITY MATRIX : right-hand-side VECTOR + coordinate - *******************************************************************************/ - -template -class FunctorEquality_PerVertexWeight_Vector : public FunctorType -{ - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - const VertexAttribute& weightTable ; - unsigned int coord ; - -public: - FunctorEquality_PerVertexWeight_Vector( - const VertexAttribute& index, - const VertexAttribute& attr, - const VertexAttribute& weight, - unsigned int c - ) : indexTable(index), attrTable(attr), weightTable(weight), coord(c) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[d])[coord]) ; - nlRowParameterd(NL_ROW_SCALING, weightTable[d]) ; - nlBegin(NL_ROW) ; - nlCoefficient(indexTable[d], 1) ; - nlEnd(NL_ROW) ; - return false ; - } -} ; - -template -class FunctorEquality_UniformWeight_Vector : public FunctorType -{ - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - float weight ; - unsigned int coord ; - -public: - FunctorEquality_UniformWeight_Vector( - const VertexAttribute& index, - const VertexAttribute& attr, - float w, - unsigned int c - ) : indexTable(index), attrTable(attr), weight(w), coord(c) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, (attrTable[d])[coord]) ; - nlRowParameterd(NL_ROW_SCALING, weight) ; - nlBegin(NL_ROW) ; - nlCoefficient(indexTable[d], 1) ; - nlEnd(NL_ROW) ; - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN TOPO MATRIX : right-hand-side 0 - *******************************************************************************/ - -template -class FunctorLaplacianTopo : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - -public: - FunctorLaplacianTopo( - MAP& m, - const VertexAttribute& index - ) : FunctorMap(m), indexTable(index) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; - nlBegin(NL_ROW); - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = 1 ; - aii += aij ; - nlCoefficient(indexTable[this->m_map.phi1(it)], aij) ; - } - nlCoefficient(indexTable[d], -aii) ; - nlEnd(NL_ROW) ; - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN TOPO MATRIX : right-hand-side SCALAR - *******************************************************************************/ - -template -class FunctorLaplacianTopoRHS_Scalar : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - -public: - FunctorLaplacianTopoRHS_Scalar( - MAP& m, - const VertexAttribute& index, - const VertexAttribute& attr - ) : FunctorMap(m), indexTable(index), attrTable(attr) - {} - - bool operator()(Dart d) - { - std::vector > coeffs ; - coeffs.reserve(12) ; - - REAL norm2 = 0 ; - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = 1 ; - aii += aij ; - coeffs.push_back(Coeff(indexTable[this->m_map.phi1(it)], aij)) ; - norm2 += aij * aij ; - } - coeffs.push_back(Coeff(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) ; - - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN TOPO MATRIX : right-hand-side VECTOR + coordinate - *******************************************************************************/ - -template -class FunctorLaplacianTopoRHS_Vector : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - const VertexAttribute& attrTable ; - unsigned int coord ; - -public: - FunctorLaplacianTopoRHS_Vector( - MAP& m, - const VertexAttribute& index, - const VertexAttribute& attr, - unsigned int c - ) : FunctorMap(m), indexTable(index), attrTable(attr), coord(c) - {} - - bool operator()(Dart d) - { - std::vector > coeffs ; - coeffs.reserve(12) ; - - REAL norm2 = 0 ; - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = 1 ; - aii += aij ; - coeffs.push_back(Coeff(indexTable[this->m_map.phi1(it)], aij)) ; - norm2 += aij * aij ; - } - coeffs.push_back(Coeff(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) ; - - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN COTAN MATRIX : right-hand-side 0 - *******************************************************************************/ - -template -class FunctorLaplacianCotan : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - const EdgeAttribute& edgeWeight ; - const VertexAttribute& vertexArea ; - -public: - FunctorLaplacianCotan( - MAP& m, - const VertexAttribute& index, - const EdgeAttribute& eWeight, - const VertexAttribute& vArea - ) : FunctorMap(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea) - {} - - bool operator()(Dart d) - { - nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; - nlBegin(NL_ROW); - REAL vArea = vertexArea[d] ; - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = edgeWeight[it] / vArea ; - aii += aij ; - nlCoefficient(indexTable[this->m_map.phi1(it)], aij) ; - } - nlCoefficient(indexTable[d], -aii) ; - nlEnd(NL_ROW) ; - - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN COTAN MATRIX : right-hand-side SCALAR - *******************************************************************************/ - -template -class FunctorLaplacianCotanRHS_Scalar : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - const EdgeAttribute& edgeWeight ; - const VertexAttribute& vertexArea ; - const VertexAttribute& attrTable ; - -public: - FunctorLaplacianCotanRHS_Scalar( - MAP& m, - const VertexAttribute& index, - const EdgeAttribute& eWeight, - const VertexAttribute& vArea, - const VertexAttribute& attr - ) : FunctorMap(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr) - {} - - bool operator()(Dart d) - { - std::vector > coeffs ; - coeffs.reserve(12) ; - - REAL vArea = vertexArea[d] ; - REAL norm2 = 0 ; - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = edgeWeight[it] / vArea ; - aii += aij ; - coeffs.push_back(Coeff(indexTable[this->m_map.phi1(it)], aij)) ; - norm2 += aij * aij ; - } - coeffs.push_back(Coeff(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) ; - - return false ; - } -} ; - -/******************************************************************************* - * LAPLACIAN COTAN MATRIX : right-hand-side VECTOR + coordinate - *******************************************************************************/ - -template -class FunctorLaplacianCotanRHS_Vector : public FunctorMap -{ - typedef typename PFP::MAP MAP ; - typedef typename PFP::MAP::IMPL MAP_IMPL; - typedef typename PFP::REAL REAL ; - -protected: - const VertexAttribute& indexTable ; - const EdgeAttribute& edgeWeight ; - const VertexAttribute& vertexArea ; - const VertexAttribute& attrTable ; - unsigned int coord ; - -public: - FunctorLaplacianCotanRHS_Vector( - MAP& m, - const VertexAttribute& index, - const EdgeAttribute& eWeight, - const VertexAttribute& vArea, - const VertexAttribute& attr, - unsigned int c - ) : FunctorMap(m), indexTable(index), edgeWeight(eWeight), vertexArea(vArea), attrTable(attr), coord(c) - {} - - bool operator()(Dart d) - { - std::vector > coeffs ; - coeffs.reserve(12) ; - - REAL vArea = vertexArea[d] ; - REAL norm2 = 0 ; - REAL aii = 0 ; - Traversor2VE t(this->m_map, d) ; - for(Dart it = t.begin(); it != t.end(); it = t.next()) - { - REAL aij = edgeWeight[it] / vArea ; - aii += aij ; - coeffs.push_back(Coeff(indexTable[this->m_map.phi1(it)], aij)) ; - norm2 += aij * aij ; - } - coeffs.push_back(Coeff(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) ; - - return false ; - } -} ; - -} // namespace LinearSolving - -} // namespace CGoGN - -#endif diff --git a/include/Algo/LinearSolving/variablesSetup.h b/include/Algo/LinearSolving/variablesSetup.h deleted file mode 100644 index e5956c5184859ea078526527bfca73adffedc2a7..0000000000000000000000000000000000000000 --- a/include/Algo/LinearSolving/variablesSetup.h +++ /dev/null @@ -1,158 +0,0 @@ -/******************************************************************************* -* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps * -* version 0.1 * -* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg * -* * -* 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. * -* * -* Web site: http://cgogn.unistra.fr/ * -* Contact information: cgogn@unistra.fr * -* * -*******************************************************************************/ - -#ifndef __LINEAR_SOLVING_VARIABLES_SETUP__ -#define __LINEAR_SOLVING_VARIABLES_SETUP__ - -namespace CGoGN -{ - -namespace LinearSolving -{ - -template -class FunctorMeshToSolver_Scalar : public FunctorType -{ - typedef typename PFP::MAP MAP; - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const CellMarker& freeMarker ; - const VertexAttribute& attrTable ; - bool lockedVertices ; - -public: - FunctorMeshToSolver_Scalar( - const VertexAttribute& index, - const CellMarker& fm, - const VertexAttribute& attr - ) : indexTable(index), freeMarker(fm), attrTable(attr), lockedVertices(false) - {} - - bool operator()(Dart d) - { - nlSetVariable(indexTable[d], attrTable[d]); - if(!freeMarker.isMarked(d)) - { - nlLockVariable(indexTable[d]); - lockedVertices = true ; - } - return false ; - } - - bool hasLockedVertices() { return lockedVertices ; } -} ; - -template -class FunctorMeshToSolver_Vector : public FunctorType -{ - typedef typename PFP::MAP MAP; - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - const CellMarker& freeMarker ; - const VertexAttribute& attrTable ; - unsigned int coord ; - bool lockedVertices ; - -public: - FunctorMeshToSolver_Vector( - const VertexAttribute& index, - const CellMarker& fm, - const VertexAttribute& attr, - unsigned int c - ) : indexTable(index), freeMarker(fm), attrTable(attr), coord(c), lockedVertices(false) - {} - - bool operator()(Dart d) - { - nlSetVariable(indexTable[d], (attrTable[d])[coord]); - if(!freeMarker.isMarked(d)) - { - nlLockVariable(indexTable[d]); - lockedVertices = true ; - } - return false ; - } - - bool hasLockedVertices() { return lockedVertices ; } -} ; - -template -class FunctorSolverToMesh_Scalar : public FunctorType -{ - typedef typename PFP::MAP MAP; - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - VertexAttribute& attrTable ; - -public: - FunctorSolverToMesh_Scalar( - const VertexAttribute& index, - VertexAttribute& attr - ) : indexTable(index), attrTable(attr) - {} - - bool operator()(Dart d) - { - attrTable[d] = nlGetVariable(indexTable[d]) ; - return false ; - } -} ; - -template -class FunctorSolverToMesh_Vector : public FunctorType -{ - typedef typename PFP::MAP MAP; - typedef typename PFP::MAP::IMPL MAP_IMPL; - -protected: - const VertexAttribute& indexTable ; - VertexAttribute& attrTable ; - unsigned int coord ; - -public: - FunctorSolverToMesh_Vector( - const VertexAttribute& index, - VertexAttribute& attr, - unsigned int c - ) : indexTable(index), attrTable(attr), coord(c) - {} - - bool operator()(Dart d) - { - (attrTable[d])[coord] = nlGetVariable(indexTable[d]) ; - return false ; - } -} ; - -} // namespace LinearSolving - -} // namespace CGoGN - -#endif