diff --git a/bindings/Modules/CMakeLists.txt b/bindings/Modules/CMakeLists.txt index d653248e..1527fd76 100644 --- a/bindings/Modules/CMakeLists.txt +++ b/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) diff --git a/bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_MechanicalOperations.cpp b/bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_MechanicalOperations.cpp new file mode 100644 index 00000000..1ac2f990 --- /dev/null +++ b/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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include + +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { + +void moduleAddMechanicalOperations(pybind11::module& m) +{ + using namespace sofa::simulation::common; + + py::class_ c (m, "MechanicalOperations"); + c.def(py::init(), 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(&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(&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(&MechanicalOperations::getMatrixDimension), py::arg("nrows"), py::arg("ncols"), py::arg("matrix") = nullptr); + c.def("getMatrixDimension", py::overload_cast(&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")); +} + +} diff --git a/bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_VectorOperations.cpp b/bindings/Modules/src/SofaPython3/SofaSimulationCore/Binding_VectorOperations.cpp new file mode 100644 index 00000000..9ac478ed --- /dev/null +++ b/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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include + +#include +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { + +void moduleAddVectorOperations(pybind11::module& m) +{ + using namespace sofa::simulation::common; + + py::class_ 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(&VectorOperations::v_alloc), py::arg("multi_vector_id")); + c.def("v_alloc", py::overload_cast(&VectorOperations::v_alloc), py::arg("multi_vector_id")); + + c.def("v_free", py::overload_cast(&VectorOperations::v_free), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false); + c.def("v_free", py::overload_cast(&VectorOperations::v_free), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false); + + c.def("v_realloc", py::overload_cast(&VectorOperations::v_realloc), py::arg("multi_vector_id"), py::arg("interactionForceField") = false, py::arg("propagate") = false); + c.def("v_realloc", py::overload_cast(&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(&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(&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."); +} + +} \ No newline at end of file diff --git a/bindings/Modules/src/SofaPython3/SofaSimulationCore/CMakeLists.txt b/bindings/Modules/src/SofaPython3/SofaSimulationCore/CMakeLists.txt new file mode 100644 index 00000000..bf261d97 --- /dev/null +++ b/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 +) diff --git a/bindings/Modules/src/SofaPython3/SofaSimulationCore/Module_SofaSimulationCore.cpp b/bindings/Modules/src/SofaPython3/SofaSimulationCore/Module_SofaSimulationCore.cpp new file mode 100644 index 00000000..b222b774 --- /dev/null +++ b/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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include + +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 \ No newline at end of file diff --git a/bindings/Modules/tests/CMakeLists.txt b/bindings/Modules/tests/CMakeLists.txt index e05e9c73..ebfa4f88 100644 --- a/bindings/Modules/tests/CMakeLists.txt +++ b/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}") diff --git a/bindings/Modules/tests/SofaSimulationCore/MechanicalOperations.py b/bindings/Modules/tests/SofaSimulationCore/MechanicalOperations.py new file mode 100644 index 00000000..c03a3f6a --- /dev/null +++ b/bindings/Modules/tests/SofaSimulationCore/MechanicalOperations.py @@ -0,0 +1,85 @@ +import unittest +import Sofa +from Sofa import SofaSimulationCore +import numpy as np + +""" +In order to test out a good subset of the Mechanical operations, let's do here do a Newton-Raphson solver in python. +""" + + +class StaticOdeSolver (Sofa.Core.OdeSolver): + def solve(self, _, dt, X, V): + mparams = Sofa.Core.MechanicalParams() + vop = SofaSimulationCore.VectorOperations(self.getContext()) + mop = SofaSimulationCore.MechanicalOperations(mparams, self.getContext()) + + # Allocate the solution vector + dx = Sofa.Core.VecId.dx() + vop.v_realloc(dx, interactionForceField=True, propagate=True) + vop.v_clear(dx) + + # Compute the residual + F = Sofa.Core.VecId.force() + mop.computeForce(result=F, clear=True, accumulate=True); + mop.projectResponse(F) + vop.v_dot(F, F) + F_norm = np.sqrt(vop.finish()) + + # Assemble the system + mop.m_setSystemMBKMatrix(m=0, b=0, k=-1) + mop.m_setSystemRHVector(F) + mop.m_setSystemLHVector(dx) + + # Solve the system + mop.m_solveSystem() + vop.v_dot(dx, dx) + dx_norm = np.sqrt(vop.finish()) + + # Propagate the solution + vop.v_peq(X, dx) # X += dx + + # Solve the constraints + # todo: Bind ConstraintParams + # mop.solveConstraint(X, ConstraintParams.ConstOrder.POS) + + print(f"Solved with |F| = {F_norm} and |dx| = {dx_norm}") + + +class MechanicalOperations(unittest.TestCase): + def test_static_solver(self): + root = Sofa.Core.Node() + createScene(root) + Sofa.Simulation.init(root) + for _ in range(5): + Sofa.Simulation.animate(root, root.dt.value) + + middle_node_index = 52 + solution_of_middle_node = np.array([0, 10.0953, -0.285531]) + self.assertLess(np.linalg.norm(solution_of_middle_node - root.mechanics.neumann.mo.position.array()[middle_node_index])/np.linalg.norm(solution_of_middle_node), 1e-5) + + +def createScene(root): + w, l = 2, 10 + nx, ny = 3, 9 + dx, dy = w/nx/2, l/ny/2 + root.addObject('RequiredPlugin', pluginName='SofaLoader SofaBoundaryCondition SofaEngine SofaSimpleFem SofaSparseSolver SofaTopologyMapping') + root.addObject('MeshObjLoader', name='surface', filename='mesh/cylinder.obj') + root.addObject('SparseGridTopology', src='@surface', name='grid', n=[nx, ny, nx]) + root.addChild('mechanics') + root.mechanics.addObject(StaticOdeSolver(name='solver')) + root.mechanics.addObject('SparseLDLSolver') + root.mechanics.addObject('MechanicalObject', name='mo', src='@../grid') + root.mechanics.addObject('HexahedronSetTopologyContainer', name='hexa_topology', hexahedra='@../grid.hexahedra') + root.mechanics.addObject('HexahedronFEMForceField', youngModulus=3000, poissonRatio=0.3) + root.mechanics.addObject('BoxROI', name='fixed_roi', box=[-1-dx,0-dx,-1-dx, 1+dx, 0+dx, 1+dx], drawBoxes=True) + root.mechanics.addObject('FixedConstraint', indices='@fixed_roi.indices') + root.mechanics.addChild('neumann') + root.mechanics.neumann.addObject('MechanicalObject', name='mo') + root.mechanics.neumann.addObject('QuadSetTopologyContainer', name='quad_topology') + root.mechanics.neumann.addObject('QuadSetTopologyModifier') + root.mechanics.neumann.addObject('QuadSetGeometryAlgorithms') + root.mechanics.neumann.addObject('Hexa2QuadTopologicalMapping', input='@../hexa_topology', output='@quad_topology') + root.mechanics.neumann.addObject('SubsetMapping', applyRestPosition=True, input='@../mo', output='@./mo', indices='@quad_topology.points') + root.mechanics.neumann.addObject('BoxROI', name='top_roi', quad='@quad_topology.quads', src='@mo', box=[-1-dx,10-dx,-1-dx, 1+dx, 10+dx, 1+dx], drawBoxes=True) + root.mechanics.neumann.addObject('QuadPressureForceField', pressure=[0, 0, -1], quadList='@top_roi.quadIndices') diff --git a/bindings/Modules/tests/SofaSimulationCore/VectorOperations.py b/bindings/Modules/tests/SofaSimulationCore/VectorOperations.py new file mode 100644 index 00000000..e92973aa --- /dev/null +++ b/bindings/Modules/tests/SofaSimulationCore/VectorOperations.py @@ -0,0 +1,103 @@ +import unittest +import Sofa +from Sofa import SofaSimulationCore +import numpy as np + + +class VectorOperations(unittest.TestCase): + def test_alloc_coord(self): + root = Sofa.Core.Node() + create_scene(root) + Sofa.Simulation.init(root) + + vop = SofaSimulationCore.VectorOperations(root) + + # Create a first multivector + multi_vec_id_1 = Sofa.Core.MultiVecCoordId() + self.assertEqual(multi_vec_id_1.defaultId.index, 0) + self.assertTrue(multi_vec_id_1.isNull()) + vop.v_alloc(multi_vec_id_1) + self.assertFalse(multi_vec_id_1.isNull()) + self.assertEqual(multi_vec_id_1.getId(root.node_1.mo).index, 5) + self.assertEqual(multi_vec_id_1.getId(root.node_2.mo).index, 5) + + # Create a second one + multi_vec_id_2 = Sofa.Core.MultiVecCoordId() + vop.v_alloc(multi_vec_id_2) + self.assertEqual(multi_vec_id_2.getId(root.node_1.mo).index, 6) + self.assertEqual(multi_vec_id_2.getId(root.node_2.mo).index, 6) + + # Create a third one, only in node_2 + vop2 = SofaSimulationCore.VectorOperations(root.node_2) + multi_vec_id_3 = Sofa.Core.MultiVecCoordId() + vop2.v_alloc(multi_vec_id_3) + self.assertEqual(multi_vec_id_3.getId(root.node_1.mo).index, 0) + self.assertEqual(multi_vec_id_3.getId(root.node_2.mo).index, 7) + + # Free + vop.v_free(multi_vec_id_1) + vop.v_free(multi_vec_id_2) + vop2.v_free(multi_vec_id_3) + + def test_alloc_deriv(self): + root = Sofa.Core.Node() + create_scene(root) + Sofa.Simulation.init(root) + + vop = SofaSimulationCore.VectorOperations(root) + + # Create a first multivector + multi_vec_id_1 = Sofa.Core.MultiVecDerivId() + self.assertEqual(multi_vec_id_1.defaultId.index, 0) + self.assertTrue(multi_vec_id_1.isNull()) + vop.v_alloc(multi_vec_id_1) + self.assertFalse(multi_vec_id_1.isNull()) + self.assertEqual(multi_vec_id_1.getId(root.node_1.mo).index, 9) + self.assertEqual(multi_vec_id_1.getId(root.node_2.mo).index, 9) + + # Create a second one + multi_vec_id_2 = Sofa.Core.MultiVecDerivId() + vop.v_alloc(multi_vec_id_2) + self.assertEqual(multi_vec_id_2.getId(root.node_1.mo).index, 10) + self.assertEqual(multi_vec_id_2.getId(root.node_2.mo).index, 10) + + # Create a third one, only in node_2 + vop2 = SofaSimulationCore.VectorOperations(root.node_2) + multi_vec_id_3 = Sofa.Core.MultiVecDerivId() + vop2.v_alloc(multi_vec_id_3) + self.assertEqual(multi_vec_id_3.getId(root.node_1.mo).index, 0) + self.assertEqual(multi_vec_id_3.getId(root.node_2.mo).index, 11) + + # Free + vop.v_free(multi_vec_id_1) + vop.v_free(multi_vec_id_2) + vop2.v_free(multi_vec_id_3) + + def test_op(self): + root = Sofa.Core.Node() + create_scene(root) + Sofa.Simulation.init(root) + + vop = SofaSimulationCore.VectorOperations(root) + + # Create a temporary multivector + v = Sofa.Core.MultiVecCoordId() + vop.v_alloc(v) + + x_id = Sofa.Core.VecCoordId.position() + vop.v_eq(v, x_id) # v = x + vop.v_peq(v, x_id, 2) # v += 2*x + vop.v_dot(v, v) # n = sqrt(v dot v) + n = np.sqrt(vop.finish()) + + x = np.concatenate((root.node_1.mo.position.array(), root.node_2.mo.position.array()), axis=None).flatten() + self.assertEqual(n, np.linalg.norm(x + (2*x))) + + +def create_scene(root): + root.addObject('RequiredPlugin', pluginName='SofaBaseMechanics') + root.addChild('node_1') + root.node_1.addObject('MechanicalObject', name='mo', template='Vec3', position=[[1,1,1], [2,2,2], [3,3,3]]) + + root.addChild('node_2') + root.node_2.addObject('MechanicalObject', name='mo', template='Vec2', position=[[1,1], [2,2], [3,3], [4, 4]]) diff --git a/bindings/Modules/tests/main.cpp b/bindings/Modules/tests/main.cpp index 98eaa5d6..1dae83c1 100644 --- a/bindings/Modules/tests/main.cpp +++ b/bindings/Modules/tests/main.cpp @@ -38,8 +38,8 @@ static struct Tests : public sofapython3::PythonTestExtractor MessageDispatcher::addHandler(&MainPerComponentLoggingMessageHandler::getInstance()) ; const std::string executable_directory = sofa::helper::Utils::getExecutableDirectory(); - addTestDirectory(executable_directory+"/Bindings.Modules.Tests.d/SofaBaseTopology", "SofaBaseTopology_"); addTestDirectory(executable_directory+"/Bindings.Modules.Tests.d/SofaDeformable", "SofaDeformable_"); + addTestDirectory(executable_directory+"/Bindings.Modules.Tests.d/SofaSimulationCore", "SofaSimulationCore_"); } } python_tests; diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.cpp new file mode 100644 index 00000000..312b2fa9 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.cpp @@ -0,0 +1,45 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include +#include + +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { + +void moduleAddBaseState(py::module& m) { + using namespace sofa::core; + py::class_> (m, "BaseState") + .def("getSize", &BaseState::getSize) + .def("resize", &BaseState::getSize) + ; + + /// register the BaseState binding in the downcasting subsystem + PythonFactory::registerType([](sofa::core::objectmodel::Base* object) + { + return py::cast(dynamic_cast(object)); + }); +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.h new file mode 100644 index 00000000..be1fde0d --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_BaseState.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddBaseState(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.cpp new file mode 100644 index 00000000..b0f04793 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.cpp @@ -0,0 +1,53 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { + +void moduleAddExecParams(py::module& m) { + using namespace sofa::core; + py::class_ c (m, "ExecParams", "Class gathering parameters use by most components methods, and transmitted by all visitors"); + c.def(py::init()); + + py::enum_ (c, "ExecMode") + .value("EXEC_NONE", ExecParams::ExecMode::EXEC_NONE) + .value("EXEC_DEFAULT", ExecParams::ExecMode::EXEC_DEFAULT) + .value("EXEC_DEBUG", ExecParams::ExecMode::EXEC_DEBUG) + .value("EXEC_GPU", ExecParams::ExecMode::EXEC_GPU) + .value("EXEC_GRAPH", ExecParams::ExecMode::EXEC_GRAPH) + .export_values() + ; + + // Public properties/methods + c.def("checkValidStorage", &ExecParams::checkValidStorage); + c.def("execMode", &ExecParams::execMode, "Mode of execution requested"); + c.def("threadID", &ExecParams::threadID, "Index of current thread (0 corresponding to the only thread in sequential mode, or first thread in parallel mode)"); + c.def("nbThreads", &ExecParams::nbThreads, "Number of threads currently known to Sofa"); + c.def("update", &ExecParams::update, "Make sure this instance is up-to-date relative to the current thread"); + c.def("setExecMode", &ExecParams::setExecMode, "Request a specific mode of execution"); + c.def("setThreadID", &ExecParams::setThreadID, "Specify the index of the current thread"); + +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.h new file mode 100644 index 00000000..c4cc3ef8 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_ExecParams.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddExecParams(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.cpp new file mode 100644 index 00000000..ee56c14a --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.cpp @@ -0,0 +1,63 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { + +void moduleAddMechanicalParams(py::module& m) { + using namespace sofa::core; + py::class_ c (m, "MechanicalParams", doc::mechanicalparams::MechanicalParamsClass); + c.def(py::init()); + + // Dt + c.def_property("dt", &MechanicalParams::dt, &MechanicalParams::setDt); + c.def("setDt", &MechanicalParams::setDt, py::arg("dt")); + + // Implicit + c.def_property("implicit", &MechanicalParams::implicit, &MechanicalParams::setImplicit, "Specify if the time integration scheme is implicit."); + c.def("setImplicit", &MechanicalParams::setImplicit, py::arg("v"), "Specify if the time integration scheme is implicit."); + + // Mass factor + c.def_property("mFactor", &MechanicalParams::mFactor, &MechanicalParams::setMFactor, "Set Mass matrix contributions factor (for implicit schemes)."); + c.def("setMFactor", &MechanicalParams::setMFactor, py::arg("m"), "Set Mass matrix contributions factor (for implicit schemes)."); + + // Damping factor + c.def_property("bFactor", &MechanicalParams::bFactor, &MechanicalParams::setBFactor, "Set Damping matrix contributions factor (for implicit schemes)."); + c.def("setBFactor", &MechanicalParams::setBFactor, py::arg("b"), "Set Damping matrix contributions factor (for implicit schemes)."); + + // Stiffness factor + c.def_property("kFactor", &MechanicalParams::kFactor, &MechanicalParams::setKFactor, "Set Stiffness matrix contributions factor (for implicit schemes)."); + c.def("setKFactor", &MechanicalParams::setKFactor, py::arg("k"), "Set Stiffness matrix contributions factor (for implicit schemes)."); + + // Symmetric + c.def_property("symmetricMatrix", &MechanicalParams::symmetricMatrix, &MechanicalParams::setSymmetricMatrix, "Set the symmetric matrix flag (for implicit schemes), for solvers specialized on symmetric matrices"); + c.def("setSymmetricMatrix", &MechanicalParams::setSymmetricMatrix, "Set the symmetric matrix flag (for implicit schemes), for solvers specialized on symmetric matrices"); + + // Energy + c.def_property("energy", &MechanicalParams::energy, &MechanicalParams::setEnergy, "Specify if the potential and kinematic energies should be computed."); + c.def("setEnergy", &MechanicalParams::setEnergy, "Specify if the potential and kinematic energies should be computed."); +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.h new file mode 100644 index 00000000..0cbcd507 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddMechanicalParams(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams_doc.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams_doc.h new file mode 100644 index 00000000..9868e224 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MechanicalParams_doc.h @@ -0,0 +1,32 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once + +namespace sofapython3::doc::mechanicalparams +{ + +static auto MechanicalParamsClass = +R"( +Class gathering parameters use by mechanical components methods, and transmitted by mechanical visitors +)"; + + +} diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.cpp new file mode 100644 index 00000000..91787c97 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.cpp @@ -0,0 +1,125 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include +#include + +namespace py { using namespace pybind11; } +using namespace sofa::core; + +namespace sofapython3 { + +template +void bindMultiVecId(py::module& m, const std::string & name) { + using MultiVecId = TMultiVecId; + py::class_ c (m, name.c_str()); + + // Constructors + c.def(py::init()); + + // Read-only can be constructed from both read-only and read-write multi vectors + if constexpr (vaccess == V_READ) { + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + if constexpr(vtype != V_ALL) { + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + } else { + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + } + } + + // Construct from the same access type + c.def(py::init&>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init&>()); + + // Construct from a ALL vector type + if constexpr(vtype != V_ALL) { + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + } else { + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + + c.def(py::init &>()); + py::implicitly_convertible, TMultiVecId>(); + } + + // Public methods/properties + c.def_property_readonly("defaultId", &MultiVecId::getDefaultId); + c.def("getDefaultId", &MultiVecId::getDefaultId); + + c.def("getName", &MultiVecId::getName); + c.def("isNull", &MultiVecId::isNull); + c.def("__str__", &MultiVecId::getName); + c.def("getId", &MultiVecId::getId, py::arg("state")); + + // Static methods/properties + c.def_static("null", &MultiVecId::null); +} + +void moduleAddMultiVecId(py::module& m) { + bindMultiVecId (m, "ConstMultiVecId"); + bindMultiVecId (m, "MultiVecId"); + bindMultiVecId (m, "ConstMultiVecCoordId"); + bindMultiVecId (m, "MultiVecCoordId"); + bindMultiVecId (m, "ConstMultiVecDerivId"); + bindMultiVecId (m, "MultiVecDerivId"); + bindMultiVecId (m, "ConstMultiMatrixDerivId"); + bindMultiVecId (m, "MultiMatrixDerivId"); +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.h new file mode 100644 index 00000000..05b79993 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_MultiVecId.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddMultiVecId(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.cpp new file mode 100644 index 00000000..84325136 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.cpp @@ -0,0 +1,64 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include +#include +#include +#include + +#include +#include + +namespace py { using namespace pybind11; } + +namespace sofapython3 { +using namespace sofa::core::behavior; + +class PyOdeSolver : public OdeSolver { +public: + SOFA_CLASS(PyOdeSolver, OdeSolver); + void solve(const sofa::core::ExecParams* params, SReal dt, sofa::core::MultiVecCoordId xResult, sofa::core::MultiVecDerivId vResult) final { + PythonEnvironment::gil acquire; + PYBIND11_OVERLOAD_PURE(void, OdeSolver, solve, params, dt, xResult, vResult); + } +}; + +void moduleAddOdeSolver(py::module& m) { + py::class_> (m, "OdeSolver", doc::odesolver::OdeSolverClass) + .def(py::init([](py::args & /*args*/, py::kwargs & kwargs) { + auto o = sofa::core::sptr (new PyOdeSolver()); + for (auto kv : kwargs) { + auto key = py::cast(kv.first); + auto value = py::reinterpret_borrow(kv.second); + BindingBase::SetAttr(*o.get(), key, value); + } + return o; + })) + ; + + /// register the BaseState binding in the downcasting subsystem + PythonFactory::registerType([](sofa::core::objectmodel::Base* object) + { + return py::cast(dynamic_cast(object)); + }); +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.h new file mode 100644 index 00000000..1e8b6df1 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddOdeSolver(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver_doc.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver_doc.h new file mode 100644 index 00000000..ffccde27 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_OdeSolver_doc.h @@ -0,0 +1,45 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once + +namespace sofapython3::doc::odesolver +{ + +static auto OdeSolverClass = +R"( +Component responsible for timestep integration, i.e. advancing the state from time t to t+dt. + +This class currently control both the integration scheme (explicit, +implicit, static, etc), and the linear system resolution algorithm +(conjugate gradient, matrix direct inversion, etc). Those two aspect will +propably be separated in a future version. + +While all computations required to do the integration step are handled by +this object, they should not be implemented directly in it, but instead +the solver propagates orders (or Visitor) to the other components in the +scenegraph that will locally execute them. This allow for greater +flexibility (the solver can just ask for the forces to be computed without +knowing what type of forces are present), as well as performances +(some computations can be executed in parallel). +)"; + + +} diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.cpp new file mode 100644 index 00000000..aec2d272 --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.cpp @@ -0,0 +1,107 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include +#include + +namespace py { using namespace pybind11; } +using namespace sofa::core; + +namespace sofapython3 { + +template +void bindVecId(py::module& m, const std::string & name) { + + using StandardVec = TStandardVec; + py::class_ s (m, ("Standard"+name).c_str()); + if constexpr (vtype == VecType::V_COORD || vtype == VecType::V_ALL) { + s.def_static("position", &StandardVec::position); + s.def_static("restPosition", &StandardVec::restPosition); + s.def_static("freePosition", &StandardVec::freePosition); + s.def_static("resetPosition", &StandardVec::resetPosition); + } + if constexpr (vtype == VecType::V_DERIV || vtype == VecType::V_ALL) { + s.def_static("velocity", &StandardVec::velocity); + s.def_static("resetVelocity", &StandardVec::resetVelocity); + s.def_static("freeVelocity", &StandardVec::freeVelocity); + s.def_static("normal", &StandardVec::normal); + s.def_static("force", &StandardVec::force); + s.def_static("externalForce", &StandardVec::externalForce); + s.def_static("dx", &StandardVec::dx); + s.def_static("dforce", &StandardVec::dforce); + } + if constexpr(vtype == VecType::V_MATDERIV || vtype == VecType::V_ALL) { + s.def_static("constraintJacobian", &StandardVec::constraintJacobian); + s.def_static("mappingJacobian", &StandardVec::mappingJacobian); + } + + if constexpr (vtype == VecType::V_ALL) { + s.def_static("getFirstDynamicIndex", &StandardVec::getFirstDynamicIndex, py::arg("t")); + } + + using VecId = TVecId; + py::class_ v (m, name.c_str()); + v.def(py::init()); + if constexpr (vtype != V_ALL) { + py::implicitly_convertible, TVecId>(); + v.def(py::init()); + } + v.def("getName", &VecId::getName); + v.def("__str__", &VecId::getName); +} + +void moduleAddVecId(py::module& m) { + + // VecType + py::enum_ (m, "VecType") + .value("V_ALL", VecType::V_ALL) + .value("V_COORD", VecType::V_COORD) + .value("V_DERIV", VecType::V_DERIV) + .value("V_MATDERIV", VecType::V_MATDERIV) + .export_values() + ; + + // VecAccess + py::enum_ (m, "VecAccess") + .value("V_READ", VecAccess::V_READ) + .value("V_WRITE", VecAccess::V_WRITE) + .export_values() + ; + + // BaseVecId + py::class_ (m, "BaseVecId") + .def_property_readonly("index", &BaseVecId::getIndex) + .def("getIndex", &BaseVecId::getIndex) + .def_property_readonly("type", &BaseVecId::getType) + .def("getType", &BaseVecId::getType) + ; + + // VecId + bindVecId (m, "ConstVecId"); + bindVecId (m, "VecId"); + bindVecId (m, "ConstVecCoordId"); + bindVecId(m, "VecCoordId"); + bindVecId (m, "ConstVecDerivId"); + bindVecId(m, "VecDerivId"); + bindVecId (m, "ConstMatrixDerivId"); + bindVecId(m, "MatrixDerivId"); +} + +} // namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.h b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.h new file mode 100644 index 00000000..d518acdc --- /dev/null +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Binding_VecId.h @@ -0,0 +1,28 @@ +/****************************************************************************** +* 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 . * +******************************************************************************* +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#pragma once +#include + +namespace sofapython3 { + +void moduleAddVecId(pybind11::module &m); + +} /// namespace sofapython3 diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/CMakeLists.txt b/bindings/Sofa/src/SofaPython3/Sofa/Core/CMakeLists.txt index 2019180d..976726cb 100644 --- a/bindings/Sofa/src/SofaPython3/Sofa/Core/CMakeLists.txt +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/CMakeLists.txt @@ -10,6 +10,7 @@ set(HEADER_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseObject_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseCamera.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseContext.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseState.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ContactListener.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ContactListener_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Context.h @@ -17,10 +18,16 @@ set(HEADER_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Controller_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_DataEngine.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_DataEngine_doc.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ExecParams.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ForceField.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ForceField_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ObjectFactory.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_OdeSolver.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_OdeSolver_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ObjectFactory_doc.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MechanicalParams.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MechanicalParams_doc.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MultiVecId.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Node.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Node_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_NodeIterator.h @@ -39,6 +46,7 @@ set(HEADER_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseData_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseCamera_doc.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Topology.h + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_VecId.h ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseMeshTopology.h ) @@ -49,6 +57,7 @@ set(SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseObject.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseCamera.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseContext.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseState.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ContactListener.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Context.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Controller.cpp @@ -57,8 +66,12 @@ set(SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Data/Binding_DataString.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Data/Binding_DataLink.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Data/Binding_DataVectorString.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ExecParams.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ForceField.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_ObjectFactory.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_OdeSolver.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MechanicalParams.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_MultiVecId.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Node.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_NodeIterator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Prefab.cpp @@ -66,6 +79,7 @@ set(SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Binding_PythonScriptEvent.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseLink.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_Topology.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Binding_VecId.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Binding_BaseMeshTopology.cpp ) diff --git a/bindings/Sofa/src/SofaPython3/Sofa/Core/Submodule_Core.cpp b/bindings/Sofa/src/SofaPython3/Sofa/Core/Submodule_Core.cpp index 2e8d6bac..c3b803ac 100644 --- a/bindings/Sofa/src/SofaPython3/Sofa/Core/Submodule_Core.cpp +++ b/bindings/Sofa/src/SofaPython3/Sofa/Core/Submodule_Core.cpp @@ -28,18 +28,24 @@ using sofa::helper::logging::Message; #include #include #include +#include #include #include #include #include #include +#include #include +#include +#include +#include #include #include #include #include #include #include +#include #include #include @@ -126,17 +132,23 @@ PYBIND11_MODULE(Core, core) moduleAddDataVectorString(core); moduleAddBaseObject(core); moduleAddBaseCamera(core); + moduleAddBaseState(core); moduleAddContactListener(core); moduleAddContext(core); moduleAddController(core); moduleAddDataEngine(core); + moduleAddExecParams(core); moduleAddForceField(core); moduleAddObjectFactory(core); + moduleAddOdeSolver(core); + moduleAddMechanicalParams(core); + moduleAddMultiVecId(core); moduleAddNode(core); moduleAddNodeIterator(core); moduleAddPrefab(core); moduleAddBaseLink(core); moduleAddTopology(core); + moduleAddVecId(core); moduleAddBaseMeshTopology(core); } diff --git a/bindings/Sofa/tests/CMakeLists.txt b/bindings/Sofa/tests/CMakeLists.txt index 5bacc9b3..de97d554 100644 --- a/bindings/Sofa/tests/CMakeLists.txt +++ b/bindings/Sofa/tests/CMakeLists.txt @@ -18,8 +18,11 @@ set(PYTHON_FILES ${CMAKE_CURRENT_SOURCE_DIR}/Core/ForceField.py ${CMAKE_CURRENT_SOURCE_DIR}/Core/DataEngine.py ${CMAKE_CURRENT_SOURCE_DIR}/Core/MyRestShapeForceField.py + ${CMAKE_CURRENT_SOURCE_DIR}/Core/MechanicalParams.py + ${CMAKE_CURRENT_SOURCE_DIR}/Core/MultiVecId.py ${CMAKE_CURRENT_SOURCE_DIR}/Core/PythonRestShapeForceField.py ${CMAKE_CURRENT_SOURCE_DIR}/Core/BaseLink.py + ${CMAKE_CURRENT_SOURCE_DIR}/Core/VecId.py ${CMAKE_CURRENT_SOURCE_DIR}/Helper/Message.py ${CMAKE_CURRENT_SOURCE_DIR}/Types/RGBAColor.py ${CMAKE_CURRENT_SOURCE_DIR}/Types/Vec3.py diff --git a/bindings/Sofa/tests/Core/MechanicalParams.py b/bindings/Sofa/tests/Core/MechanicalParams.py new file mode 100644 index 00000000..828009df --- /dev/null +++ b/bindings/Sofa/tests/Core/MechanicalParams.py @@ -0,0 +1,49 @@ +import Sofa +import unittest + + +class Test(unittest.TestCase): + def test_creation(self): + mparams = Sofa.Core.MechanicalParams() + + # dt + mparams.dt = 1.25 + self.assertEqual(mparams.dt, 1.25) + mparams.setDt(1.50) + self.assertEqual(mparams.dt, 1.50) + + # Implicit + mparams.implicit = True + self.assertTrue(mparams.implicit) + mparams.setImplicit(False) + self.assertFalse(mparams.implicit) + + # Mass factor + mparams.mFactor = 1.30 + self.assertEqual(mparams.mFactor, 1.30) + mparams.setMFactor(1.35) + self.assertEqual(mparams.mFactor, 1.35) + + # Damping factor + mparams.bFactor = 1.35 + self.assertEqual(mparams.bFactor, 1.35) + mparams.setBFactor(1.40) + self.assertEqual(mparams.bFactor, 1.40) + + # Stiffness factor + mparams.kFactor = 1.45 + self.assertEqual(mparams.kFactor, 1.45) + mparams.setKFactor(1.50) + self.assertEqual(mparams.kFactor, 1.50) + + # Symmetric + mparams.symmetricMatrix = True + self.assertTrue(mparams.symmetricMatrix) + mparams.setSymmetricMatrix(False) + self.assertFalse(mparams.symmetricMatrix) + + # Energy + mparams.energy = True + self.assertTrue(mparams.energy) + mparams.setEnergy(False) + self.assertFalse(mparams.energy) diff --git a/bindings/Sofa/tests/Core/MultiVecId.py b/bindings/Sofa/tests/Core/MultiVecId.py new file mode 100644 index 00000000..5cf84aff --- /dev/null +++ b/bindings/Sofa/tests/Core/MultiVecId.py @@ -0,0 +1,38 @@ +import Sofa +import unittest + + +class MultiVecId(unittest.TestCase): + """ + Test the creation of MultiVecId + + Note that this test case will not test the allocation inside mechanical object. This is instead tested inside + the vector operations (Sofa.SimulationCore.VectorOperations) test case. + """ + def test_coord(self): + # Const Coord + vec = Sofa.Core.ConstMultiVecCoordId() + self.assertEqual(vec.defaultId.index, 0) + self.assertEqual(vec.getDefaultId().getIndex(), 0) + self.assertEqual(vec.defaultId.type, Sofa.Core.VecType.V_COORD) + + vec = Sofa.Core.ConstMultiVecCoordId(Sofa.Core.ConstVecCoordId(15)) + self.assertEqual(vec.defaultId.index, 15) + self.assertEqual(vec.getDefaultId().getIndex(), 15) + self.assertEqual(vec.defaultId.type, Sofa.Core.VecType.V_COORD) + + vec = Sofa.Core.ConstMultiVecCoordId(Sofa.Core.VecCoordId(15)) + self.assertEqual(vec.defaultId.index, 15) + self.assertEqual(vec.getDefaultId().getIndex(), 15) + self.assertEqual(vec.defaultId.type, Sofa.Core.VecType.V_COORD) + + # Coord + vec = Sofa.Core.MultiVecCoordId() + self.assertEqual(vec.defaultId.index, 0) + self.assertEqual(vec.getDefaultId().getIndex(), 0) + self.assertEqual(vec.defaultId.type, Sofa.Core.VecType.V_COORD) + + vec = Sofa.Core.MultiVecCoordId(Sofa.Core.VecCoordId(15)) + self.assertEqual(vec.defaultId.index, 15) + self.assertEqual(vec.getDefaultId().getIndex(), 15) + self.assertEqual(vec.defaultId.type, Sofa.Core.VecType.V_COORD) diff --git a/bindings/Sofa/tests/Core/VecId.py b/bindings/Sofa/tests/Core/VecId.py new file mode 100644 index 00000000..6e9f4e90 --- /dev/null +++ b/bindings/Sofa/tests/Core/VecId.py @@ -0,0 +1,149 @@ +import Sofa +import unittest + + +class VecId(unittest.TestCase): + def test_coord(self): + # Const Coord + vec = Sofa.Core.ConstVecCoordId() + self.assertEqual(vec.index, 0) + self.assertEqual(vec.getIndex(), 0) + self.assertEqual(vec.type, Sofa.Core.VecType.V_COORD) + + vec = Sofa.Core.ConstVecCoordId(15) + self.assertEqual(vec.index, 15) + self.assertEqual(vec.getIndex(), 15) + self.assertEqual(vec.type, Sofa.Core.VecType.V_COORD) + + # Coord + vec = Sofa.Core.VecCoordId() + self.assertEqual(vec.index, 0) + self.assertEqual(vec.getIndex(), 0) + self.assertEqual(vec.type, Sofa.Core.VecType.V_COORD) + + vec = Sofa.Core.VecCoordId(15) + self.assertEqual(vec.index, 15) + self.assertEqual(vec.getIndex(), 15) + self.assertEqual(vec.type, Sofa.Core.VecType.V_COORD) + + def test_deriv(self): + # Const Deriv + vec = Sofa.Core.ConstVecDerivId() + self.assertEqual(vec.index, 0) + self.assertEqual(vec.getIndex(), 0) + self.assertEqual(vec.type, Sofa.Core.VecType.V_DERIV) + + vec = Sofa.Core.ConstVecDerivId(15) + self.assertEqual(vec.index, 15) + self.assertEqual(vec.getIndex(), 15) + self.assertEqual(vec.type, Sofa.Core.VecType.V_DERIV) + + # Deriv + vec = Sofa.Core.VecDerivId() + self.assertEqual(vec.index, 0) + self.assertEqual(vec.getIndex(), 0) + self.assertEqual(vec.type, Sofa.Core.VecType.V_DERIV) + + vec = Sofa.Core.VecDerivId(15) + self.assertEqual(vec.index, 15) + self.assertEqual(vec.getIndex(), 15) + self.assertEqual(vec.type, Sofa.Core.VecType.V_DERIV) + + def test_default_coord_vectors(self): + # Position + self.assertEqual(Sofa.Core.VecCoordId.position().index, 1) + self.assertEqual(Sofa.Core.VecCoordId.position().getName(), "position(V_COORD)") + self.assertEqual(Sofa.Core.ConstVecCoordId.position().index, 1) + self.assertEqual(Sofa.Core.ConstVecCoordId.position().getName(), "position(V_COORD)") + + # restPosition + self.assertEqual(Sofa.Core.VecCoordId.restPosition().index, 2) + self.assertEqual(Sofa.Core.VecCoordId.restPosition().getName(), "restPosition(V_COORD)") + self.assertEqual(Sofa.Core.ConstVecCoordId.restPosition().index, 2) + self.assertEqual(Sofa.Core.ConstVecCoordId.restPosition().getName(), "restPosition(V_COORD)") + + # freePosition + self.assertEqual(Sofa.Core.VecCoordId.freePosition().index, 3) + self.assertEqual(Sofa.Core.VecCoordId.freePosition().getName(), "freePosition(V_COORD)") + self.assertEqual(Sofa.Core.ConstVecCoordId.freePosition().index, 3) + self.assertEqual(Sofa.Core.ConstVecCoordId.freePosition().getName(), "freePosition(V_COORD)") + + # resetPosition + self.assertEqual(Sofa.Core.VecCoordId.resetPosition().index, 4) + self.assertEqual(Sofa.Core.VecCoordId.resetPosition().getName(), "resetPosition(V_COORD)") + self.assertEqual(Sofa.Core.ConstVecCoordId.resetPosition().index, 4) + self.assertEqual(Sofa.Core.ConstVecCoordId.resetPosition().getName(), "resetPosition(V_COORD)") + + # Dynamic offset + self.assertEqual(Sofa.Core.VecId.getFirstDynamicIndex(Sofa.Core.VecType.V_COORD), 5) + self.assertEqual(Sofa.Core.ConstVecId.getFirstDynamicIndex(Sofa.Core.VecType.V_COORD), 5) + + def test_default_deriv_vectors(self): + # velocity + self.assertEqual(Sofa.Core.VecDerivId.velocity().index, 1) + self.assertEqual(Sofa.Core.VecDerivId.velocity().getName(), "velocity(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.velocity().index, 1) + self.assertEqual(Sofa.Core.ConstVecDerivId.velocity().getName(), "velocity(V_DERIV)") + + # resetVelocity + self.assertEqual(Sofa.Core.VecDerivId.resetVelocity().index, 2) + self.assertEqual(Sofa.Core.VecDerivId.resetVelocity().getName(), "resetVelocity(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.resetVelocity().index, 2) + self.assertEqual(Sofa.Core.ConstVecDerivId.resetVelocity().getName(), "resetVelocity(V_DERIV)") + + # freeVelocity + self.assertEqual(Sofa.Core.VecDerivId.freeVelocity().index, 3) + self.assertEqual(Sofa.Core.VecDerivId.freeVelocity().getName(), "freeVelocity(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.freeVelocity().index, 3) + self.assertEqual(Sofa.Core.ConstVecDerivId.freeVelocity().getName(), "freeVelocity(V_DERIV)") + + # normal + self.assertEqual(Sofa.Core.VecDerivId.normal().index, 4) + self.assertEqual(Sofa.Core.VecDerivId.normal().getName(), "normal(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.normal().index, 4) + self.assertEqual(Sofa.Core.ConstVecDerivId.normal().getName(), "normal(V_DERIV)") + + # force + self.assertEqual(Sofa.Core.VecDerivId.force().index, 5) + self.assertEqual(Sofa.Core.VecDerivId.force().getName(), "force(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.force().index, 5) + self.assertEqual(Sofa.Core.ConstVecDerivId.force().getName(), "force(V_DERIV)") + + # externalForce + self.assertEqual(Sofa.Core.VecDerivId.externalForce().index, 6) + self.assertEqual(Sofa.Core.VecDerivId.externalForce().getName(), "externalForce(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.externalForce().index, 6) + self.assertEqual(Sofa.Core.ConstVecDerivId.externalForce().getName(), "externalForce(V_DERIV)") + + # dx + self.assertEqual(Sofa.Core.VecDerivId.dx().index, 7) + self.assertEqual(Sofa.Core.VecDerivId.dx().getName(), "dx(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.dx().index, 7) + self.assertEqual(Sofa.Core.ConstVecDerivId.dx().getName(), "dx(V_DERIV)") + + # dforce + self.assertEqual(Sofa.Core.VecDerivId.dforce().index, 8) + self.assertEqual(Sofa.Core.VecDerivId.dforce().getName(), "dforce(V_DERIV)") + self.assertEqual(Sofa.Core.ConstVecDerivId.dforce().index, 8) + self.assertEqual(Sofa.Core.ConstVecDerivId.dforce().getName(), "dforce(V_DERIV)") + + # Dynamic offset + self.assertEqual(Sofa.Core.VecId.getFirstDynamicIndex(Sofa.Core.VecType.V_DERIV), 9) + self.assertEqual(Sofa.Core.ConstVecId.getFirstDynamicIndex(Sofa.Core.VecType.V_DERIV), 9) + + def test_default_matderiv_vectors(self): + # constraintJacobian + self.assertEqual(Sofa.Core.MatrixDerivId.constraintJacobian().index, 1) + self.assertEqual(Sofa.Core.MatrixDerivId.constraintJacobian().getName(), "holonomic(V_MATDERIV)") + self.assertEqual(Sofa.Core.ConstMatrixDerivId.constraintJacobian().index, 1) + self.assertEqual(Sofa.Core.ConstMatrixDerivId.constraintJacobian().getName(), "holonomic(V_MATDERIV)") + + # mappingJacobian + self.assertEqual(Sofa.Core.MatrixDerivId.mappingJacobian().index, 2) + self.assertEqual(Sofa.Core.MatrixDerivId.mappingJacobian().getName(), "nonHolonomic(V_MATDERIV)") + self.assertEqual(Sofa.Core.ConstMatrixDerivId.mappingJacobian().index, 2) + self.assertEqual(Sofa.Core.ConstMatrixDerivId.mappingJacobian().getName(), "nonHolonomic(V_MATDERIV)") + + # Dynamic offset + self.assertEqual(Sofa.Core.VecId.getFirstDynamicIndex(Sofa.Core.VecType.V_MATDERIV), 3) + self.assertEqual(Sofa.Core.ConstVecId.getFirstDynamicIndex(Sofa.Core.VecType.V_MATDERIV), 3) diff --git a/examples/StaticSolver.py b/examples/StaticSolver.py new file mode 100644 index 00000000..25b95682 --- /dev/null +++ b/examples/StaticSolver.py @@ -0,0 +1,69 @@ +import Sofa +from Sofa import SofaSimulationCore +import numpy as np + +class StaticOdeSolver (Sofa.Core.OdeSolver): + def solve(self, _, dt, X, V): + mparams = Sofa.Core.MechanicalParams() + vop = SofaSimulationCore.VectorOperations(self.getContext()) + mop = SofaSimulationCore.MechanicalOperations(mparams, self.getContext()) + + # Allocate the solution vector + dx = Sofa.Core.VecId.dx() + vop.v_realloc(dx, interactionForceField=True, propagate=True) + vop.v_clear(dx) + + # Compute the residual + F = Sofa.Core.VecId.force() + mop.computeForce(result=F, clear=True, accumulate=True); + mop.projectResponse(F) + vop.v_dot(F, F) + F_norm = np.sqrt(vop.finish()) + + # Assemble the system + mop.m_setSystemMBKMatrix(m=0, b=0, k=-1) + mop.m_setSystemRHVector(F) + mop.m_setSystemLHVector(dx) + + # Solve the system + mop.m_solveSystem() + vop.v_dot(dx, dx) + dx_norm = np.sqrt(vop.finish()) + + # Propagate the solution + vop.v_peq(X, dx) # X += dx + + # Solve the constraints + # todo: Bind ConstraintParams + # mop.solveConstraint(X, ConstraintParams.ConstOrder.POS) + + print(f"Solved with |F| = {F_norm} and |dx| = {dx_norm}") + +def createScene(root): + w, l = 2, 10 + nx, ny = 3, 9 + dx, dy = w/nx/2, l/ny/2 + root.addObject('RequiredPlugin', pluginName='SofaLoader SofaBoundaryCondition SofaEngine SofaSimpleFem SofaImplicitOdeSolver SofaOpenglVisual SofaSparseSolver SofaTopologyMapping') + root.addObject('MeshObjLoader', name='surface', filename='mesh/cylinder.obj') + root.addObject('SparseGridTopology', src='@surface', name='grid', n=[nx, ny, nx]) + root.addChild('mechanics') +# root.mechanics.addObject('StaticSolver', name='solver', printLog=False) + root.mechanics.addObject(StaticOdeSolver(name='solver')) + root.mechanics.addObject('SparseLDLSolver') + root.mechanics.addObject('MechanicalObject', name='mo', src='@../grid') + root.mechanics.addObject('HexahedronSetTopologyContainer', name='hexa_topology', hexahedra='@../grid.hexahedra') + root.mechanics.addObject('HexahedronFEMForceField', youngModulus=3000, poissonRatio=0.3) + root.mechanics.addObject('BoxROI', name='fixed_roi', box=[-1-dx,0-dx,-1-dx, 1+dx, 0+dx, 1+dx], drawBoxes=True) + root.mechanics.addObject('FixedConstraint', indices='@fixed_roi.indices') + root.mechanics.addChild('visual') + root.mechanics.visual.addObject('OglModel', name='vm', src='@/surface') + root.mechanics.visual.addObject('BarycentricMapping', applyRestPosition=True) + root.mechanics.addChild('neumann') + root.mechanics.neumann.addObject('MechanicalObject', name='mo') + root.mechanics.neumann.addObject('QuadSetTopologyContainer', name='quad_topology') + root.mechanics.neumann.addObject('QuadSetTopologyModifier') + root.mechanics.neumann.addObject('QuadSetGeometryAlgorithms') + root.mechanics.neumann.addObject('Hexa2QuadTopologicalMapping', input='@../hexa_topology', output='@quad_topology') + root.mechanics.neumann.addObject('SubsetMapping', applyRestPosition=True, input='@../mo', output='@./mo', indices='@quad_topology.points') + root.mechanics.neumann.addObject('BoxROI', name='top_roi', quad='@quad_topology.quads', src='@mo', box=[-1-dx,10-dx,-1-dx, 1+dx, 10+dx, 1+dx], drawBoxes=True) + root.mechanics.neumann.addObject('QuadPressureForceField', pressure=[0, 0, -1], quadList='@top_roi.quadIndices')