ODE solver in python

bindings/Modules/CMakeLists.txt

@@ -1,6 +1,6 @@
 project(Bindings.Modules)
 
-set(MODULEBINDINGS_MODULE_LIST SofaBaseTopology SofaDeformable)
+set(MODULEBINDINGS_MODULE_LIST SofaBaseTopology SofaDeformable SofaSimulationCore)
 
 find_package(Sofa.GL QUIET)
 if(Sofa.GL_FOUND)

bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_MechanicalOperations.cpp

@@ -0,0 +1,84 @@
+/******************************************************************************
+*                              SofaPython3 plugin                             *
+*                  (c) 2021 CNRS, University of Lille, INRIA                  *
+*                                                                             *
+* This program 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 program 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 program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Contact information: [email protected]                             *
+******************************************************************************/
+
+#include <pybind11/pybind11.h>
+
+#include <sofa/simulation/MechanicalOperations.h>
+
+namespace py { using namespace pybind11; }
+
+namespace sofapython3 {
+
+void moduleAddMechanicalOperations(pybind11::module& m)
+{
+    using namespace sofa::simulation::common;
+
+    py::class_<MechanicalOperations> c (m, "MechanicalOperations");
+    c.def(py::init<const sofa::core::MechanicalParams*, sofa::core::objectmodel::BaseContext*, bool>(), py::arg("mparams"), py::arg("ctx"), py::arg("precomputedTraversalOrder")=false);
+
+    // Mechanical Vector operations
+    c.def("propagateDx", &MechanicalOperations::propagateDx, py::arg("dx"), py::arg("ignore_flag") = false, "Propagate the given displacement through all mappings");
+    c.def("propagateDxAndResetDf", &MechanicalOperations::propagateDxAndResetDf, py::arg("dx"), py::arg("df"), "Propagate the given displacement through all mappings and reset the current force delta");
+    c.def("propagateX", &MechanicalOperations::propagateX, py::arg("x"), "Propagate the given position through all mappings");
+    c.def("propagateV", &MechanicalOperations::propagateV, py::arg("v"), "Propagate the given velocity through all mappings");
+    c.def("propagateXAndV", &MechanicalOperations::propagateXAndV, py::arg("x"), py::arg("v"), "Propagate the given position and velocity through all mappings");
+    c.def("propagateXAndResetF", &MechanicalOperations::propagateXAndResetF, py::arg("x"), py::arg("f"), "Propagate the given position through all mappings and reset the current force delta");
+    c.def("projectPosition", &MechanicalOperations::projectPosition, py::arg("x"), py::arg("time") = 0., "Apply projective constraints to the given position vector");
+    c.def("projectVelocity", &MechanicalOperations::projectVelocity, py::arg("v"), py::arg("time") = 0., "Apply projective constraints to the given velocity vector");
+    c.def("projectResponse", [](MechanicalOperations & self, sofa::core::MultiVecDerivId dx) {
+        self.projectResponse(dx);
+    }, py::arg("dx"), "Apply projective constraints to the given vector");
+    c.def("projectPositionAndVelocity", &MechanicalOperations::projectPositionAndVelocity, py::arg("x"), py::arg("v"), py::arg("time") = 0., "Apply projective constraints to the given position and velocity vectors");
+    c.def("addMdx", &MechanicalOperations::addMdx, py::arg("res"), py::arg("dx"), py::arg("factor") = 1.0, "res += factor M.dx");
+    c.def("integrateVelocity", &MechanicalOperations::integrateVelocity, py::arg("res"), py::arg("x"), py::arg("v"), py::arg("dt"), "res = x + v.dt");
+    c.def("accFromF", &MechanicalOperations::accFromF, py::arg("a"), py::arg("f"), "a = M^-1 . f");
+    c.def("computeEnergy", &MechanicalOperations::computeEnergy, py::arg("kineticEnergy"), py::arg("potentialEnergy"), "Compute Energy");
+    c.def("computeForce", py::overload_cast<sofa::core::MultiVecDerivId, bool, bool, bool>(&MechanicalOperations::computeForce), py::arg("result"), py::arg("clear")=true, py::arg("accumulate")=true, py::arg("neglectingCompliance")=true, "Compute the current force (given the latest propagated position and velocity)");
+    c.def("computeForce", py::overload_cast<SReal, sofa::core::MultiVecDerivId, sofa::core::MultiVecCoordId, sofa::core::MultiVecDerivId, bool>(&MechanicalOperations::computeForce), py::arg("t"), py::arg("f"), py::arg("x"), py::arg("v"), py::arg("neglectingCompliance")=true, "Compute f(x,v) at time t. Parameters x and v not const due to propagation through mappings.");
+    c.def("computeDf", &MechanicalOperations::computeDf, py::arg("df"), py::arg("clear")=true, py::arg("accumulate")=true, "Compute the current force delta (given the latest propagated displacement)");
+    c.def("computeDfV", &MechanicalOperations::computeDfV, py::arg("df"), py::arg("clear")=true, py::arg("accumulate")=true, "Compute the current force delta (given the latest propagated velocity)");
+    c.def("addMBKdx", &MechanicalOperations::addMBKdx, py::arg("df"), py::arg("m"), py::arg("b"), py::arg("k"), py::arg("clear")=true, py::arg("accumulate")=true, "accumulate $ df += (m M + b B + k K) dx $ (given the latest propagated displacement)");
+    c.def("addMBKv", &MechanicalOperations::addMBKv, py::arg("df"), py::arg("m"), py::arg("b"), py::arg("k"), py::arg("clear")=true, py::arg("accumulate")=true, "accumulate $ df += (m M + b B + k K) velocity $");
+    c.def("addSeparateGravity", &MechanicalOperations::addSeparateGravity, py::arg("dt"), py::arg("result") = sofa::core::VecDerivId::velocity(), "Add dt*Gravity to the velocity");
+    c.def("computeContactForce", &MechanicalOperations::computeContactForce, py::arg("result"));
+    c.def("computeContactDf", &MechanicalOperations::computeContactDf, py::arg("df"));
+    c.def("computeAcc", &MechanicalOperations::computeAcc, py::arg("t"), py::arg("a"), py::arg("x"), py::arg("v"), "Compute a(x,v) at time t. Parameters x and v not const due to propagation through mappings.");
+    c.def("computeContactAcc", &MechanicalOperations::computeContactAcc, py::arg("t"), py::arg("a"), py::arg("x"), py::arg("v"), "Parameters x and v not const due to propagation through mappings.");
+
+    // Matrix operations using LinearSolver components
+    c.def("m_resetSystem", &MechanicalOperations::m_resetSystem);
+    c.def("m_setSystemMBKMatrix", &MechanicalOperations::m_setSystemMBKMatrix, py::arg("m"), py::arg("b"), py::arg("k"));
+    c.def("m_setSystemRHVector", &MechanicalOperations::m_setSystemRHVector, py::arg("f"));
+    c.def("m_setSystemLHVector", &MechanicalOperations::m_setSystemLHVector, py::arg("dx"));
+    c.def("m_solveSystem", &MechanicalOperations::m_solveSystem);
+
+    // Constraints
+    c.def("solveConstraint", &MechanicalOperations::solveConstraint, py::arg("id"), py::arg("order"));
+
+    // Matrix operations
+    c.def("getMatrixDimension", py::overload_cast<sofa::Size*, sofa::Size*, sofa::core::behavior::MultiMatrixAccessor*>(&MechanicalOperations::getMatrixDimension), py::arg("nrows"), py::arg("ncols"), py::arg("matrix") = nullptr);
+    c.def("getMatrixDimension", py::overload_cast<sofa::core::behavior::MultiMatrixAccessor*>(&MechanicalOperations::getMatrixDimension), py::arg("matrix"));
+    c.def("addMBK_ToMatrix", &MechanicalOperations::addMBK_ToMatrix, py::arg("matrix"), py::arg("mFact"), py::arg("bFact"), py::arg("kFact"));
+    c.def("multiVector2BaseVector", &MechanicalOperations::multiVector2BaseVector, py::arg("src"), py::arg("dest"), py::arg("matrix"));
+    c.def("baseVector2MultiVector", &MechanicalOperations::baseVector2MultiVector, py::arg("src"), py::arg("dest"), py::arg("matrix"));
+    c.def("multiVectorPeqBaseVector", &MechanicalOperations::multiVectorPeqBaseVector, py::arg("dest"), py::arg("src"), py::arg("matrix"));
+}
+
+}

bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_VectorOperations.cpp

@@ -0,0 +1,67 @@
+/******************************************************************************
+*                              SofaPython3 plugin                             *
+*                  (c) 2021 CNRS, University of Lille, INRIA                  *
+*                                                                             *
+* This program 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 program 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 program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Contact information: [email protected]                             *
+******************************************************************************/
+
+#include <pybind11/pybind11.h>
+
+#include <sofa/simulation/config.h>
+#include <sofa/simulation/VectorOperations.h>
+
+namespace py { using namespace pybind11; }
+
+namespace sofapython3 {
+
+void moduleAddVectorOperations(pybind11::module& m)
+{
+    using namespace sofa::simulation::common;
+
+    py::class_<VectorOperations> c (m, "VectorOperations");
+    c.def(py::init([](sofa::core::objectmodel::BaseContext* context, bool precomputedTraversalOrder = false) {
+        return VectorOperations(sofa::core::ExecParams::defaultInstance(), context, precomputedTraversalOrder);
+    }), py::arg("context"), py::arg("precomputedTraversalOrder") = false);
+
+    c.def("v_alloc", py::overload_cast<sofa::core::MultiVecCoordId&>(&VectorOperations::v_alloc), py::arg("multi_vector_id"));
+    c.def("v_alloc", py::overload_cast<sofa::core::MultiVecDerivId&>(&VectorOperations::v_alloc), py::arg("multi_vector_id"));
+
+    c.def("v_free", py::overload_cast<sofa::core::MultiVecCoordId&, bool, bool>(&VectorOperations::v_free), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false);
+    c.def("v_free", py::overload_cast<sofa::core::MultiVecDerivId&, bool, bool>(&VectorOperations::v_free), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false);
+
+    c.def("v_realloc", py::overload_cast<sofa::core::MultiVecCoordId&, bool, bool>(&VectorOperations::v_realloc), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false);
+    c.def("v_realloc", py::overload_cast<sofa::core::MultiVecDerivId&, bool, bool>(&VectorOperations::v_realloc), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false);
+
+    c.def("v_clear", &VectorOperations::v_clear, py::arg("multi_vector_id"));
+    c.def("v_size", &VectorOperations::v_size, py::arg("multi_vector_id"));
+
+    c.def("v_eq",  py::overload_cast<sofa::core::MultiVecId, sofa::core::ConstMultiVecId>(&VectorOperations::v_eq), py::arg("v").noconvert(false), py::arg("a").noconvert(false), "v = a");
+    c.def("v_eq", [](VectorOperations & self, sofa::core::MultiVecId &v, sofa::core::ConstMultiVecId & a){
+        self.v_eq(v, a);
+    });
+    c.def("v_eq",  py::overload_cast<sofa::core::MultiVecId, sofa::core::ConstMultiVecId, SReal>(&VectorOperations::v_eq), py::arg("v").noconvert(false), py::arg("a").noconvert(false), py::arg("f").noconvert(false), "v = f*a");
+    c.def("v_peq",  &VectorOperations::v_peq, py::arg("v"), py::arg("a"), py::arg("f")=1.0, "v += f*a");
+    c.def("v_teq",  &VectorOperations::v_teq, py::arg("v"), py::arg("f"), "v *= f");
+    c.def("v_op",  &VectorOperations::v_op, py::arg("v"), py::arg("a"), py::arg("b"), py::arg("f"), "v=a+b*f");
+
+    c.def("v_dot", &VectorOperations::v_dot, py::arg("a"), py::arg("b"), "a dot b ( get result using finish )");
+    c.def("v_norm", &VectorOperations::v_norm, py::arg("a"), py::arg("l"), "Compute the norm of a vector ( get result using finish ). The type of norm is set by parameter l. Use 0 for the infinite norm. Note that the 2-norm is more efficiently computed using the square root of the dot product.");
+    c.def("finish", &VectorOperations::finish);
+
+    c.def("v_threshold", &VectorOperations::v_threshold, py::arg("a"), py::arg("threshold"), "Nullify the values below the given threshold.");
+}
+
+}

bindings/Modules/src/SofaPython3/SofaSimulationCore/CMakeLists.txt

@@ -0,0 +1,22 @@
+project(Bindings.Modules.SofaSimulationCore)
+
+set(SOURCE_FILES
+    ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MechanicalOperations.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Binding_VectorOperations.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Module_SofaSimulationCore.cpp
+)
+
+set(HEADER_FILES
+)
+
+find_package(Sofa.SimulationCore REQUIRED)
+
+SP3_add_python_module(
+    TARGET       ${PROJECT_NAME}
+    PACKAGE      Bindings
+    MODULE       SofaSimulationCore
+    DESTINATION  Sofa
+    SOURCES      ${SOURCE_FILES}
+    HEADERS      ${HEADER_FILES}
+    DEPENDS      SofaSimulationCore SofaPython3::Plugin
+)

bindings/Modules/src/SofaPython3/SofaSimulationCore/Module_SofaSimulationCore.cpp

@@ -0,0 +1,38 @@
+/******************************************************************************
+*                              SofaPython3 plugin                             *
+*                  (c) 2021 CNRS, University of Lille, INRIA                  *
+*                                                                             *
+* This program 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 program 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 program. If not, see <http://www.gnu.org/licenses/>.        *
+*******************************************************************************
+* Contact information: [email protected]                             *
+******************************************************************************/
+
+#include <pybind11/pybind11.h>
+#include <sofa/simulation/init.h>
+
+namespace py { using namespace pybind11; }
+
+namespace sofapython3
+{
+
+void moduleAddMechanicalOperations(pybind11::module& m);
+void moduleAddVectorOperations(pybind11::module& m);
+
+PYBIND11_MODULE(SofaSimulationCore, m)
+{
+    moduleAddMechanicalOperations(m);
+    moduleAddVectorOperations(m);
+}
+
+} // namespace sofapython3

bindings/Modules/tests/CMakeLists.txt

@@ -7,6 +7,8 @@ set(SOURCE_FILES
 set(PYTHON_FILES
     ${CMAKE_CURRENT_SOURCE_DIR}/SofaDeformable/LinearSpring.py
     ${CMAKE_CURRENT_SOURCE_DIR}/SofaDeformable/SpringForceField.py
+    ${CMAKE_CURRENT_SOURCE_DIR}/SofaSimulationCore/MechanicalOperations.py
+    ${CMAKE_CURRENT_SOURCE_DIR}/SofaSimulationCore/VectorOperations.py
 )
 
 find_package(SofaGTestMain REQUIRED)
@@ -26,7 +28,7 @@ sofa_auto_set_target_rpath(
 
 add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME})
 
-set(DIR_BINDING_LIST SofaBaseTopology SofaDeformable)
+set(DIR_BINDING_LIST SofaBaseTopology SofaDeformable SofaSimulationCore)
 get_property(_isMultiConfig GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG)
 foreach(dir_binding ${DIR_BINDING_LIST})
     if (NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${dir_binding}")