COMP: Clears compilation errors in force model

Clears compilation errors in force model. Implements getForceAndMatrix in force models.
Adds get type in dynamical model. Adds get and set dynamic model in dynamic object.
Removes raw pointers in nonLinearSystem class.

Base/ForceModel/CMakeLists.txt

@@ -11,6 +11,7 @@ imstk_add_library( ForceModel
 	VegaFEM::forceModel
     VegaFEM::stvk
     VegaFEM::graph
+	VegaFEM::volumetricMesh
   )
 
 #-----------------------------------------------------------------------------

Base/ForceModel/imstkCorotationalFEMForceModel.h

@@ -37,11 +37,14 @@ namespace imstk
 class CorotationalFEMForceModel : virtual public InternalForceModel
 {
 public:
-    CorotationalFEMForceModel(std::shared_ptr<vega::VolumetricMesh> mesh, int warp = 1)
+    CorotationalFEMForceModel(std::shared_ptr<vega::VolumetricMesh> mesh, int warp = 1) : InternalForceModel()
     {
-        m_corotationalLinearFEM = std::make_shared<vega::CorotationalLinearFEM>(mesh.get());
+        auto tetMesh = std::dynamic_pointer_cast<vega::TetMesh>(mesh);
+
+        m_corotationalLinearFEM = std::make_shared<vega::CorotationalLinearFEM>(tetMesh.get());
         m_warp = warp;
     }
+    CorotationalFEMForceModel() = delete;
 
     virtual ~CorotationalFEMForceModel();
 

Base/ForceModel/imstkForceModelConfig.cpp

@@ -121,6 +121,7 @@ ForceModelConfig::getForceModelType()
 HyperElasticMaterialType
 ForceModelConfig::getHyperelasticMaterialType()
 {
+    // Set up some variables
     if (this->m_stringsOptionMap["invertibleMaterial"] == "StVK")
     {
         return HyperElasticMaterialType::StVK;
@@ -135,7 +136,7 @@ ForceModelConfig::getHyperelasticMaterialType()
     }
     else
     {
-        LOG(INFO) << "Hyperelastic model type not assigned";
+        LOG(INFO) << "Force model type not assigned";
         return HyperElasticMaterialType::none;
     }
 }

Base/ForceModel/imstkIsotropicHyperelasticFEMForceModel.h

@@ -45,27 +45,29 @@ public:
         bool withGravity = true,
         double gravity = 10.0) : InternalForceModel()
     {
+        auto tetMesh = std::dynamic_pointer_cast<vega::TetMesh>(mesh);
+
         int enableCompressionResistance = 1;
         double compressionResistance = 500;
         switch (materialType)
         {
             case HyperElasticMaterialType::StVK:
                 m_isotropicMaterial = std::make_shared<vega::StVKIsotropicMaterial>(
-                    mesh.get(),
+                    tetMesh.get(),
                     enableCompressionResistance,
                     compressionResistance);
                 break;
 
             case HyperElasticMaterialType::NeoHookean:
                 m_isotropicMaterial = std::make_shared<vega::NeoHookeanIsotropicMaterial>(
-                    mesh.get(),
+                    tetMesh.get(),
                     enableCompressionResistance,
                     compressionResistance);
                 break;
 
             case HyperElasticMaterialType::MooneyRivlin:
                 m_isotropicMaterial = std::make_shared<vega::MooneyRivlinIsotropicMaterial>(
-                    mesh.get(),
+                    tetMesh.get(),
                     enableCompressionResistance,
                     compressionResistance);
                 break;
@@ -75,13 +77,15 @@ public:
         }
 
         m_isotropicHyperelasticFEM = std::make_shared<vega::IsotropicHyperelasticFEM>(
-            mesh.get(),
-            m_isotropicMaterial,
+            tetMesh.get(),
+            m_isotropicMaterial.get(),
             inversionThreshold,
             withGravity,
             gravity);
     }
 
+    IsotropicHyperelasticFEForceModel() = delete;
+
     virtual ~IsotropicHyperelasticFEForceModel();
 
     void getInternalForce(const Vectord& u, Vectord& internalForce)

Base/ForceModel/imstkLinearFEMForceModel.h

@@ -40,7 +40,9 @@ class LinearFEMForceModel : virtual public InternalForceModel
 public:
     LinearFEMForceModel(std::shared_ptr<vega::VolumetricMesh> mesh, bool withGravity = true, double gravity = -9.81) : InternalForceModel()
     {
-        m_stVKInternalForces = std::make_shared<vega::StVKInternalForces>(mesh.get(), 0, withGravity, gravity);
+        auto tetMesh = std::dynamic_pointer_cast<vega::TetMesh>(mesh);
+
+        m_stVKInternalForces = std::make_shared<vega::StVKInternalForces>(tetMesh.get(), nullptr, withGravity, gravity);
 
         auto stVKStiffnessMatrix = std::make_shared<vega::StVKStiffnessMatrix>(m_stVKInternalForces.get());
         auto s = m_stiffnessMatrix.get();
@@ -59,11 +61,21 @@ public:
         m_stiffnessMatrix->MultiplyVector(data, internalForce.data());
     }
 
-    void getTangentStiffnessMatrix(const Vectord& u, SparseMatrixd tangentStiffnessMatrix)
+    void getTangentStiffnessMatrix(const Vectord& u, SparseMatrixd& tangentStiffnessMatrix)
     {
         InternalForceModel::updateValuesFromMatrix(m_stiffnessMatrix, tangentStiffnessMatrix.valuePtr());
     }
 
+    void getTangentStiffnessMatrixTopology(vega::SparseMatrix** tangentStiffnessMatrix)
+    {
+        *tangentStiffnessMatrix = new vega::SparseMatrix(*m_stiffnessMatrix.get());
+    }
+
+    void GetForceAndMatrix(const Vectord& u, Vectord& internalForce, SparseMatrixd& tangentStiffnessMatrix)
+    {
+        getInternalForce(u, internalForce);
+        getTangentStiffnessMatrix(u, tangentStiffnessMatrix);
+    }
 protected:
     std::shared_ptr<vega::SparseMatrix> m_stiffnessMatrix;
     std::shared_ptr<vega::StVKInternalForces> m_stVKInternalForces;

Base/ForceModel/imstkMassSpringForceModel.h

@@ -39,6 +39,7 @@ public:
     {
         m_massSpringSystem = massSpringSystem;
     }
+    MassSpringForceModel() = delete;
     virtual ~MassSpringForceModel();
 
     void getInternalForce(const Vectord& u, Vectord& internalForce)

Base/ForceModel/imstkStVKForceModel.h

@@ -30,6 +30,7 @@
 //vega
 #include "StVKInternalForces.h"
 #include "StVKStiffnessMatrix.h"
+#include "tetMesh.h"
 
 namespace imstk
 {
@@ -39,7 +40,9 @@ class StVKForceModel : virtual public InternalForceModel
 public:
     StVKForceModel(std::shared_ptr<vega::VolumetricMesh> mesh, bool withGravity = true, double gravity = 10.0) : InternalForceModel()
     {
-        m_stVKInternalForces = std::make_shared<vega::StVKInternalForces>(mesh.get(), 0, withGravity, gravity);
+        auto tetMesh = std::dynamic_pointer_cast<vega::TetMesh>(mesh);
+
+        m_stVKInternalForces = std::make_shared<vega::StVKInternalForces>(tetMesh.get(), nullptr, withGravity, gravity);
 
         m_vegaStVKStiffnessMatrix = std::make_shared<vega::StVKStiffnessMatrix>(m_stVKInternalForces.get());
     }
@@ -64,6 +67,12 @@ public:
         InternalForceModel::updateValuesFromMatrix(m_vegaTangentStiffnessMatrix, tangentStiffnessMatrix.valuePtr());
     }
 
+    void GetForceAndMatrix(const Vectord& u, Vectord& internalForce, SparseMatrixd& tangentStiffnessMatrix)
+    {
+        getInternalForce(u, internalForce);
+        getTangentStiffnessMatrix(u, tangentStiffnessMatrix);
+    }
+
 protected:
 
     std::shared_ptr<vega::StVKInternalForces> m_stVKInternalForces;

Base/Geometry/Reader/imstkVegaMeshReader.cpp

@@ -100,4 +100,5 @@ VegaMeshReader::getVolumeMeshFromVegaVolumeMesh(std::shared_ptr<vega::Volumetric
         return nullptr;
     }
 }
+
 }

Base/SceneElements/Objects/imstkDeformableBodyModel.cpp

@@ -21,6 +21,7 @@
 
 #include <fstream>
 
+//imstk
 #include "imstkDeformableBodyModel.h"
 
 // vega
@@ -30,7 +31,7 @@
 namespace imstk
 {
 
-DeformableBodyModel::DeformableBodyModel() : DynamicalModel(DynamicalModel::Type::elastoDynamics), m_damped(false){}
+DeformableBodyModel::DeformableBodyModel(DynamicalModel::Type type) : DynamicalModel(type), m_damped(false){}
 
 void
 DeformableBodyModel::setForceModelConfiguration(std::shared_ptr<ForceModelConfig> fmConfig)
@@ -98,10 +99,10 @@ DeformableBodyModel::initialize()
 
     m_vegaPhysicsMesh = VegaMeshReader::getVegaVolumeMeshFromVolumeMesh(std::static_pointer_cast<VolumetricMesh>(m_forceModelGeometry));
 
-    initializeForceModel();//c
-    initializeMassMatrix();//c
-    initializeDampingMatrix();//c
-    initializeTangentStiffness();//c
+    initializeForceModel();
+    initializeMassMatrix();
+    initializeDampingMatrix();
+    initializeTangentStiffness();
     loadInitialStates();
     loadBoundaryConditions();
     initializeGravityForce();
@@ -175,12 +176,17 @@ DeformableBodyModel::initializeForceModel()
 
     case ForceModelType::Corotational:
 
-        this->m_internalForceModel = std::make_shared<CorotationalFEMForceModel>(m_vegaPhysicsMesh, isGravityPresent, g);
+        this->m_internalForceModel = std::make_shared<CorotationalFEMForceModel>(m_vegaPhysicsMesh);
         break;
 
     case ForceModelType::Invertible:
 
-        this->m_internalForceModel = std::make_shared<IsotropicHyperelasticFEForceModel>(m_vegaPhysicsMesh, -MAX_D, isGravityPresent, g);
+        this->m_internalForceModel = std::make_shared<IsotropicHyperelasticFEForceModel>(
+                                        m_forceModelConfiguration->getHyperelasticMaterialType(),
+                                        m_vegaPhysicsMesh,
+                                        -MAX_D,
+                                        isGravityPresent,
+                                        g);
         break;
 
     default:
@@ -431,23 +437,25 @@ DeformableBodyModel::updateBodyStates(const Vectord& delataV)
 }
 
 NonLinearSystem::VectorFunctionType
-DeformableBodyModel::getFunction(const Vectord& q)
+DeformableBodyModel::getFunction()
 {
+    //const Vectord& q
     // Function to evaluate the nonlinear objective function given the current state
     return [&, this](const Vectord&) -> const Vectord&
     {
-        computeImplicitSystemRHS(state, newState);
+        computeImplicitSystemRHS(*m_previousState.get(), *m_currentState.get());
         return this->m_Feff;
     };
 }
 
 NonLinearSystem::MatrixFunctionType
-DeformableBodyModel::getFunctionGradient(const Vectord& q)
+DeformableBodyModel::getFunctionGradient()
 {
+    //const Vectord& q
     // Gradient of the nonlinear objective function given the current state
     return [&, this](const Vectord&) -> const SparseMatrixd&
     {
-        computeImplicitSystemLHS(state, newState);
+        computeImplicitSystemLHS(*m_previousState.get(), *m_currentState.get());
         return m_Keff;
     };
 }

Base/SceneElements/Objects/imstkDeformableBodyModel.h

@@ -61,7 +61,7 @@ public:
     ///
     /// \brief Constructor
     ///
-    DeformableBodyModel();
+    DeformableBodyModel(DynamicalModel::Type type = DynamicalModel::Type::elastoDynamics);
 
     ///
     /// \brief Destructor
@@ -93,7 +93,7 @@ public:
     std::shared_ptr<Geometry> getModelGeometry();
 
     ///
-    /// \brief Returns the tangent linear system given curent state
+    /// \brief Returns the tangent linear system given current state
     ///
     void getTangent(Vectord& q);
 
@@ -183,17 +183,16 @@ public:
     /// \brief Returns the "function" that evaluates the nonlinear function given
     /// the state vector
     ///
-    NonLinearSystem::VectorFunctionType getFunction(const Vectord& q);
+    NonLinearSystem::VectorFunctionType getFunction();
 
     ///
     /// \brief Returns the "function" that evaluates the gradient of the nonlinear
     /// function given the state vector
     ///
-    NonLinearSystem::MatrixFunctionType getFunctionGradient(const Vectord& q);
-
+    NonLinearSystem::MatrixFunctionType getFunctionGradient();
 
     ///
-    /// \brief
+    /// \brief Initialize the Eigen matrix with data inside vega sparse matrix
     ///
     static void initializeEigenMatrixFromVegaMatrix(const vega::SparseMatrix& vegaMatrix, SparseMatrixd& eigenMatrix);
 
@@ -212,7 +211,7 @@ protected:
     std::shared_ptr<ForceModelConfig>   m_forceModelConfiguration;  ///> Store the configuration here
     std::shared_ptr<Geometry>           m_forceModelGeometry;       ///> Geometry used by force model
 
-    bool m_damped;
+    bool m_damped;      ///> Viscous or structurally damped system
 
     /// Matrices typical to a elastodynamics and 2nd order analogous systems
     SparseMatrixd m_M;    ///> Mass matrix
@@ -233,12 +232,12 @@ protected:
     Vectord m_explicitExternalForce;   ///> Vector of explicitly defined external forces
 
     // Dirichlet boundary conditions
-    std::vector<std::size_t> m_fixedNodeIds;
+    std::vector<std::size_t> m_fixedNodeIds;    ///> Nodal IDs of the nodes that are fixed
 
-    std::shared_ptr<vega::VolumetricMesh> m_vegaPhysicsMesh;
-    std::shared_ptr<vega::SparseMatrix> m_vegaMassMatrix;///> Vega mass matrix
-    std::shared_ptr<vega::SparseMatrix> m_vegaTangentStiffnessMatrix;///> Vega Tangent stiffness matrix
-    std::shared_ptr<vega::SparseMatrix> m_vegaDampingMatrix;///> Vega Laplacian damping matrix
+    std::shared_ptr<vega::VolumetricMesh> m_vegaPhysicsMesh;            ///> Mesh used for Physics
+    std::shared_ptr<vega::SparseMatrix> m_vegaMassMatrix;               ///> Vega mass matrix
+    std::shared_ptr<vega::SparseMatrix> m_vegaTangentStiffnessMatrix;   ///> Vega Tangent stiffness matrix
+    std::shared_ptr<vega::SparseMatrix> m_vegaDampingMatrix;            ///> Vega Laplacian damping matrix
 };
 
 } // imstk

Base/SceneElements/Objects/imstkDeformableObject.cpp

@@ -24,5 +24,27 @@
 namespace imstk
 {
 
+Vectord&
+DeformableObject::getContactForce()
+{
+    auto defModel = std::dynamic_pointer_cast<imstk::DeformableBodyModel>(m_dynamicalModel);
+
+    if (!defModel)
+    {
+        LOG(WARNING) << "Dynamics pointer cast failure in DeformableObject::getContactForce()";
+    }
+
+    return defModel->getContactForce();
+}
+
+void
+DeformableObject::setDynamicalModel(std::shared_ptr<DynamicalModel> dynaDefModel)
+{
+    if (!dynaDefModel || dynaDefModel->getType() != DynamicalModel::Type::elastoDynamics)
+    {
+        LOG(WARNING) << "Dynamic model set is not of expected type (elastodynamics)!";
+    }
+    m_dynamicalModel = dynaDefModel;
+}
 
 }

Base/SceneElements/Objects/imstkDeformableObject.h

@@ -26,6 +26,7 @@
 #include <string>
 
 #include "imstkDynamicObject.h"
+#include "imstkDeformableBodyModel.h"
 
 #include "imstkProblemState.h"
 #include "imstkMath.h"
@@ -60,35 +61,45 @@ public:
     void initializeState();
     void initializeState(const Vectord& p, const Vectord& v);
 
+    ///
+    /// \brief Set/Get dynamical model
+    ///
+    void setDynamicalModel(std::shared_ptr<DynamicalModel> dynaDefModel) override;
+
+    ///
+    /// \brief Get the vector that holds the contact forces
+    ///
+    Vectord& getContactForce();
+
     // Get/Set States of the body
     const Vectord& getDisplacements() const
     {
-        m_dynamicalModel->getCurrentState()->getQ();
+        return m_dynamicalModel->getCurrentState()->getQ();
     }
 
     const Vectord& getPrevDisplacements() const
     {
-        m_dynamicalModel->getPreviousState()->getQ();
+        return m_dynamicalModel->getPreviousState()->getQ();
     }
 
     const Vectord& getVelocities() const
     {
-        m_dynamicalModel->getCurrentState()->getQDot();
+        return m_dynamicalModel->getCurrentState()->getQDot();
     }
 
     const Vectord& getPrevVelocities() const
     {
-        m_dynamicalModel->getPreviousState()->getQDot();
+        return m_dynamicalModel->getPreviousState()->getQDot();
     }
 
     const Vectord& getAccelerations() const
     {
-        m_dynamicalModel->getCurrentState()->getQDotDot();
+        return m_dynamicalModel->getCurrentState()->getQDotDot();
     }
 
     const Vectord& getPrevAccelerations() const
     {
-        m_dynamicalModel->getPreviousState()->getQDotDot();
+        return m_dynamicalModel->getPreviousState()->getQDotDot();
     }
 
 protected:

Base/SceneElements/Objects/imstkDynamicObject.cpp

@@ -60,6 +60,18 @@ DynamicObject::setPhysicsToVisualMap(std::shared_ptr<GeometryMap> map)
     m_physicsToVisualGeomMap = map;
 }
 
+std::shared_ptr<DynamicalModel>
+DynamicObject::getDynamicalModel() const
+{
+    return m_dynamicalModel;
+}
+
+void
+DynamicObject::setDynamicalModel(std::shared_ptr<DynamicalModel> dynaModel)
+{
+    m_dynamicalModel = dynaModel;
+}
+
 size_t
 DynamicObject::getNumOfDOF() const
 {

Base/SceneElements/Objects/imstkDynamicObject.h

@@ -59,6 +59,12 @@ public:
     std::shared_ptr<GeometryMap> getPhysicsToVisualMap() const;
     void setPhysicsToVisualMap(std::shared_ptr<GeometryMap> map);
 
+    ///
+    /// \brief Set/Get dynamical model
+    ///
+    std::shared_ptr<DynamicalModel> getDynamicalModel() const;
+    virtual void setDynamicalModel(std::shared_ptr<DynamicalModel> dynaModel);
+
     ///
     /// \brief Returns the number of degree of freedom
     ///

Base/SceneElements/Objects/imstkDynamicalModel.h

@@ -103,6 +103,14 @@ public:
     ///
     virtual void updateBodyStates(const Vectord& q) = 0;
 
+    ///
+    /// \brief Get the type of the object
+    ///
+    const Type& getType() const
+    {
+        return m_type;
+    }
+
 protected:
 
     Type m_type; ///> Mathematical model type

Base/Solvers/imstkNonLinearSolver.cpp

@@ -164,12 +164,12 @@ NonLinearSolver::getArmijoMax() const
     return this->m_armijoMax;
 }
 
-void NonLinearSolver::setSystem(NonLinearSystem* newSystem)
+void NonLinearSolver::setSystem(std::shared_ptr<NonLinearSystem> newSystem)
 {
     this->m_nonLinearSystem = newSystem;
 }
 
-NonLinearSystem*
+std::shared_ptr<NonLinearSystem>
 NonLinearSolver::getSystem() const
 {
     return this->m_nonLinearSystem;

Base/Solvers/imstkNonLinearSolver.h

@@ -89,7 +89,7 @@ public:
     double getAlpha() const;
 
     ///
-    /// \brief Set/Get ArmijoMax. Maximum number of steplength reductions.
+    /// \brief Set/Get ArmijoMax. Maximum number of step length reductions.
     ///
     /// \param newArmijoMax New iteration parameter.
     ///
@@ -101,8 +101,8 @@ public:
     ///
     /// \param newSystem Non-linear system replacement.
     ///
-    void setSystem(NonLinearSystem* newSystem);
-    NonLinearSystem* getSystem() const;
+    void setSystem(std::shared_ptr<NonLinearSystem> newSystem);
+    std::shared_ptr<NonLinearSystem> getSystem() const;
 
     ///
     /// \brief Set a customized iterate update function.
@@ -115,8 +115,9 @@ protected:
     std::array<double, 2> m_sigma;      ///< Safeguarding bounds for the line search
     double m_alpha;                     ///< Parameter to measure decrease
     size_t m_armijoMax;                 ///< Maximum number of step length reductions
-    NonLinearSystem *m_nonLinearSystem; ///< System of non-linear equations
-    UpdateIterateType m_updateIterate;  ///< Update iteration function
+
+    std::shared_ptr<NonLinearSystem> m_nonLinearSystem; ///< System of non-linear equations
+    UpdateIterateType m_updateIterate;                  ///< Update iteration function
 };
 
 }

Base/TimeIntegrators/imstkBackwardEuler.h

@@ -40,7 +40,13 @@ public:
     ///
     /// \brief Constructor
     ///
-    BackwardEuler(const double dT) : TimeIntegrator(dT), m_type(Type::BackwardEuler){};
+    BackwardEuler(const double dT = 0.01) : TimeIntegrator(dT)
+    {
+        m_type = Type::BackwardEuler;
+        m_alpha = { { 1, 0, 0 } };
+        m_beta = { { 1, -1, 0 } };
+        m_gamma = { { 1, -2, -1 } };
+    };
 
     ///
     /// \brief Destructor
@@ -73,9 +79,9 @@ public:
 protected:
 
     // Coefficients of the time integrator
-    std::array<double, 3> m_alpha = { { 1, 0, 0 } };
-    std::array<double, 3> m_beta = { { 1, -1, 0 } };
-    std::array<double, 3> m_gamma = { { 1, -2, -1 } };
+    std::array<double, 3> m_alpha;
+    std::array<double, 3> m_beta;
+    std::array<double, 3> m_gamma;
 
     double m_dT; ///> Delta T
 };