BUG: Adds unknown vector to NL system

Adds unknown vector to nonlinear system.
Adds simplified nonlinear solver for testing

Base/SceneElements/Objects/imstkDeformableBodyModel.cpp

@@ -112,7 +112,9 @@ DeformableBodyModel::initialize(std::shared_ptr<Geometry> geometry)
 
     m_Feff.resize(m_numDOF);
     m_Finternal.resize(m_numDOF);
+    m_Finternal.setConstant(0.0);//?
     m_Fcontact.resize(m_numDOF);
+    m_qSol.resize(m_numDOF);
 }
 
 void

Base/SceneElements/Objects/imstkDeformableBodyModel.h

@@ -195,9 +195,11 @@ public:
     ///
     /// \brief Get the contact force vector
     ///
-    Vectord& getContactForce()
+    Vectord& getContactForce();
+
+    Vectord& getUnknownVec()
     {
-        return m_Fcontact;
+        return m_qSol;
     }
 
 protected:
@@ -220,6 +222,7 @@ protected:
     Vectord m_Finternal;    ///> Vector of gravity forces
     Vectord m_Feff;         ///> Vector of effective forces
     Vectord m_Fcontact;     ///> Vector of contact forces
+    Vectord m_qSol;         ///> Vector to maintain solution at each iteration of nonlinear solver
 
     // External field forces
     Vectord m_gravityForce;   ///> Vector of gravity forces

Base/Solvers/imstkNewtonMethod.cpp

@@ -43,7 +43,7 @@ NewtonMethod::solve(Vectord& x)
 {
     if(!this->m_nonLinearSystem)
     {
-        // TODO: log this
+        LOG(WARNING) << "NewtonMethod::solve - nonlinear system is not set to the nonlinear solver";
         return;
     }
 
@@ -84,6 +84,37 @@ NewtonMethod::solve(Vectord& x)
     }
 }
 
+void
+NewtonMethod::solveSimple()
+{
+    auto x = this->m_nonLinearSystem->getUnknownVector();
+
+    if (!this->m_nonLinearSystem)
+    {
+        LOG(WARNING) << "NewtonMethod::solve - nonlinear system is not set to the nonlinear solver";
+        return;
+    }
+
+    // Compute norms, set tolerances and other temporaries
+    //double fnorm = this->m_nonLinearSystem->evaluateF(x).norm();
+    double stopTolerance = 1.0e-6;
+
+    this->linearSolver->setTolerance(stopTolerance);
+
+    //Vectord dx = x;
+
+    for (size_t i = 0; i < this->maxIterations; ++i)
+    {
+        /*if (fnorm < stopTolerance)
+        {
+        return;
+        }*/
+        this->updateJacobian(x);
+        this->linearSolver->solve(x);
+        //this->m_updateIterate(-dx, x);
+    }
+}
+
 
 void
 NewtonMethod::updateForcingTerm(const double ratio, const double stopTolerance, const double fnorm)
@@ -121,19 +152,20 @@ NewtonMethod::updateJacobian(const Vectord& x)
 {
     // Evaluate the Jacobian and sets the matrix
     if(!this->m_nonLinearSystem)
-    {// TODO: Print message or log error.
+    {
+        LOG(WARNING) << "NewtonMethod::updateJacobian - nonlinear system is not set to the nonlinear solver";
         return;
     }
 
+    auto &b = this->m_nonLinearSystem->m_F(x);
     auto &A = this->m_nonLinearSystem->m_dF(x);
 
     if(A.innerSize() == 0)
-    {// TODO: Print message and/or log error.
+    {
+        // TODO: Print message and/or log error.
         return;
     }
 
-    auto &b = this->m_nonLinearSystem->m_F(x);
-
     auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(A, b);
     this->linearSolver->setSystem(linearSystem);
 }

Base/Solvers/imstkNewtonMethod.h

@@ -53,6 +53,7 @@ public:
     /// \param x Current iterate.
     ///
     void solve(Vectord &x) override;
+    void solveSimple() override;
 
     ///
     /// \brief Update forcing term according to Eisenstat-Walker criteria.

Base/Solvers/imstkNonLinearSolver.h

@@ -50,6 +50,7 @@ public:
     /// \brief Main solve routine.
     ///
     virtual void solve(Vectord& x) = 0;
+    virtual void solveSimple() = 0;
 
     ///
     /// \brief Backtracking line search method based on the Armijo-Goldstein condition

Base/Solvers/imstkNonlinearSystem.h

@@ -24,7 +24,8 @@
 
 #include "imstkMath.h"
 
-namespace imstk {
+namespace imstk
+{
 
 ///
 /// \class NonLinearSystem
@@ -35,7 +36,7 @@ class NonLinearSystem
 {
 public:
     using VectorFunctionType = std::function < const Vectord&(const Vectord&) >;
-    using MatrixFunctionType = std::function < const SparseMatrixd&(const Vectord&) > ;
+    using MatrixFunctionType = std::function < const SparseMatrixd&(const Vectord&) >;
 
 public:
     ///
@@ -66,9 +67,26 @@ public:
     ///
     virtual const SparseMatrixd& evaluateJacobian(const Vectord& x);
 
+    ///
+    /// \brief
+    ///
+    void setUnknownVector(Vectord& v)
+    {
+        m_unknown = v;
+    }
+
+    ///
+    /// \brief
+    ///
+    Vectord& getUnknownVector()
+    {
+        return m_unknown;
+    }
+
 public:
     VectorFunctionType m_F;  ///> Nonlinear function
-    MatrixFunctionType m_dF;  ///> Gradient of the Nonlinear function with respect to the unknown vector
+    MatrixFunctionType m_dF; ///> Gradient of the Nonlinear function with respect to the unknown vector
+    Vectord m_unknown;
 };
 
 }