Initial commit

CMakeLists.txt

@@ -0,0 +1,19 @@
+cmake_minimum_required(VERSION 3.15)
+project(SofaSimpleForcefield VERSION 1.0)
+
+find_package(Sofa.Core REQUIRED)
+find_package(Eigen3 REQUIRED)
+
+set(CMAKE_CXX_STANDARD 17)
+
+add_library(${PROJECT_NAME} SHARED init.cpp SVKElasticForcefield.cpp)
+target_link_libraries(${PROJECT_NAME} PUBLIC Sofa.Core Eigen3::Eigen)
+
+sofa_create_package_with_targets(
+    PACKAGE_NAME ${PROJECT_NAME}
+    PACKAGE_VERSION ${PROJECT_VERSION}
+    TARGETS ${PROJECT_NAME} AUTO_SET_TARGET_PROPERTIES
+    INCLUDE_SOURCE_DIR "src"
+    INCLUDE_INSTALL_DIR ${PROJECT_NAME}
+    RELOCATABLE "plugins"
+)

README.md

@@ -0,0 +1,38 @@
+# Simple implementation of a Saint-Venant-Kirchhoff force field
+
+This plugin can be used as a simple demonstration of how SOFA 
+"force fields" components are built up. In this case, the force
+field API will provide the internal potential energy,
+residual and jacobian of a Saint-Venant-Kirchhoff material assembled
+on linear tetrahedral meshes. These three function will be used 
+by the system assembler and ODE solvers of SOFA.
+
+## Compiling with Ubuntu LTS 20.04
+### SOFA (skip this if you already compiled SOFA)
+
+```console
+user@host:~$ sudo apt install qtbase5-dev libqt5charts5-dev libqt5opengl5-dev libopengl0 libeigen3-dev libglew-dev zlib1g-dev libboost-dev libboost-filesystem-dev g++ cmake git
+user@host:~$ export SOFA_SRC=/opt/sofa_src
+user@host:~$ export SOFA_ROOT=/opt/sofa
+user@host:~$ git clone https://github.com/sofa-framework/sofa.git $SOFA_SRC
+user@host:~$ cmake -DCMAKE_INSTALL_PREFIX=$SOFA_ROOT -DCMAKE_BUILD_TYPE=Release -S $SOFA_SRC -B $SOFA_SRC/build
+user@host:~$ cmake --build $SOFA_SRC/build -j4
+user@host:~$ cmake --install $SOFA_SRC/build
+```
+
+### SofaSimpleForcefield
+```console
+user@host:~$ export SSFF_SRC=/opt/SofaSimpleForceField_src
+user@host:~$ git clone https://github.com/jnbrunet/SofaSimpleForcefield.git $SSFF_SRC
+user@host:~$ cmake -DCMAKE_PREFIX_PATH=$SOFA_ROOT/lib/cmake -DCMAKE_INSTALL_PREFIX=$SOFA_ROOT/plugins/SofaSimpleForceField -DCMAKE_BUILD_TYPE=Release -S $SSFF_SRC -B $SSFF_SRC/build
+user@host:~$ cmake --build $SSFF_SRC/build -j4
+user@host:~$ cmake --install $SSFF_SRC/build
+```
+
+## Running the cantilever example scene
+```console
+user@host:~$ $SOFA_ROOT/bin/runSofa -l SofaSimpleForceField $SSFF_SRC/cantilever_beam.scn
+```
+
+## Result
+![image](https://user-images.githubusercontent.com/6951981/127413110-76fb452e-723b-4e74-b3ac-5952d54d663d.png)

SVKElasticForcefield.h

@@ -0,0 +1,75 @@
+#pragma once
+
+#include <sofa/core/behavior/ForceField.h>
+#include <sofa/core/topology/BaseTopology.h>
+#include <Eigen/Eigen>
+#include "Tetrahedron.h"
+
+
+class SVKElasticForcefield : public sofa::core::behavior::ForceField<sofa::defaulttype::Vec3Types> {
+public:
+    SOFA_CLASS(SVKElasticForcefield, SOFA_TEMPLATE(sofa::core::behavior::ForceField, sofa::defaulttype::Vec3Types));
+
+    // Aliases
+    using Element = Tetrahedron;
+    using Coord = sofa::type::Vec3;
+
+    template <class T>
+    using Data = sofa::core::objectmodel::Data<T>;
+
+    template <typename ObjectType>
+    using Link = sofa::core::objectmodel::SingleLink<SVKElasticForcefield, ObjectType, sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK>;
+
+    // Data structures
+
+    struct GaussNode {
+        Real weight = 0;
+        Real jacobian_determinant = 0;
+        Eigen::Matrix<double, Element::NumberOfNodes, 3> dN_dx = Eigen::Matrix<double, Element::NumberOfNodes, 3>::Zero();
+    };
+
+    // public methods
+    SVKElasticForcefield();
+
+    /**
+     * Initialize the forcefield by pre-computing the derivative of the shape function at each
+     * Gauss point w.r.t the initial (undeformed) position of the mesh.
+     */
+    void init() override;
+
+    /** StvK elasticity potential, i.e. W(x) = 1/2 lambda Tr^2(e) + mu e:e */
+    double getPotentialEnergy (
+            const sofa::core::MechanicalParams * mparams,
+            const Data<sofa::type::vector<sofa::type::Vec3>> & d_x) const override;
+
+    /** Linear elasticity residual, i.e. R(x) = dW(x)/dE */
+    void addForce(
+            const sofa::core::MechanicalParams * mparams,
+            Data<sofa::type::vector<sofa::type::Vec3>> & d_f,
+            const Data<sofa::type::vector<sofa::type::Vec3>>& d_x,
+            const Data<sofa::type::vector<sofa::type::Vec3>>& d_v) override;
+
+    /** Jacobian of the elastic residual, i.e. K(x) = dR(x)/dE */
+    void addKToMatrix(
+            sofa::defaulttype::BaseMatrix * matrix,
+            double kFact,
+            unsigned int & offset) override;
+
+    /** Performs df = K*dx where K is the tangent stiffness matrix (see addKToMatrix) */
+    void addDForce(
+            const sofa::core::MechanicalParams* mparams,
+            Data<sofa::type::vector<sofa::type::Vec3>> & d_df,
+            const Data<sofa::type::vector<sofa::type::Vec3>> & d_dx) override;
+
+    /** Draw the computational mesh as the simulation advance. Only for debug purposes. */
+    void draw(const sofa::core::visual::VisualParams* vparams) override;
+
+    /** Used to automatically define a minimal view box for the GUI */
+    void computeBBox(const sofa::core::ExecParams* params, bool onlyVisible) override;
+
+private:
+    Data<double> d_youngModulus;
+    Data<double> d_poissonRatio;
+    Link<sofa::core::topology::BaseMeshTopology>   d_topology_container;          ///< Container of node indices per element
+    std::vector<std::array<GaussNode, Element::NumberOfGaussNode>> p_gauss_nodes; ///< Set of Gauss nodes per elements
+};

SofaSimpleForcefieldConfig.cmake.in

@@ -0,0 +1,13 @@
+# CMake package configuration file for the plugin @PROJECT_NAME@
+
+@PACKAGE_GUARD@
+@PACKAGE_INIT@
+
+find_package(Sofa.Core REQUIRED)
+find_package(Eigen3 REQUIRED)
+
+if(NOT TARGET @PROJECT_NAME@)
+include("${CMAKE_CURRENT_LIST_DIR}/@[email protected]")
+endif()
+
+check_required_components(@PROJECT_NAME@)

Tetrahedron.h

@@ -0,0 +1,43 @@
+#pragma once
+
+#include <Eigen/Eigen>
+#include <array>
+
+struct Tetrahedron {
+    static constexpr auto NumberOfNodes = 4;
+    static constexpr auto NumberOfGaussNode = 1;
+
+    struct GaussNode {
+        Eigen::Vector3d position;
+        double weight;
+    };
+
+    /// Shape function of each nodes evaluated at local coordinates (u, v, w)
+    static inline auto N(double u, double v, double w) -> Eigen::Matrix<double, NumberOfNodes, 1> {
+        return {
+                1 - u - v - w,  // Node 0
+                u,              // Node 1
+                v,              // Node 1
+                w               // Node 3
+        };
+    }
+
+    /// Derivatives of the shape function at each nodes w.r.t local coordinates and evaluated at (u, v, w)
+    static inline auto dN(double /*u*/, double /*v*/, double /*w*/) -> Eigen::Matrix<double, NumberOfNodes, 3> {
+        Eigen::Matrix<double, NumberOfNodes, 3> m;
+        //  dL/du  dL/dv  dL/dw
+        m <<   -1,    -1,    -1,   // Node 0
+                1,     0,     0,   // Node 1
+                0,     1,     0,   // Node 1
+                0,     0,     1;   // Node 3
+        return m;
+    }
+
+    /// Numerical integration points of the element
+    static inline auto gauss_nodes() -> const std::array<GaussNode, NumberOfGaussNode> & {
+        static const std::array<GaussNode, NumberOfGaussNode> gauss_nodes {
+                GaussNode {Eigen::Vector3d(1/4., 1/4., 1/4.), double(1/6.)} // Node 0
+        };
+        return gauss_nodes;
+    }
+};

cantilever_beam.scn

@@ -0,0 +1,23 @@
+<Node>
+    <RequiredPlugin pluginName="SofaSimpleForcefield SofaBaseMechanics SofaBaseTopology SofaSparseSolver SofaImplicitOdeSolver SofaTopologyMapping SofaBoundaryCondition SofaEngine" />
+    <VisualStyle displayFlags="showForceFields showBehaviorModels" />
+    <RegularGridTopology name="grid" min="-7.5 -7.5 0" max="7.5 7.5 80" n="9 9 21" />
+
+    <StaticSolver newton_iterations="25" relative_correction_tolerance_threshold="1e-15" relative_residual_tolerance_threshold="1e-10" printLog="1" />
+    <SparseLDLSolver template="CompressedRowSparseMatrixMat3x3d"/>
+
+    <MechanicalObject name="mo" src="@grid" />
+    <TetrahedronSetTopologyContainer name="topology" />
+    <TetrahedronSetTopologyModifier/>
+    <Hexa2TetraTopologicalMapping input="@grid" output="@topology" swapping="1" />
+
+    <SVKElasticForcefield youngModulus="3000" poissonRatio="0.3" topology="@topology" />
+
+    <BoxROI name="fixed_roi" box="-7.5 -7.5 -0.9 7.5 7.5 0.1" />
+    <FixedConstraint indices="@fixed_roi.indices" />
+
+    <BoxROI name="top_roi" box="-7.5 -7.5 79.9 7.5 7.5 80.1" />
+    <TriangleSetGeometryAlgorithms />
+    <TrianglePressureForceField pressure="0 -10 0" topology="@topology" triangleList="@top_roi.triangleIndices" showForces="1" />
+
+</Node>

init.cpp

@@ -0,0 +1,44 @@
+// Here are just several convenient functions to help user to know what contains the plugin
+
+extern "C" {
+void        initExternalModule();
+const char* getModuleName();
+const char* getModuleVersion();
+const char* getModuleLicense();
+const char* getModuleDescription();
+const char* getModuleComponentList();
+}
+
+void initExternalModule()
+{
+    static bool first = true;
+    if (first)
+    {
+        first = false;
+    }
+}
+
+const char* getModuleName()
+{
+    return "SofaSimpleForcefield";
+}
+
+const char* getModuleVersion()
+{
+    return "1.0";
+}
+
+const char* getModuleLicense()
+{
+    return "LGPL 2.1";
+}
+
+const char* getModuleDescription()
+{
+    return "Simple implementation of a Saint-Venant-Kirchhoff force field";
+}
+
+const char* getModuleComponentList()
+{
+    return "SVKElasticForcefield";
+}