Skip to content

Commit af98611

Browse files
alxbilgeralxbilger
authored
Support more types of linear systems (#542)
* Support more types of linear systems * add support for BTD
1 parent ae6b330 commit af98611

File tree

5 files changed

+262
-80
lines changed

5 files changed

+262
-80
lines changed

‎bindings/Modules/src/SofaPython3/SofaLinearSystem/Binding_LinearSystem.cpp‎

Lines changed: 88 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -17,15 +17,12 @@
1717
*******************************************************************************
1818
* Contact information: [email protected] *
1919
******************************************************************************/
20-
#include <SofaPython3/SofaLinearSystem/Binding_LinearSystem.h>
21-
#include <SofaPython3/SofaLinearSystem/Binding_LinearSystem_doc.h>
22-
#include <pybind11/pybind11.h>
20+
#include <SofaPython3/SofaLinearSystem/Binding_LinearSystem.inl>
2321

24-
#include <SofaPython3/Sofa/Core/Binding_Base.h>
25-
#include <SofaPython3/PythonFactory.h>
2622

27-
#include <sofa/component/linearsystem/TypedMatrixLinearSystem.h>
2823
#include <pybind11/eigen.h>
24+
#include <sofa/linearalgebra/BTDMatrix.h>
25+
#include <sofa/linearalgebra/BlockVector.h>
2926
#include <sofa/linearalgebra/CompressedRowSparseMatrix.h>
3027

3128

@@ -39,12 +36,6 @@ using EigenSparseMatrix = Eigen::SparseMatrix<Real, Eigen::RowMajor>;
3936
template<class Real>
4037
using EigenMatrixMap = Eigen::Map<Eigen::SparseMatrix<Real, Eigen::RowMajor> >;
4138

42-
template<class Real>
43-
using Vector = Eigen::Matrix<Real,Eigen::Dynamic, 1>;
44-
45-
template<class Real>
46-
using EigenVectorMap = Eigen::Map<Vector<Real>>;
47-
4839
template<class TBlock>
4940
EigenSparseMatrix<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>::Real>
5041
toEigen(sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>& matrix)
@@ -70,62 +61,104 @@ toEigen(sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>& matrix)
7061
}
7162
}
7263

73-
template<class TBlock>
74-
void bindLinearSystems(py::module &m)
64+
template<class TMatrix, class TVector>
65+
void bindMatrixAccessCRS(LinearSystemClass<TMatrix, TVector>& c)
7566
{
76-
using CRS = sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>;
77-
using Real = typename CRS::Real;
78-
using CRSLinearSystem = sofa::component::linearsystem::TypedMatrixLinearSystem<CRS, sofa::linearalgebra::FullVector<Real> >;
67+
using Real = typename TMatrix::Real;
68+
using CRSLinearSystem = sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>;
7969

80-
const std::string typeName = CRSLinearSystem::GetClass()->className + CRSLinearSystem::GetCustomTemplateName();
81-
82-
py::class_<CRSLinearSystem,
83-
sofa::core::objectmodel::BaseObject,
84-
sofapython3::py_shared_ptr<CRSLinearSystem> > c(m, typeName.c_str(), sofapython3::doc::linearsystem::linearSystemClass);
85-
86-
c.def("A", [](CRSLinearSystem& self) -> EigenSparseMatrix<Real>
87-
{
88-
if (CRS* matrix = self.getSystemMatrix())
70+
const auto matrixAccess =
71+
[](CRSLinearSystem& self) -> EigenSparseMatrix<Real>
8972
{
90-
if (matrix->colsValue.empty()) //null matrix: contains no entries
73+
if (TMatrix* matrix = self.getSystemMatrix())
9174
{
92-
return EigenSparseMatrix<Real>{matrix->rows(), matrix->cols()};
75+
if (matrix->colsValue.empty()) //null matrix: contains no entries
76+
{
77+
return EigenSparseMatrix<Real>{matrix->rows(), matrix->cols()};
78+
}
79+
return toEigen(*matrix);
9380
}
94-
return toEigen(*matrix);
95-
}
96-
return {};
97-
}, sofapython3::doc::linearsystem::linearSystem_A);
81+
return {};
82+
};
9883

99-
c.def("b", [](CRSLinearSystem& self) -> Vector<Real>
100-
{
101-
if (auto* vector = self.getRHSVector())
102-
{
103-
return EigenVectorMap<Real>(vector->ptr(), vector->size());
104-
}
105-
return {};
106-
}, sofapython3::doc::linearsystem::linearSystem_b);
84+
c.def("A", matrixAccess, sofapython3::doc::linearsystem::linearSystem_CRS_A);
85+
c.def("get_system_matrix", matrixAccess, sofapython3::doc::linearsystem::linearSystem_CRS_get_system_matrix);
86+
}
10787

108-
c.def("x", [](CRSLinearSystem& self) -> Vector<Real>
109-
{
110-
if (auto* vector = self.getSolutionVector())
111-
{
112-
return EigenVectorMap<Real>(vector->ptr(), vector->size());
113-
}
114-
return {};
115-
}, sofapython3::doc::linearsystem::linearSystem_x);
88+
template<>
89+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<SReal>, sofa::linearalgebra::FullVector<SReal>>& c)
90+
{
91+
bindMatrixAccessCRS(c);
92+
}
11693

94+
template<>
95+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<2,2,SReal>>, sofa::linearalgebra::FullVector<SReal>>& c)
96+
{
97+
bindMatrixAccessCRS(c);
98+
}
11799

118-
/// register the binding in the downcasting subsystem
119-
PythonFactory::registerType<CRSLinearSystem>([](sofa::core::objectmodel::Base* object)
120-
{
121-
return py::cast(dynamic_cast<CRSLinearSystem*>(object));
122-
});
100+
template<>
101+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<3,3,SReal>>, sofa::linearalgebra::FullVector<SReal>>& c)
102+
{
103+
bindMatrixAccessCRS(c);
104+
}
105+
106+
template<>
107+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<4,4,SReal>>, sofa::linearalgebra::FullVector<SReal>>& c)
108+
{
109+
bindMatrixAccessCRS(c);
110+
}
111+
112+
template<>
113+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<6,6,SReal>>, sofa::linearalgebra::FullVector<SReal>>& c)
114+
{
115+
bindMatrixAccessCRS(c);
116+
}
117+
118+
template<>
119+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<8,8,SReal>>, sofa::linearalgebra::FullVector<SReal>>& c)
120+
{
121+
bindMatrixAccessCRS(c);
122+
}
123+
124+
template<class Real>
125+
using DenseMatrix = Eigen::Matrix<Real,Eigen::Dynamic, Eigen::Dynamic>;
126+
127+
template<class Real>
128+
using EigenDenseMatrixMap = Eigen::Map<DenseMatrix<Real>>;
129+
130+
template<>
131+
void bindMatrixAccess(LinearSystemClass<sofa::linearalgebra::FullMatrix<SReal>, sofa::linearalgebra::FullVector<SReal>>& c)
132+
{
133+
using CRSLinearSystem = sofa::component::linearsystem::TypedMatrixLinearSystem<sofa::linearalgebra::FullMatrix<SReal>, sofa::linearalgebra::FullVector<SReal>>;
134+
135+
const auto matrixAccess =
136+
[](CRSLinearSystem& self) -> EigenDenseMatrixMap<SReal>
137+
{
138+
sofa::linearalgebra::FullMatrix<SReal>* matrix = self.getSystemMatrix();
139+
const auto row = matrix ? matrix->rows() : 0;
140+
const auto col = matrix ? matrix->cols() : 0;
141+
return EigenDenseMatrixMap<SReal>(matrix ? matrix->ptr() : nullptr, row, col);
142+
};
143+
c.def("A", matrixAccess, sofapython3::doc::linearsystem::linearSystem_FullMatrix_A);
144+
c.def("get_system_matrix", matrixAccess, sofapython3::doc::linearsystem::linearSystem_FullMatrix_get_system_matrix);
123145
}
124146

125147
void moduleAddLinearSystem(py::module &m)
126148
{
127-
bindLinearSystems<SReal>(m);
128-
bindLinearSystems<sofa::type::Mat<3, 3, SReal> >(m);
149+
bindLinearSystems<sofa::linearalgebra::FullMatrix<SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
150+
bindLinearSystems<sofa::linearalgebra::SparseMatrix<SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
151+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
152+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<2,2,SReal> >, sofa::linearalgebra::FullVector<SReal> >(m);
153+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<3,3,SReal> >, sofa::linearalgebra::FullVector<SReal> >(m);
154+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<4,4,SReal> >, sofa::linearalgebra::FullVector<SReal> >(m);
155+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<6,6,SReal> >, sofa::linearalgebra::FullVector<SReal> >(m);
156+
bindLinearSystems<sofa::linearalgebra::CompressedRowSparseMatrix<sofa::type::Mat<8,8,SReal> >, sofa::linearalgebra::FullVector<SReal> >(m);
157+
bindLinearSystems<sofa::linearalgebra::DiagonalMatrix<SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
158+
bindLinearSystems<sofa::linearalgebra::BlockDiagonalMatrix<3,SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
159+
bindLinearSystems<sofa::linearalgebra::RotationMatrix<SReal>, sofa::linearalgebra::FullVector<SReal> >(m);
160+
161+
bindLinearSystems<sofa::linearalgebra::BTDMatrix<6, SReal>, sofa::linearalgebra::BlockVector<6, SReal> >(m);
129162
}
130163

131164
}

‎bindings/Modules/src/SofaPython3/SofaLinearSystem/Binding_LinearSystem.inl‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
/******************************************************************************
2+
* SofaPython3 plugin *
3+
* (c) 2021 CNRS, University of Lille, INRIA *
4+
* *
5+
* This program is free software; you can redistribute it and/or modify it *
6+
* under the terms of the GNU Lesser General Public License as published by *
7+
* the Free Software Foundation; either version 2.1 of the License, or (at *
8+
* your option) any later version. *
9+
* *
10+
* This program is distributed in the hope that it will be useful, but WITHOUT *
11+
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
12+
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
13+
* for more details. *
14+
* *
15+
* You should have received a copy of the GNU Lesser General Public License *
16+
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
17+
*******************************************************************************
18+
* Contact information: [email protected] *
19+
******************************************************************************/
20+
#pragma once
21+
#include <SofaPython3/SofaLinearSystem/Binding_LinearSystem.h>
22+
#include <SofaPython3/SofaLinearSystem/Binding_LinearSystem_doc.h>
23+
#include <sofa/component/linearsystem/TypedMatrixLinearSystem.h>
24+
#include <pybind11/pybind11.h>
25+
26+
#include <SofaPython3/PythonFactory.h>
27+
#include <SofaPython3/Sofa/Core/Binding_Base.h>
28+
#include <sofa/linearalgebra/EigenSparseMatrix.h>
29+
30+
namespace py { using namespace pybind11; }
31+
32+
namespace sofapython3
33+
{
34+
35+
template<class Real>
36+
using Vector = Eigen::Matrix<Real,Eigen::Dynamic, 1>;
37+
38+
template<class Real>
39+
using EigenVectorMap = Eigen::Map<Vector<Real>>;
40+
41+
template<class TVector>
42+
Vector<typename TVector::Real>
43+
getVector(TVector* vector)
44+
{
45+
using Real = typename TVector::Real;
46+
if (vector)
47+
{
48+
return EigenVectorMap<Real>(vector->ptr(), vector->size());
49+
}
50+
return {};
51+
}
52+
53+
template<class TMatrix, class TVector>
54+
Vector<typename TVector::Real> getRHSVector(sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>& linearSystem)
55+
{
56+
return getVector(linearSystem.getRHSVector());
57+
}
58+
59+
template<class TMatrix, class TVector>
60+
Vector<typename TVector::Real> getSolutionVector(sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>& linearSystem)
61+
{
62+
return getVector(linearSystem.getSolutionVector());
63+
}
64+
65+
template<class TMatrix, class TVector>
66+
using LinearSystemClass =
67+
py::class_<sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>,
68+
sofa::core::objectmodel::BaseObject,
69+
sofapython3::py_shared_ptr<sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>> >;
70+
71+
template<class TMatrix, class TVector>
72+
void bindMatrixAccess(LinearSystemClass<TMatrix, TVector>& c)
73+
{}
74+
75+
template<class TMatrix, class TVector>
76+
void bindLinearSystems(py::module &m)
77+
{
78+
using LinearSystem = sofa::component::linearsystem::TypedMatrixLinearSystem<TMatrix, TVector>;
79+
80+
const std::string typeName =
81+
LinearSystem::GetClass()->className
82+
+ LinearSystem::GetCustomTemplateName();
83+
84+
LinearSystemClass<TMatrix, TVector> c(m, typeName.c_str(), sofapython3::doc::linearsystem::linearSystemClass);
85+
86+
c.def("b", getRHSVector<TMatrix, TVector>, sofapython3::doc::linearsystem::linearSystem_b);
87+
c.def("get_rhs_vector", getRHSVector<TMatrix, TVector>, sofapython3::doc::linearsystem::linearSystem_b);
88+
89+
c.def("x", getSolutionVector<TMatrix, TVector>, sofapython3::doc::linearsystem::linearSystem_x);
90+
c.def("get_solution_vector", getSolutionVector<TMatrix, TVector>, sofapython3::doc::linearsystem::linearSystem_x);
91+
92+
bindMatrixAccess<TMatrix, TVector>(c);
93+
94+
/// register the binding in the downcasting subsystem
95+
PythonFactory::registerType<LinearSystem>([](sofa::core::objectmodel::Base* object)
96+
{
97+
return py::cast(dynamic_cast<LinearSystem*>(object));
98+
});
99+
}
100+
101+
}

‎bindings/Modules/src/SofaPython3/SofaLinearSystem/Binding_LinearSystem_doc.h‎

Lines changed: 31 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ Linear system. Supports only CompressedRowSparseMatrix.
3232
linear_system = root.addObject('MatrixLinearSystem', template='CompressedRowSparseMatrixd')
3333
)";
3434

35-
static auto linearSystem_A =
35+
static auto linearSystem_CRS_A =
3636
R"(
3737
Returns the global system matrix as a scipy sparse matrix
3838
@@ -42,6 +42,36 @@ linear_system = root.addObject('MatrixLinearSystem', template='CompressedRowSpar
4242
matrix = linear_system.A()
4343
)";
4444

45+
static auto linearSystem_CRS_get_system_matrix =
46+
R"(
47+
Returns the global system matrix as a scipy sparse matrix
48+
49+
example:
50+
------------
51+
linear_system = root.addObject('MatrixLinearSystem', template='CompressedRowSparseMatrixd')
52+
matrix = linear_system.get_system_matrix()
53+
)";
54+
55+
static auto linearSystem_FullMatrix_A =
56+
R"(
57+
Returns the global system matrix as a numpy array matrix
58+
59+
example:
60+
------------
61+
linear_system = root.addObject('MatrixLinearSystem', template='FullMatrix')
62+
matrix = linear_system.A()
63+
)";
64+
65+
static auto linearSystem_FullMatrix_get_system_matrix =
66+
R"(
67+
Returns the global system matrix as a numpy array matrix
68+
69+
example:
70+
------------
71+
linear_system = root.addObject('MatrixLinearSystem', template='FullMatrix')
72+
matrix = linear_system.get_system_matrix()
73+
)";
74+
4575
static auto linearSystem_b =
4676
R"(
4777
Returns the global system right hand side as a numpy array

‎bindings/Modules/src/SofaPython3/SofaLinearSystem/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ project(Bindings.Modules.SofaLinearSystem)
22

33
set(HEADER_FILES
44
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSystem.h
5+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSystem.inl
56
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSystem_doc.h
67
)
78

0 commit comments

Comments
 (0)