From c3aec22dc46c1e869f16936b044a6805a412b33b Mon Sep 17 00:00:00 2001
From: Sreekanth Arikatla <sreekanth.arikatla@kitware.com>
Date: Tue, 5 Jul 2016 11:20:24 -0400
Subject: [PATCH] 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                |  1 +
 .../imstkCorotationalFEMForceModel.h          |  7 +++--
 Base/ForceModel/imstkForceModelConfig.cpp     |  3 +-
 .../imstkIsotropicHyperelasticFEMForceModel.h | 14 +++++----
 Base/ForceModel/imstkLinearFEMForceModel.h    | 16 ++++++++--
 Base/ForceModel/imstkMassSpringForceModel.h   |  1 +
 Base/ForceModel/imstkStVKForceModel.h         | 11 ++++++-
 Base/Geometry/Reader/imstkVegaMeshReader.cpp  |  1 +
 .../Objects/imstkDeformableBodyModel.cpp      | 30 ++++++++++++-------
 .../Objects/imstkDeformableBodyModel.h        | 23 +++++++-------
 .../Objects/imstkDeformableObject.cpp         | 22 ++++++++++++++
 .../Objects/imstkDeformableObject.h           | 23 ++++++++++----
 .../Objects/imstkDynamicObject.cpp            | 12 ++++++++
 .../Objects/imstkDynamicObject.h              |  6 ++++
 .../Objects/imstkDynamicalModel.h             |  8 +++++
 Base/Solvers/imstkNonLinearSolver.cpp         |  4 +--
 Base/Solvers/imstkNonLinearSolver.h           | 11 +++----
 Base/TimeIntegrators/imstkBackwardEuler.h     | 14 ++++++---
 18 files changed, 156 insertions(+), 51 deletions(-)

diff --git a/Base/ForceModel/CMakeLists.txt b/Base/ForceModel/CMakeLists.txt
index 31fa6b2c8..920e9497c 100644
--- a/Base/ForceModel/CMakeLists.txt
+++ b/Base/ForceModel/CMakeLists.txt
@@ -11,6 +11,7 @@ imstk_add_library( ForceModel
 	VegaFEM::forceModel
     VegaFEM::stvk
     VegaFEM::graph
+	VegaFEM::volumetricMesh
   )
 
 #-----------------------------------------------------------------------------
diff --git a/Base/ForceModel/imstkCorotationalFEMForceModel.h b/Base/ForceModel/imstkCorotationalFEMForceModel.h
index 4f69cb50b..636103ce8 100644
--- a/Base/ForceModel/imstkCorotationalFEMForceModel.h
+++ b/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();
 
diff --git a/Base/ForceModel/imstkForceModelConfig.cpp b/Base/ForceModel/imstkForceModelConfig.cpp
index d95544785..d55afbf8c 100644
--- a/Base/ForceModel/imstkForceModelConfig.cpp
+++ b/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;
     }
 }
diff --git a/Base/ForceModel/imstkIsotropicHyperelasticFEMForceModel.h b/Base/ForceModel/imstkIsotropicHyperelasticFEMForceModel.h
index ea977bf78..52bdbb967 100644
--- a/Base/ForceModel/imstkIsotropicHyperelasticFEMForceModel.h
+++ b/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)
diff --git a/Base/ForceModel/imstkLinearFEMForceModel.h b/Base/ForceModel/imstkLinearFEMForceModel.h
index 6a1128439..b86db915d 100644
--- a/Base/ForceModel/imstkLinearFEMForceModel.h
+++ b/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;
diff --git a/Base/ForceModel/imstkMassSpringForceModel.h b/Base/ForceModel/imstkMassSpringForceModel.h
index bb3e778fb..fce6d7da2 100644
--- a/Base/ForceModel/imstkMassSpringForceModel.h
+++ b/Base/ForceModel/imstkMassSpringForceModel.h
@@ -39,6 +39,7 @@ public:
     {
         m_massSpringSystem = massSpringSystem;
     }
+    MassSpringForceModel() = delete;
     virtual ~MassSpringForceModel();
 
     void getInternalForce(const Vectord& u, Vectord& internalForce)
diff --git a/Base/ForceModel/imstkStVKForceModel.h b/Base/ForceModel/imstkStVKForceModel.h
index fc563ce2e..b8c687184 100644
--- a/Base/ForceModel/imstkStVKForceModel.h
+++ b/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;
diff --git a/Base/Geometry/Reader/imstkVegaMeshReader.cpp b/Base/Geometry/Reader/imstkVegaMeshReader.cpp
index cc26b8ffa..02986581a 100644
--- a/Base/Geometry/Reader/imstkVegaMeshReader.cpp
+++ b/Base/Geometry/Reader/imstkVegaMeshReader.cpp
@@ -100,4 +100,5 @@ VegaMeshReader::getVolumeMeshFromVegaVolumeMesh(std::shared_ptr<vega::Volumetric
         return nullptr;
     }
 }
+
 }
diff --git a/Base/SceneElements/Objects/imstkDeformableBodyModel.cpp b/Base/SceneElements/Objects/imstkDeformableBodyModel.cpp
index d4f38afdb..42b40fe86 100644
--- a/Base/SceneElements/Objects/imstkDeformableBodyModel.cpp
+++ b/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;
     };
 }
diff --git a/Base/SceneElements/Objects/imstkDeformableBodyModel.h b/Base/SceneElements/Objects/imstkDeformableBodyModel.h
index 9a2ef0d39..738f165a9 100644
--- a/Base/SceneElements/Objects/imstkDeformableBodyModel.h
+++ b/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
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.cpp b/Base/SceneElements/Objects/imstkDeformableObject.cpp
index 180eb6d05..9325d8fc5 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.cpp
+++ b/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;
+}
 
 }
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.h b/Base/SceneElements/Objects/imstkDeformableObject.h
index f127594c7..0743390f4 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.h
+++ b/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:
diff --git a/Base/SceneElements/Objects/imstkDynamicObject.cpp b/Base/SceneElements/Objects/imstkDynamicObject.cpp
index c9a8d1e7a..1c041d83f 100644
--- a/Base/SceneElements/Objects/imstkDynamicObject.cpp
+++ b/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
 {
diff --git a/Base/SceneElements/Objects/imstkDynamicObject.h b/Base/SceneElements/Objects/imstkDynamicObject.h
index a207ac275..ba27ee704 100644
--- a/Base/SceneElements/Objects/imstkDynamicObject.h
+++ b/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
     ///
diff --git a/Base/SceneElements/Objects/imstkDynamicalModel.h b/Base/SceneElements/Objects/imstkDynamicalModel.h
index f4cd42fb4..d1f370088 100644
--- a/Base/SceneElements/Objects/imstkDynamicalModel.h
+++ b/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
diff --git a/Base/Solvers/imstkNonLinearSolver.cpp b/Base/Solvers/imstkNonLinearSolver.cpp
index 1eb8dc4d5..7296445cc 100644
--- a/Base/Solvers/imstkNonLinearSolver.cpp
+++ b/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;
diff --git a/Base/Solvers/imstkNonLinearSolver.h b/Base/Solvers/imstkNonLinearSolver.h
index a9bded2c5..90fe234d9 100644
--- a/Base/Solvers/imstkNonLinearSolver.h
+++ b/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
 };
 
 }
diff --git a/Base/TimeIntegrators/imstkBackwardEuler.h b/Base/TimeIntegrators/imstkBackwardEuler.h
index 6a9df6d9a..e1abeaf92 100644
--- a/Base/TimeIntegrators/imstkBackwardEuler.h
+++ b/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
 };
-- 
GitLab