Bindings

CMakeLists.txt

@@ -77,3 +77,8 @@ endif()
 
 sofa_add_subdirectory(plugin extensions/CUDA Elasticity.CUDA)
 sofa_add_subdirectory(application benchmark/cpp ElasticityBenchmark)
+
+sofa_find_package(SofaPython3 QUIET)
+if(SofaPython3_FOUND)
+    add_subdirectory(extensions/Python)
+endif()

examples/python/stvenant_kirchhoff.py

@@ -0,0 +1,117 @@
+import Sofa
+import Elasticity
+import Sofa.SofaDeformable
+import SofaRuntime
+from SofaTypes import Mat3x3
+
+def kronecker(i, j):
+    return 1.0 if i == j else 0.0
+
+
+class Material(Elasticity.HyperElasticMaterialVec3d):
+    def __init__(self, young_modulus : float, poisson_ratio : float, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+        self.addData(name="young_modulus", value=young_modulus, type="float", help="Young's modulus", group="Material properties")
+        self.addData(name="poisson_ratio", value=poisson_ratio, type="float", help="Poisson's ratio", group="Material properties")
+
+    def secondPiolaKirchhoffStress(self, deformation_gradient):
+        F = deformation_gradient
+        C = F.transposed() * F
+        E = (1.0/2.0) * (C - Mat3x3.Identity())
+
+        lambda_param, mu_param = Sofa.SofaDeformable.toLameParameters3D(self.young_modulus.value,
+                                                                        self.poisson_ratio.value)
+
+        return lambda_param * E.trace() * Mat3x3.Identity() + 2.0 * mu_param * E
+
+    def firstPiolaKirchhoffStress(self, deformation_gradient):
+        return deformation_gradient * self.secondPiolaKirchhoffStress(deformation_gradient)
+
+
+    def materialTangentModulus(self, deformation_gradient):
+
+        F = deformation_gradient
+        S = self.secondPiolaKirchhoffStress(deformation_gradient)
+
+        lambda_param, mu_param = Sofa.SofaDeformable.toLameParameters3D(self.young_modulus.value,
+                                                                        self.poisson_ratio.value)
+
+        def elasticity_tensor(i, j, k, l):
+            return mu_param * (kronecker(i, k) * kronecker(j, l) + kronecker(i, l) * kronecker(j, k)) + \
+                lambda_param * kronecker(i, j) * kronecker(k, l)
+
+        def tensor_generation(i, j, k, l):
+            A_ijkl = kronecker(i, k) * S[l][j]
+            for q in range(3):
+                for r in range(3):
+                    A_ijkl = A_ijkl + F[i][q] * elasticity_tensor(q, j, l, r) * F[k][r]
+            return A_ijkl
+
+        return tensor_generation
+
+
+def createScene(root_node):
+
+    with root_node.addChild("plugins") as plugins:
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Constraint.Projective")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Engine.Select")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.LinearSolver.Direct")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.LinearSolver.Iterative")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.LinearSolver.Ordering")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.LinearSystem")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Mass")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.ODESolver.Backward")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.SolidMechanics.FEM.Elastic")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.StateContainer")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Topology.Container.Dynamic")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Topology.Container.Grid")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Topology.Mapping")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.Component.Visual")
+        plugins.addObject('RequiredPlugin', pluginName="Sofa.GL.Component.Rendering3D")
+        plugins.addObject('RequiredPlugin', pluginName="Elasticity")
+
+    root_node.addObject('DefaultAnimationLoop')
+
+    root_node.addObject('VisualStyle', displayFlags="showBehaviorModels showForceFields")
+    root_node.addObject('VisualGrid', size="0.1")
+    root_node.addObject('LineAxis', size="0.1")
+    root_node.addObject('OglSceneFrame')
+
+    # root_node.addObject('EulerImplicitSolver', name="backward_Euler", rayleighStiffness="0.1", rayleighMass="0.1")
+    root_node.addObject('NewtonRaphsonSolver')
+    root_node.addObject('StaticSolver', name="solver")
+    root_node.addObject('SparseLDLSolver', name="linear_solver", template="CompressedRowSparseMatrixMat3x3d")
+
+    root_node.addObject('RegularGridTopology', name="grid", min="-0.01 -0.01 0", max="0.01 0.01 0.2", n="5 5 30")
+    root_node.addObject('MechanicalObject', template="Vec3", name="state", showObject="true")
+    root_node.addObject('NodalMassDensity', property="1100")
+    root_node.addObject('HexahedronFEMMass')
+    root_node.addObject('BoxROI', template="Vec3", name="box_roi", box="-0.011 -0.011 -0.0001   0.011 0.011 0.0001",
+                        drawBoxes="1")
+    root_node.addObject('FixedProjectiveConstraint', template="Vec3", indices="@box_roi.indices")
+    root_node.addObject('HyperelasticityFEMForceField', template="Vec3,Hexahedron", name="fem",
+                        computeForceStrategy="sequenced", computeForceDerivStrategy="sequenced")
+
+    young_modulus = 2e5
+    poisson_ratio = 0.45
+    material = Material(young_modulus, poisson_ratio)
+    root_node.addObject(material)
+
+def main():
+    import SofaRuntime
+    import SofaImGui
+    import Sofa.Gui
+
+    root = Sofa.Core.Node("root")
+    createScene(root)
+    Sofa.Simulation.initRoot(root)
+
+    Sofa.Gui.GUIManager.Init("myscene", "imgui")
+    Sofa.Gui.GUIManager.createGUI(root, __file__)
+    Sofa.Gui.GUIManager.SetDimension(1600, 900)
+    Sofa.Gui.GUIManager.MainLoop(root)
+    Sofa.Gui.GUIManager.closeGUI()
+
+if __name__ == "__main__":
+    main()

extensions/Python/CMakeLists.txt

@@ -0,0 +1,23 @@
+project(Bindings_Elasticity)
+
+set(SOURCE_FILES
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Module_Elasticity.cpp
+
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Binding_HyperelasticMaterial.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Binding_HyperelasticMaterial.cpp
+)
+
+if (NOT TARGET SofaPython3::Plugin)
+    find_package(SofaPython3 REQUIRED COMPONENTS Elasticity Bindings.Sofa)
+endif()
+
+SP3_add_python_module(
+    TARGET       ${PROJECT_NAME}
+    PACKAGE      Elasticity
+    MODULE       Elasticity
+    DESTINATION  .
+    SOURCES      ${SOURCE_FILES}
+    DEPENDS      SofaPython3::Plugin SofaPython3::Bindings.Sofa Elasticity
+
+)
+message("-- SofaPython3 bindings for Elasticity will be created.")

extensions/Python/src/Binding_HyperelasticMaterial.cpp

@@ -0,0 +1,137 @@
+#include <Binding_HyperelasticMaterial.h>
+#include <Elasticity/component/HyperelasticMaterial.h>
+#include <SofaPython3/PythonEnvironment.h>
+#include <SofaPython3/PythonFactory.h>
+#include <SofaPython3/Sofa/Core/Binding_Base.h>
+#include <string>
+
+namespace elasticity::python3
+{
+namespace py { using namespace pybind11; }
+using namespace std::string_literals;
+
+/**
+ * Trampoline class for HyperelasticMaterial
+ */
+template <class TDataTypes>
+class PyHyperelasticMaterial : public HyperelasticMaterial<TDataTypes>
+{
+public:
+    SOFA_CLASS(PyHyperelasticMaterial, HyperelasticMaterial<TDataTypes>);
+    using HyperelasticMaterial<TDataTypes>::HyperelasticMaterial;
+    using StressTensor = HyperelasticMaterial<TDataTypes>::StressTensor;
+    using TangentModulus = HyperelasticMaterial<TDataTypes>::TangentModulus;
+
+    void init() override;
+    StressTensor firstPiolaKirchhoffStress(Strain<TDataTypes>& strain) override;
+    TangentModulus materialTangentModulus(Strain<TDataTypes>& strain) override;
+};
+
+template <class TDataTypes>
+void PyHyperelasticMaterial<TDataTypes>::init()
+{
+    HyperelasticMaterial<TDataTypes>::init();
+    sofapython3::PythonEnvironment::gil acquire;
+    PYBIND11_OVERRIDE(void, HyperelasticMaterial<TDataTypes>, init, );
+}
+
+template <class TDataTypes>
+auto PyHyperelasticMaterial<TDataTypes>::firstPiolaKirchhoffStress(Strain<TDataTypes>& strain) -> StressTensor
+{
+    const auto& F = strain.deformationGradient();
+    PYBIND11_OVERLOAD_PURE(StressTensor, HyperelasticMaterial<TDataTypes>, firstPiolaKirchhoffStress, F);
+}
+
+template <class TDataTypes>
+auto PyHyperelasticMaterial<TDataTypes>::materialTangentModulus(Strain<TDataTypes>& strain) -> TangentModulus
+{
+    using Real = sofa::Real_t<TDataTypes>;
+    using Callable = std::function<Real(sofa::Size, sofa::Size, sofa::Size, sofa::Size)>;
+
+    sofapython3::PythonEnvironment::gil acquire;
+
+    py::function override = py::get_override(
+        static_cast<const HyperelasticMaterial<TDataTypes>*>(this),
+        "materialTangentModulus");
+
+    if (!override)
+    {
+        pybind11::pybind11_fail(
+            "Tried to call pure virtual function "
+            "\"HyperelasticMaterial::materialTangentModulus\"");
+    }
+
+    const auto& F = strain.deformationGradient();
+    py::object pythonCallable = override(F);
+
+    if (!PyCallable_Check(pythonCallable.ptr()))
+    {
+        throw py::type_error("materialTangentModulus must return a callable accepting 4 parameters");
+    }
+
+    auto callable = [pythonCallable](sofa::Size i, sofa::Size j, sofa::Size k, sofa::Size l) -> Real
+    {
+        return pythonCallable(i, j, k, l).template cast<Real>();
+    };
+
+    return TangentModulus(callable);
+}
+
+template <class TDataTypes>
+void declareHyperElasticMaterial(pybind11::module &m)
+{
+    const std::string pyclassName = "HyperElasticMaterial"s + TDataTypes::Name();
+
+    using Class = HyperelasticMaterial<TDataTypes>;
+    using Trampoline = PyHyperelasticMaterial<TDataTypes>;
+    using Base = sofa::core::objectmodel::BaseComponent;
+
+    py::class_<
+        Class,
+        Base,
+        Trampoline,
+        sofapython3::py_shared_ptr<Class>
+    > c(m, pyclassName.c_str(), py::dynamic_attr());
+
+    // Python constructor
+    c.def(py::init([](py::args &args, py::kwargs &kwargs)
+    {
+        auto material = sofa::core::sptr<PyHyperelasticMaterial<TDataTypes>> (new PyHyperelasticMaterial<TDataTypes>());
+        material->f_listening.setValue(true);
+
+        if (args.size() == 1)
+            material->setName(py::cast<std::string>(args[0]));
+
+        py::object cc = py::cast(material);
+        for (auto kv : kwargs)
+        {
+            std::string key = py::cast<std::string>(kv.first);
+            py::object value = py::reinterpret_borrow<py::object>(kv.second);
+            if (key == "name")
+            {
+                if (!args.empty())
+                {
+                    throw py::type_error("The name is set twice as a "
+                                        "named argument='" + py::cast<std::string>(value) + "' and as a"
+                                                                                            "positional argument='" +
+                                        py::cast<std::string>(args[0]) + "'.");
+                }
+            }
+            sofapython3::BindingBase::SetAttr(cc, key, value);
+        }
+        return material;
+    }));
+
+    sofapython3::PythonFactory::registerType<HyperelasticMaterial<TDataTypes>>([](sofa::core::objectmodel::Base* object)
+    {
+        return py::cast(dynamic_cast<HyperelasticMaterial<TDataTypes>*>(object));
+    });
+}
+
+void moduleAddHyperelasticMaterial(pybind11::module& m)
+{
+    declareHyperElasticMaterial<sofa::defaulttype::Vec2Types>(m);
+    declareHyperElasticMaterial<sofa::defaulttype::Vec3Types>(m);
+}
+
+}  // namespace elasticity::python3

extensions/Python/src/Binding_HyperelasticMaterial.h

@@ -0,0 +1,7 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+namespace elasticity::python3
+{
+void moduleAddHyperelasticMaterial(pybind11::module &m);
+}

extensions/Python/src/Module_Elasticity.cpp

@@ -0,0 +1,38 @@
+/******************************************************************************
+*                 SOFA, Simulation Open-Framework Architecture                *
+*                    (c) 2021 INRIA, USTL, UJF, CNRS, MGH                     *
+*                                                                             *
+* 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: contact@sofa-framework.org                             *
+******************************************************************************/
+
+#include <pybind11/pybind11.h>
+
+#include <Elasticity/init.h>
+#include <Binding_HyperelasticMaterial.h>
+
+
+namespace py { using namespace pybind11; }
+
+namespace elasticity::python3
+{
+
+PYBIND11_MODULE(Elasticity, m)
+{
+    elasticity::initializePlugin();
+    moduleAddHyperelasticMaterial(m);
+}
+
+} // namespace elasticity::python3

src/Elasticity/component/HyperelasticityFEMForceField.h

@@ -84,9 +84,6 @@ class HyperelasticityFEMForceField :
     DeformationGradient computeDeformationGradient(
         const sofa::type::Mat<spatial_dimensions, TopologicalDimension, Real>& J_q,
         const sofa::type::Mat<TopologicalDimension, spatial_dimensions, Real>& J_Q_inv);
-    DeformationGradient computeDeformationGradient2(
-        const std::array<Coord, NumberOfNodesInElement>& elementNodesCoordinates,
-        const sofa::type::Mat<NumberOfNodesInElement, spatial_dimensions, Real>& dN_dQ);
 
     struct PrecomputedData
     {

src/Elasticity/component/HyperelasticityFEMForceField.inl

@@ -181,9 +181,7 @@ void HyperelasticityFEMForceField<DataTypes, ElementType>::computeHessian(const
             // gradient of the shape functions in the physical element evaluated at the quadrature point
             const sofa::type::Mat<NumberOfNodesInElement, spatial_dimensions, Real>& dN_dQ = precomputedData.dN_dQ;
 
-            // both ways to compute the deformation gradient are equivalent
             const DeformationGradient F = computeDeformationGradient(J_q, J_Q_inv);
-            // const DeformationGradient F = computeDeformationGradient2(elementNodesCoordinates, dN_dQ);
 
             Strain<DataTypes> strain(deformationGradient, F);
 
@@ -232,19 +230,6 @@ auto HyperelasticityFEMForceField<DataTypes, ElementType>::computeDeformationGra
     return J_q * J_Q_inv;
 }
 
-template <class DataTypes, class ElementType>
-auto HyperelasticityFEMForceField<DataTypes, ElementType>::computeDeformationGradient2(
-    const std::array<Coord, NumberOfNodesInElement>& elementNodesCoordinates,
-    const sofa::type::Mat<NumberOfNodesInElement, spatial_dimensions, Real>& dN_dQ)  -> DeformationGradient
-{
-    DeformationGradient F;
-
-    for (sofa::Size i = 0; i < NumberOfNodesInElement; ++i)
-        F += sofa::type::dyad(elementNodesCoordinates[i], dN_dQ[i]);
-
-    return F;
-}
-
 template <class TDataTypes, class TElementType>
 void HyperelasticityFEMForceField<TDataTypes, TElementType>::precomputeData()
 {
@@ -285,6 +270,7 @@ void HyperelasticityFEMForceField<TDataTypes, TElementType>::beforeElementForce(
     const sofa::core::MechanicalParams* mparams, sofa::type::vector<ElementGradient>& f,
     const sofa::VecCoord_t<DataTypes>& x)
 {
+    //store coordinates to use it later when computing the Hessian
     m_coordinates = &x;
 
     const auto& elements = FiniteElement::getElementSequence(*this->l_topology);
@@ -340,9 +326,7 @@ void HyperelasticityFEMForceField<TDataTypes, TElementType>::computeElementsForc
             // gradient of the shape functions in the physical element evaluated at the quadrature point
             const sofa::type::Mat<NumberOfNodesInElement, spatial_dimensions, Real>& dN_dQ = precomputedData.dN_dQ;
 
-            // both ways to compute the deformation gradient are equivalent
             const DeformationGradient F = computeDeformationGradient(J_q, J_Q_inv);
-            // const DeformationGradient F = computeDeformationGradient2(elementNodesCoordinates, dN_dQ);
 
             Strain<DataTypes> strain(deformationGradient, F);
             const auto P = l_material->firstPiolaKirchhoffStress(strain);