Skip to content

Commit 64ecd52

Browse files
damienmarchaldamienmarchal
authored
Merge pull request sofa-framework#260 from alxbilger/matrixaccess
Add bindings for matrix access
2 parents 74feb01 + 1aaf67b commit 64ecd52

File tree

12 files changed

+393
-2
lines changed

12 files changed

+393
-2
lines changed

‎bindings/BindingsConfig.cmake.in‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ find_package(SofaPython3 QUIET REQUIRED COMPONENTS
1111
Bindings.SofaGui
1212
Bindings.SofaRuntime
1313
Bindings.SofaTypes
14+
Bindings.SofaLinearSolver
1415
)
1516
if(SP3_WITH_SOFAEXPORTER)
1617
find_package(SofaPython3 QUIET REQUIRED COMPONENTS Bindings.SofaExporter)

‎bindings/Modules/Bindings.ModulesConfig.cmake.in‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@ find_package(SofaPython3 QUIET REQUIRED COMPONENTS Plugin Bindings.Sofa)
99
find_package(SofaBaseTopology QUIET REQUIRED)
1010
# Required by Bindings.Modules.SofaDeformable
1111
find_package(SofaDeformable QUIET REQUIRED)
12+
# Required by Bindings.Modules.SofaLinearSolver
13+
find_package(Sofa.Component.LinearSolver.Iterative REQUIRED)
1214

1315
# If we are importing this config file and the target is not yet there this is indicating that
1416
# target is an imported one. So we include it

‎bindings/Modules/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
project(Bindings.Modules)
22

3-
set(MODULEBINDINGS_MODULE_LIST SofaBaseTopology SofaDeformable)
3+
set(MODULEBINDINGS_MODULE_LIST SofaBaseTopology SofaDeformable SofaLinearSolver)
44

55
find_package(Sofa.GL QUIET)
66
if(Sofa.GL_FOUND)

‎bindings/Modules/src/SofaPython3/SofaLinearSolver/Binding_LinearSolver.cpp‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
/******************************************************************************
2+
* SOFA, Simulation Open-Framework Architecture *
3+
* (c) 2021 INRIA, USTL, UJF, CNRS, MGH *
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+
#include <SofaPython3/SofaLinearSolver/Binding_LinearSolver.h>
21+
#include <SofaPython3/SofaLinearSolver/Binding_LinearSolver_doc.h>
22+
#include <pybind11/pybind11.h>
23+
24+
#include <SofaPython3/Sofa/Core/Binding_Base.h>
25+
#include <SofaPython3/PythonFactory.h>
26+
27+
#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
28+
#include <pybind11/eigen.h>
29+
30+
namespace py { using namespace pybind11; }
31+
32+
namespace sofapython3 {
33+
34+
using EigenSparseMatrix = Eigen::SparseMatrix<SReal, Eigen::RowMajor>;
35+
using EigenMatrixMap = Eigen::Map<Eigen::SparseMatrix<SReal, Eigen::RowMajor> >;
36+
using Vector = Eigen::Matrix<SReal,Eigen::Dynamic, 1>;
37+
using EigenVectorMap = Eigen::Map<Vector>;
38+
39+
template<class TBlock>
40+
EigenSparseMatrix toEigen(sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>& matrix)
41+
{
42+
sofa::linearalgebra::CompressedRowSparseMatrix<typename TBlock::Real> filtered;
43+
filtered.copyNonZeros(matrix);
44+
filtered.compress();
45+
return EigenMatrixMap(filtered.rows(), filtered.cols(), filtered.getColsValue().size(),
46+
(EigenMatrixMap::StorageIndex*)filtered.rowBegin.data(), (EigenMatrixMap::StorageIndex*)filtered.colsIndex.data(), filtered.colsValue.data());
47+
}
48+
49+
template<>
50+
EigenSparseMatrix toEigen<SReal>(sofa::linearalgebra::CompressedRowSparseMatrix<SReal>& matrix)
51+
{
52+
return EigenMatrixMap(matrix.rows(), matrix.cols(), matrix.getColsValue().size(),
53+
(EigenMatrixMap::StorageIndex*)matrix.rowBegin.data(), (EigenMatrixMap::StorageIndex*)matrix.colsIndex.data(), matrix.colsValue.data());
54+
}
55+
56+
template<class TBlock>
57+
void bindLinearSolvers(py::module &m)
58+
{
59+
using CRS = sofa::linearalgebra::CompressedRowSparseMatrix<TBlock>;
60+
using CRSLinearSolver = sofa::component::linearsolver::MatrixLinearSolver<CRS, sofa::linearalgebra::FullVector<SReal> >;
61+
using Real = typename CRS::Real;
62+
63+
const std::string typeName = CRSLinearSolver::GetClass()->className + CRSLinearSolver::GetCustomTemplateName();
64+
65+
py::class_<CRSLinearSolver,
66+
sofa::core::objectmodel::BaseObject,
67+
sofapython3::py_shared_ptr<CRSLinearSolver> > c(m, typeName.c_str(), sofapython3::doc::linearsolver::linearSolverClass);
68+
69+
c.def("A", [](CRSLinearSolver& self) -> EigenSparseMatrix
70+
{
71+
if (CRS* matrix = self.getSystemMatrix())
72+
{
73+
return toEigen(*matrix);
74+
}
75+
return {};
76+
}, sofapython3::doc::linearsolver::linearSolver_A);
77+
78+
c.def("b", [](CRSLinearSolver& self) -> Vector
79+
{
80+
if (auto* vector = self.getSystemRHVector())
81+
{
82+
return EigenVectorMap(vector->ptr(), vector->size());
83+
}
84+
return {};
85+
}, sofapython3::doc::linearsolver::linearSolver_b);
86+
87+
c.def("x", [](CRSLinearSolver& self) -> Vector
88+
{
89+
if (auto* vector = self.getSystemLHVector())
90+
{
91+
return EigenVectorMap(vector->ptr(), vector->size());
92+
}
93+
return {};
94+
}, sofapython3::doc::linearsolver::linearSolver_x);
95+
96+
97+
/// register the binding in the downcasting subsystem
98+
PythonFactory::registerType<CRSLinearSolver>([](sofa::core::objectmodel::Base* object)
99+
{
100+
return py::cast(dynamic_cast<CRSLinearSolver*>(object));
101+
});
102+
}
103+
104+
void moduleAddLinearSolver(py::module &m)
105+
{
106+
bindLinearSolvers<SReal>(m);
107+
bindLinearSolvers<sofa::type::Mat<3, 3, SReal> >(m);
108+
}
109+
110+
}

‎bindings/Modules/src/SofaPython3/SofaLinearSolver/Binding_LinearSolver.h‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
/******************************************************************************
2+
* SOFA, Simulation Open-Framework Architecture *
3+
* (c) 2021 INRIA, USTL, UJF, CNRS, MGH *
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+
21+
#pragma once
22+
23+
#include <pybind11/pybind11.h>
24+
25+
namespace sofapython3
26+
{
27+
28+
void moduleAddLinearSolver(pybind11::module &m);
29+
30+
} /// namespace sofapython3
31+

‎bindings/Modules/src/SofaPython3/SofaLinearSolver/Binding_LinearSolver_doc.h‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
/******************************************************************************
2+
* SOFA, Simulation Open-Framework Architecture *
3+
* (c) 2021 INRIA, USTL, UJF, CNRS, MGH *
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+
21+
#pragma once
22+
23+
namespace sofapython3::doc::linearsolver {
24+
25+
static auto linearSolverClass =
26+
R"(
27+
Linear solver. Supports only scalar CompressedRowSparseMatrix.
28+
29+
example:
30+
------------
31+
32+
linear_solver = root.addObject('SparseLDLSolver', template='CompressedRowSparseMatrixd') # supported
33+
linear_solver = root.addObject('SparseLDLSolver', template='CompressedRowSparseMatrixMat3x3d') # not supported
34+
)";
35+
36+
static auto linearSolver_A =
37+
R"(
38+
Returns the global system matrix as a scipy sparse matrix
39+
40+
example:
41+
------------
42+
linear_solver = root.addObject('SparseLDLSolver', template='CompressedRowSparseMatrixd')
43+
matrix = linear_solver.A()
44+
)";
45+
46+
static auto linearSolver_b =
47+
R"(
48+
Returns the global system right hand side as a numpy array
49+
50+
example:
51+
------------
52+
linear_solver = root.addObject('SparseLDLSolver', template='CompressedRowSparseMatrixd')
53+
matrix = linear_solver.b()
54+
)";
55+
56+
static auto linearSolver_x =
57+
R"(
58+
Returns the global system solution vector as a numpy array
59+
60+
example:
61+
------------
62+
linear_solver = root.addObject('SparseLDLSolver', template='CompressedRowSparseMatrixd')
63+
matrix = linear_solver.x()
64+
)";
65+
66+
} // namespace sofapython3::doc::baseCamera

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

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
project(Bindings.Modules.SofaLinearSolver)
2+
3+
set(HEADER_FILES
4+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSolver.h
5+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSolver_doc.h
6+
)
7+
8+
set(SOURCE_FILES
9+
${CMAKE_CURRENT_SOURCE_DIR}/Binding_LinearSolver.cpp
10+
${CMAKE_CURRENT_SOURCE_DIR}/Module_SofaLinearSolver.cpp
11+
)
12+
13+
if (NOT TARGET SofaPython3::Plugin)
14+
find_package(SofaPython3 REQUIRED COMPONENTS Plugin Bindings.Sofa)
15+
endif()
16+
17+
find_package(Sofa.Component.LinearSolver.Iterative REQUIRED)
18+
19+
SP3_add_python_module(
20+
TARGET ${PROJECT_NAME}
21+
PACKAGE Bindings.Modules
22+
MODULE SofaLinearSolver
23+
DESTINATION Sofa
24+
SOURCES ${SOURCE_FILES}
25+
HEADERS ${HEADER_FILES}
26+
DEPENDS Sofa.Component.LinearSolver.Iterative SofaPython3::Plugin SofaPython3::Bindings.Sofa.Core
27+
)

‎bindings/Modules/src/SofaPython3/SofaLinearSolver/Module_SofaLinearSolver.cpp‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
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+
21+
#include <pybind11/pybind11.h>
22+
#include <SofaPython3/SofaLinearSolver/Binding_LinearSolver.h>
23+
24+
namespace sofapython3
25+
{
26+
27+
PYBIND11_MODULE(SofaLinearSolver, m)
28+
{
29+
moduleAddLinearSolver(m);
30+
}
31+
32+
} // namespace sofapython3
33+

‎bindings/Modules/tests/CMakeLists.txt‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,8 @@ set(SOURCE_FILES
77
set(PYTHON_FILES
88
${CMAKE_CURRENT_SOURCE_DIR}/SofaDeformable/LinearSpring.py
99
${CMAKE_CURRENT_SOURCE_DIR}/SofaDeformable/SpringForceField.py
10+
11+
${CMAKE_CURRENT_SOURCE_DIR}/SofaLinearSolver/matrix_access.py
1012
)
1113

1214
find_package(SofaGTestMain REQUIRED)
@@ -26,7 +28,7 @@ sofa_auto_set_target_rpath(
2628

2729
add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME})
2830

29-
set(DIR_BINDING_LIST SofaBaseTopology SofaDeformable)
31+
set(DIR_BINDING_LIST SofaBaseTopology SofaDeformable SofaLinearSolver)
3032
get_property(_isMultiConfig GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG)
3133
foreach(dir_binding ${DIR_BINDING_LIST})
3234
if (NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${dir_binding}")

‎bindings/Modules/tests/SofaLinearSolver/matrix_access.py‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
import unittest
2+
import Sofa.Core
3+
import Sofa.Components
4+
from Sofa import SofaLinearSolver
5+
6+
class Test(unittest.TestCase):
7+
8+
def simulate_beam(self, linear_solver_template):
9+
root = Sofa.Core.Node("rootNode")
10+
11+
loop = root.addObject('DefaultAnimationLoop')
12+
13+
root.addObject('RequiredPlugin', name='Sofa.Component.ODESolver.Backward')
14+
root.addObject('RequiredPlugin', name='Sofa.Component.LinearSolver.Direct')
15+
root.addObject('RequiredPlugin', name='Sofa.Component.Engine.Select')
16+
root.addObject('RequiredPlugin', name='Sofa.Component.Constraint.Projective')
17+
root.addObject('RequiredPlugin', name='Sofa.Component.SolidMechanics.FEM.Elastic')
18+
root.addObject('RequiredPlugin', name='Sofa.Component.Mass')
19+
20+
root.addObject('EulerImplicitSolver', rayleighStiffness="0.1", rayleighMass="0.1")
21+
linear_solver = root.addObject('SparseLDLSolver', name='linear_solver', applyPermutation="false", template=linear_solver_template)
22+
23+
root.addObject('MechanicalObject', name="DoFs")
24+
root.addObject('UniformMass', name="mass", totalMass="320")
25+
root.addObject('RegularGridTopology', name="grid", nx="4", ny="4", nz="20", xmin="-9", xmax="-6", ymin="0", ymax="3", zmin="0", zmax="19")
26+
root.addObject('BoxROI', name="box", box=[-10, -1, -0.0001, -5, 4, 0.0001])
27+
root.addObject('FixedConstraint', indices="@box.indices")
28+
root.addObject('HexahedronFEMForceField', name="FEM", youngModulus="4000", poissonRatio="0.3", method="large")
29+
30+
Sofa.Simulation.init(root)
31+
Sofa.Simulation.animate(root, 0.0001)
32+
33+
return root
34+
35+
def test_matrix_access_blocks_scalar(self):
36+
root = self.simulate_beam("CompressedRowSparseMatrixd")
37+
A = root.linear_solver.A()
38+
39+
self.assertEqual(A.ndim, 2)
40+
self.assertEqual(A.shape, (960, 960))
41+
self.assertGreater(A.nnz, 960) // Should be Equal to one specific value but a different between OS has been reported here: https://github.com/sofa-framework/sofa/issues/3036
42+
43+
def test_matrix_access_blocks3x3(self):
44+
root = self.simulate_beam("CompressedRowSparseMatrixMat3x3d")
45+
A = root.linear_solver.A()
46+
47+
self.assertEqual(A.ndim, 2)
48+
self.assertEqual(A.shape, (960, 960))
49+
self.assertGreater(A.nnz, 960) // Should be Equal to one specific value but a different between OS has been reported here: https://github.com/sofa-framework/sofa/issues/3036

0 commit comments

Comments
 (0)