diff --git a/Examples/DeformableBody/DeformableBodyExample.cpp b/Examples/DeformableBody/DeformableBodyExample.cpp
index a49faf30dcbb8f0bb199bc3bf57ddbd0510369dd..76d0b99e242809282d922648e3075c49a9669e40 100644
--- a/Examples/DeformableBody/DeformableBodyExample.cpp
+++ b/Examples/DeformableBody/DeformableBodyExample.cpp
@@ -23,7 +23,7 @@
#include "imstkSimulationManager.h"
#include "imstkDeformableObject.h"
#include "imstkBackwardEuler.h"
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
#include "imstkNewtonSolver.h"
#include "imstkGaussSeidel.h"
#include "imstkPlane.h"
@@ -112,7 +112,7 @@ main()
scene->addSceneObject(planeObj);
// create a nonlinear system
- auto nlSystem = std::make_shared<NonLinearSystem>(
+ auto nlSystem = std::make_shared<NonLinearSystem<SparseMatrixd>>(
dynaModel->getFunction(),
dynaModel->getFunctionGradient());
@@ -130,7 +130,7 @@ main()
}
// create a non-linear solver and add to the scene
- auto nlSolver = std::make_shared<NewtonSolver>();
+ auto nlSolver = std::make_shared<NewtonSolver<SparseMatrixd>>();
nlSolver->setLinearSolver(linSolver);
nlSolver->setSystem(nlSystem);
scene->addNonlinearSolver(nlSolver);
diff --git a/Examples/NonlinearSolver/CMakeLists.txt b/Examples/NonlinearSolver/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a66dd3fc570afa0725ae5830d17ba193ce86e329
--- /dev/null
+++ b/Examples/NonlinearSolver/CMakeLists.txt
@@ -0,0 +1,35 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+###########################################################################
+
+project(Example-NonlinearSolver)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} nonlinearSolver.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} Solvers Common)
+
diff --git a/Examples/NonlinearSolver/nonlinearSolver.cpp b/Examples/NonlinearSolver/nonlinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1725f74c8351bfb07f27fed3fbb7c273434782d9
--- /dev/null
+++ b/Examples/NonlinearSolver/nonlinearSolver.cpp
@@ -0,0 +1,91 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include <ios>
+#include <iostream>
+#include <vector>
+#include <iomanip>
+
+#include "imstkMath.h"
+#include "imstkNewtonSolver.h"
+#include "imstkDirectLinearSolver.h"
+#include "imstkNonLinearSolver.h"
+
+using namespace imstk;
+
+int
+main(int argc, char** argv)
+{
+ const int N = 2;
+ auto x = Vectord(N);
+ auto xe = Vectord(N);
+ auto y = Vectord(N);
+ auto A = Matrixd(N, N);
+
+ x[0] = 100.0;
+ x[1] = 100.0;
+ xe[0] = 1.0;
+ xe[1] = 10.0;
+
+ auto func = [&y](const Vectord& x, const bool isSemiImplicit) -> const Vectord& {
+ // auto y = Vectord(x.size());
+ y[0] = x[0] * x[0] - 1.0;
+ y[1] = x[1] * x[1] - 100.0;
+
+ return y;
+ };
+
+ auto jac = [&A](const Vectord& x) -> const Matrixd& {
+ // auto A = Matrixd(x.size(), x.size());
+ A(0, 0) = 2 * x[0];
+ A(0, 1) = 0.0;
+ A(1, 0) = 0.0;
+ A(1, 1) = 2 * x[1];
+
+ return A;
+ };
+
+ auto updateX = [&x](const Vectord& du, const bool isSemiImplicit)
+ {
+ x -= du;
+ return;
+ };
+
+ auto updateXold = [](void) {};
+
+ auto nlSystem = std::make_shared<NonLinearSystem<Matrixd>>(func, jac);
+ nlSystem->setUnknownVector(x);
+ nlSystem->setUpdateFunction(updateX);
+ nlSystem->setUpdatePreviousStatesFunction(updateXold);
+
+ auto linSolver = std::make_shared<DirectLinearSolver<Matrixd>>();
+ auto nlSolver = std::make_shared<NewtonSolver<Matrixd>>();
+ nlSolver->setMaxIterations(100);
+ nlSolver->setRelativeTolerance(1e-8);
+ nlSolver->setAbsoluteTolerance(1e-10);
+ nlSolver->setSystem(nlSystem);
+ nlSolver->setLinearSolver(linSolver);
+
+ std::cout << "init_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl;
+ nlSolver->solve();
+
+ std::cout << "final_error = " << std::setprecision(12) << std::scientific << (x-xe).norm() << std::endl;
+}
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index fd2f00405633c7fa323f53093af627720964cf38..02d42bdfcfbfad82ad859c45fbcf1eb9547f8c2e 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -261,7 +261,7 @@ FEMDeformableBodyModel::loadBoundaryConditions()
if (file.is_open())
{
size_t index;
- auto maxAllowed = m_vegaPhysicsMesh->getNumVertices();
+ size_t maxAllowed = m_vegaPhysicsMesh->getNumVertices();
while (!file.eof())
{
file >> index;
@@ -689,7 +689,7 @@ FEMDeformableBodyModel::updateBodyIntermediateStates(
m_qSol = m_currentState->getQ();
}
-NonLinearSystem::VectorFunctionType
+NonLinearSystem<SparseMatrixd>::VectorFunctionType
FEMDeformableBodyModel::getFunction()
{
#pragma warning( push )
@@ -711,7 +711,7 @@ FEMDeformableBodyModel::getFunction()
#pragma warning( pop )
}
-NonLinearSystem::MatrixFunctionType
+NonLinearSystem<SparseMatrixd>::MatrixFunctionType
FEMDeformableBodyModel::getFunctionGradient()
{
#pragma warning( push )
@@ -731,7 +731,7 @@ FEMDeformableBodyModel::getFunctionGradient()
#pragma warning( pop )
}
-NonLinearSystem::UpdateFunctionType
+NonLinearSystem<SparseMatrixd>::UpdateFunctionType
FEMDeformableBodyModel::getUpdateFunction()
{
// Function to evaluate the nonlinear objective function given the current state
@@ -743,7 +743,7 @@ FEMDeformableBodyModel::getUpdateFunction()
};
}
-NonLinearSystem::UpdatePrevStateFunctionType
+NonLinearSystem<SparseMatrixd>::UpdatePrevStateFunctionType
FEMDeformableBodyModel::getUpdatePrevStateFunction()
{
// Function to evaluate the nonlinear objective function given the current state
diff --git a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index ce911e5ef0d97f7573ec98b683ecdbd2ec185924..9641e0cba7dfb8409c037f238f0619cd75070c82 100644
--- a/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Source/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -24,7 +24,7 @@
// imstk
#include "imstkDynamicalModel.h"
#include "imstkInternalForceModelTypes.h"
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
// vega
#include "sparseMatrix.h"
@@ -38,7 +38,9 @@ namespace imstk
{
class TimeIntegrator;
class VegaMeshIO;
-class NewtonSolver;
+// class NewtonSolver<SparseMatrixd>;
+// template <typename SystemMatrix>
+// class NewtonSolver;
class InternalForceModel;
struct FEMModelConfig
@@ -72,6 +74,8 @@ struct FEMModelConfig
class FEMDeformableBodyModel : public DynamicalModel<FeDeformBodyState>
{
using kinematicState = FeDeformBodyState;
+using System = NonLinearSystem<SparseMatrixd>;
+
public:
///
/// \brief Constructor
@@ -209,19 +213,19 @@ public:
/// \brief Returns the "function" that evaluates the nonlinear function given
/// the state vector
///
- NonLinearSystem::VectorFunctionType getFunction();
+ System::VectorFunctionType getFunction();
///
/// \brief Get the function that updates the model given the solution
///
- NonLinearSystem::UpdateFunctionType getUpdateFunction();
- NonLinearSystem::UpdatePrevStateFunctionType getUpdatePrevStateFunction();
+ System::UpdateFunctionType getUpdateFunction();
+ System::UpdatePrevStateFunctionType getUpdatePrevStateFunction();
///
/// \brief Returns the "function" that evaluates the gradient of the nonlinear
/// function given the state vector
///
- NonLinearSystem::MatrixFunctionType getFunctionGradient();
+ System::MatrixFunctionType getFunctionGradient();
///
/// \brief Get the contact force vector
@@ -268,7 +272,7 @@ public:
protected:
std::shared_ptr<InternalForceModel> m_internalForceModel; ///> Mathematical model for intenal forces
std::shared_ptr<TimeIntegrator> m_timeIntegrator; ///> Time integrator
- std::shared_ptr<NonLinearSystem> m_nonLinearSystem; ///> Nonlinear system resulting from TI and force model
+ std::shared_ptr<System> m_nonLinearSystem; ///> Nonlinear system resulting from TI and force model
std::shared_ptr<FEMModelConfig> m_FEModelConfig;
diff --git a/Source/Solvers/imstkDirectLinearSolver.cpp b/Source/Solvers/imstkDirectLinearSolver.cpp
index e8ee18238bf9700fa08a86e478a024e876223c80..77b4018ee1db170783bed698d08968edaaf57dce 100644
--- a/Source/Solvers/imstkDirectLinearSolver.cpp
+++ b/Source/Solvers/imstkDirectLinearSolver.cpp
@@ -21,6 +21,7 @@
#include "imstkDirectLinearSolver.h"
#include "imstkLinearSystem.h"
+#include "imstkLogger.h"
#include <iostream>
namespace imstk
@@ -58,17 +59,14 @@ setSystem(std::shared_ptr<LinearSystem<SparseMatrixd>> newSystem)
void
DirectLinearSolver<SparseMatrixd>::solve(const Vectord& rhs, Vectord& x)
{
+ if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
x = m_solver.solve(rhs);
}
void
DirectLinearSolver<SparseMatrixd>::solve(Vectord& x)
{
- if (!m_linearSystem)
- {
- return;
- }
-
+ if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
x.setZero();
auto b = m_linearSystem->getRHSVector();
@@ -78,21 +76,17 @@ DirectLinearSolver<SparseMatrixd>::solve(Vectord& x)
void
DirectLinearSolver<Matrixd>::solve(const Vectord& rhs, Vectord& x)
{
+ if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
x = m_solver.solve(rhs);
}
void
DirectLinearSolver<Matrixd>::solve(Vectord& x)
{
- if (!m_linearSystem)
- {
- /// \todo Log this
- return;
- }
-
+ if(!m_linearSystem) LOG(FATAL) << "Linear system has not been set";
x.setZero();
auto b = m_linearSystem->getRHSVector();
x = m_solver.solve(b);
}
-} // imstk
\ No newline at end of file
+} // imstk
diff --git a/Source/Solvers/imstkDirectLinearSolver.h b/Source/Solvers/imstkDirectLinearSolver.h
index 848fd52bf19692657ee8b89bda0982f476ffb24c..e586d811f03ec1e8fb0b21dd0cdc961d23616017 100644
--- a/Source/Solvers/imstkDirectLinearSolver.h
+++ b/Source/Solvers/imstkDirectLinearSolver.h
@@ -46,7 +46,8 @@ public:
///
/// \brief Default constructor/destructor.
///
- DirectLinearSolver() = delete;
+ // DirectLinearSolver() = delete;
+ DirectLinearSolver() {}
~DirectLinearSolver() {};
///
@@ -69,6 +70,11 @@ public:
///
void setSystem(std::shared_ptr<LinearSystemType> newSystem) override;
+ ///
+ /// \brief Sets the matrix
+ ///
+ void setSystem(std::shared_ptr<Matrixd> matrix);
+
///
/// \brief Returns true if the solver is iterative
///
diff --git a/Source/Solvers/imstkLinearSolver.h b/Source/Solvers/imstkLinearSolver.h
index e49596b25b1e3709cd5e7dbf3ebe36a753be2c37..7b279c1db9f4828b5996159d90acf2d98b0d1c16 100644
--- a/Source/Solvers/imstkLinearSolver.h
+++ b/Source/Solvers/imstkLinearSolver.h
@@ -52,7 +52,7 @@ public:
///
/// \brief Default constructor/destructor
///
- LinearSolver() {};
+ LinearSolver() : m_linearSystem(nullptr) {};
virtual ~LinearSolver() {};
///
diff --git a/Source/Solvers/imstkLinearSystem.h b/Source/Solvers/imstkLinearSystem.h
index 660b3ba0c8eacf14a32ba47d1808d941895ad5c5..5c0814258a80dd0e654bbfcdaf44e4e519d12b3b 100644
--- a/Source/Solvers/imstkLinearSystem.h
+++ b/Source/Solvers/imstkLinearSystem.h
@@ -21,7 +21,7 @@
#pragma once
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
namespace imstk
{
@@ -30,7 +30,7 @@ namespace imstk
/// \brief Represents the linear system of the form \f$ Ax = b \f$
///
template<typename SystemMatrixType>
-class LinearSystem : public NonLinearSystem
+class LinearSystem : public NonLinearSystem<SystemMatrixType>
{
public:
///
@@ -91,7 +91,7 @@ public:
///
/// \brief Set the system matrix corresponding to this ODE system.
///
- void setMatrix(const SparseMatrixd& newMatrix)
+ void setMatrix(const SystemMatrixType& newMatrix)
{
m_A = newMatrix;
}
diff --git a/Source/Solvers/imstkNewtonSolver.cpp b/Source/Solvers/imstkNewtonSolver.cpp
index 70fbb6074914cf7ca725f93f4f804d5f097c2280..7632722119b1635cff0efc2ee94d6dd36e8faf57 100644
--- a/Source/Solvers/imstkNewtonSolver.cpp
+++ b/Source/Solvers/imstkNewtonSolver.cpp
@@ -20,15 +20,18 @@
=========================================================================*/
#include <iostream>
+#include <memory>
#include "imstkNewtonSolver.h"
#include "imstkIterativeLinearSolver.h"
#include "imstkConjugateGradient.h"
+#include "imstkDirectLinearSolver.h"
#include "imstkLogger.h"
namespace imstk
{
-NewtonSolver::NewtonSolver() :
+template <>
+NewtonSolver<SparseMatrixd>::NewtonSolver() :
m_linearSolver(std::make_shared<ConjugateGradient>()),
m_forcingTerm(0.9),
m_absoluteTolerance(1e-3),
@@ -40,17 +43,32 @@ NewtonSolver::NewtonSolver() :
{
}
+template <>
+NewtonSolver<Matrixd>::NewtonSolver() :
+ m_linearSolver(std::make_shared<DirectLinearSolver<Matrixd>>()),
+ m_forcingTerm(0.9),
+ m_absoluteTolerance(1e-3),
+ m_relativeTolerance(1e-6),
+ m_gamma(0.9),
+ m_etaMax(0.9),
+ m_maxIterations(1),
+ m_useArmijo(true)
+{
+}
+
+
+template <typename SystemMatrix>
void
-NewtonSolver::solveGivenState(Vectord& x)
+NewtonSolver<SystemMatrix>::solveGivenState(Vectord& x)
{
- if (!m_nonLinearSystem)
+ 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 = m_nonLinearSystem->evaluateF(x, m_isSemiImplicit).norm();
+ double fnorm = this->m_nonLinearSystem->evaluateF(x, this->m_isSemiImplicit).norm();
double stopTolerance = m_absoluteTolerance + m_relativeTolerance * fnorm;
m_linearSolver->setTolerance(stopTolerance);
@@ -65,7 +83,7 @@ NewtonSolver::solveGivenState(Vectord& x)
}
this->updateJacobian(x);
m_linearSolver->solve(dx);
- m_updateIterate(-dx, x);
+ this->m_updateIterate(-dx, x);
double newNorm = this->armijo(dx, x, fnorm);
@@ -83,23 +101,24 @@ NewtonSolver::solveGivenState(Vectord& x)
}
}
+template <typename SystemMatrix>
void
-NewtonSolver::solve()
+NewtonSolver<SystemMatrix>::solve()
{
- if (!m_nonLinearSystem)
+ if (!this->m_nonLinearSystem)
{
LOG(WARNING) << "NewtonMethod::solve - nonlinear system is not set to the nonlinear solver";
return;
}
size_t iterNum;
- const auto& u = m_nonLinearSystem->getUnknownVector();
+ const auto& u = this->m_nonLinearSystem->getUnknownVector();
Vectord du = u; // make this a class member in future
+ double error0 = MAX_D;
double epsilon = m_relativeTolerance * m_relativeTolerance;
for (iterNum = 0; iterNum < m_maxIterations; ++iterNum)
{
- double error0 = MAX_D;
double error = updateJacobian(u);
if (iterNum == 0)
@@ -109,49 +128,52 @@ NewtonSolver::solve()
if (error / error0 < epsilon && iterNum > 0)
{
- //std::cout << "Num. of Newton Iterations: " << i << "\tError ratio: " << error/error0 << std::endl;
+ // std::cout << "Num. of Newton Iterations: " << iterNum << "\tError ratio: " << error/error0 << ", " << error << " " << error0 << std::endl;
break;
}
m_linearSolver->solve(du);
- m_nonLinearSystem->m_FUpdate(du, m_isSemiImplicit);
+ this->m_nonLinearSystem->m_FUpdate(du, this->m_isSemiImplicit);
}
- m_nonLinearSystem->m_FUpdatePrevState();
- if (iterNum == m_maxIterations && !m_isSemiImplicit)
+ this->m_nonLinearSystem->m_FUpdatePrevState();
+
+ if (iterNum == m_maxIterations && !this->m_isSemiImplicit)
{
LOG(WARNING) << "NewtonMethod::solve - The solver did not converge after max. iterations";
}
}
+template <typename SystemMatrix>
double
-NewtonSolver::updateJacobian(const Vectord& x)
+NewtonSolver<SystemMatrix>::updateJacobian(const Vectord& x)
{
// Evaluate the Jacobian and sets the matrix
- if (!m_nonLinearSystem)
+ if (!this->m_nonLinearSystem)
{
LOG(WARNING) << "NewtonMethod::updateJacobian - nonlinear system is not set to the nonlinear solver";
return -1;
}
- auto& A = m_nonLinearSystem->m_dF(x);
+ auto& A = this->m_nonLinearSystem->m_dF(x);
if (A.innerSize() == 0)
{
LOG(WARNING) << "NewtonMethod::updateJacobian - Size of matrix is 0!";
return -1;
}
- auto& b = m_nonLinearSystem->m_F(x, m_isSemiImplicit);
+ auto& b = this->m_nonLinearSystem->m_F(x, this->m_isSemiImplicit);
- auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(A, b);
- //linearSystem->setLinearProjectors(m_nonLinearSystem->getLinearProjectors()); /// \todo Left for near future reference. Clear in future.
+ auto linearSystem = std::make_shared<typename LinearSolverType::LinearSystemType>(A, b);
+ //linearSystem->setLinearProjectors(this->m_nonLinearSystem->getLinearProjectors()); /// \todo Left for near future reference. Clear in future.
m_linearSolver->setSystem(linearSystem);
- return b.dot(b);
+ return std::sqrt(b.dot(b));
}
+template <typename SystemMatrix>
void
-NewtonSolver::updateForcingTerm(const double ratio, const double stopTolerance, const double fnorm)
+NewtonSolver<SystemMatrix>::updateForcingTerm(const double ratio, const double stopTolerance, const double fnorm)
{
double eta = m_gamma * ratio * ratio;
double forcingTermSqr = m_forcingTerm * m_forcingTerm;
@@ -166,27 +188,36 @@ NewtonSolver::updateForcingTerm(const double ratio, const double stopTolerance,
m_forcingTerm = std::max(std::min(eta, m_etaMax), 0.5 * stopTolerance / fnorm);
}
+template <typename SystemMatrix>
void
-NewtonSolver::setLinearSolver(std::shared_ptr<NewtonSolver::LinearSolverType> newLinearSolver)
+NewtonSolver<SystemMatrix>::setLinearSolver(std::shared_ptr<NewtonSolver::LinearSolverType> newLinearSolver)
{
m_linearSolver = newLinearSolver;
}
-std::shared_ptr<NewtonSolver::LinearSolverType>
-NewtonSolver::getLinearSolver() const
+template <typename SystemMatrix>
+// std::shared_ptr<NewtonSolver<SystemMatrix>::LinearSolverType>
+auto
+NewtonSolver<SystemMatrix>::getLinearSolver() const -> std::shared_ptr<LinearSolverType>
{
return m_linearSolver;
}
+template <typename SystemMatrix>
void
-NewtonSolver::setAbsoluteTolerance(const double aTolerance)
+NewtonSolver<SystemMatrix>::setAbsoluteTolerance(const double aTolerance)
{
m_absoluteTolerance = aTolerance;
}
+template <typename SystemMatrix>
double
-NewtonSolver::getAbsoluteTolerance() const
+NewtonSolver<SystemMatrix>::getAbsoluteTolerance() const
{
return m_absoluteTolerance;
}
+
+
+template class NewtonSolver<SparseMatrixd>;
+template class NewtonSolver<Matrixd>;
} // imstk
diff --git a/Source/Solvers/imstkNewtonSolver.h b/Source/Solvers/imstkNewtonSolver.h
index 047f58392e255641d9034d534b8a4d69fc9c677d..b06e0d977f2441602f3347abe80838107a1d1f30 100644
--- a/Source/Solvers/imstkNewtonSolver.h
+++ b/Source/Solvers/imstkNewtonSolver.h
@@ -33,10 +33,11 @@ namespace imstk
/// root of F or fails. Global convergence is achieved using a line search sub-process
/// and the Armijo rule
///
-class NewtonSolver : public NonLinearSolver
+template <typename SystemMatrix>
+class NewtonSolver : public NonLinearSolver<SystemMatrix>
{
public:
- using LinearSolverType = LinearSolver<SparseMatrixd>;
+ using LinearSolverType = LinearSolver<SystemMatrix>;
public:
///
@@ -78,7 +79,8 @@ public:
///
/// \brief Get LinearSolver
///
- std::shared_ptr<LinearSolverType> getLinearSolver() const;
+ // std::shared_ptr<LinearSolverType> getLinearSolver() const;
+ auto getLinearSolver() const -> std::shared_ptr<LinearSolverType>;
///
/// \brief Update jacobians
@@ -90,7 +92,7 @@ public:
///
/// \brief Get JacobianMatrix. Returns jacobian matrix
///
- SparseMatrixd& getJacobianMatrix();
+ SystemMatrix& getJacobianMatrix();
///
/// \brief Set AbsoluteTolerance
@@ -169,7 +171,7 @@ public:
void setUseArmijo(const bool value)
{
m_useArmijo = value;
- (value) ? m_armijoMax = 30 : m_armijoMax = 0;
+ (value) ? this->m_armijoMax = 30 : this->m_armijoMax = 0;
}
///
@@ -193,7 +195,7 @@ public:
///
virtual void setToFullyImplicit() override
{
- m_isSemiImplicit = false;
+ this->m_isSemiImplicit = false;
}
///
@@ -201,7 +203,7 @@ public:
///
virtual void setToSemiImplicit() override
{
- m_isSemiImplicit = true;
+ this->m_isSemiImplicit = true;
m_maxIterations = 1;
}
diff --git a/Source/Solvers/imstkNonLinearSolver.cpp b/Source/Solvers/imstkNonLinearSolver.cpp
index efad583aba8c50c47e9cfbcf16f19054e04e2c60..c0f7f8b80454b93ae0af5609d1b68650c8b8b073 100644
--- a/Source/Solvers/imstkNonLinearSolver.cpp
+++ b/Source/Solvers/imstkNonLinearSolver.cpp
@@ -20,10 +20,12 @@
=========================================================================*/
#include "imstkNonLinearSolver.h"
+#include "imstkMath.h"
namespace imstk
{
-NonLinearSolver::NonLinearSolver() : m_sigma(std::array<double, 2>
+template <typename SystemMatrix>
+NonLinearSolver<SystemMatrix>::NonLinearSolver() : m_sigma(std::array<double, 2>
{
{ 0.1, 0.5 }
}),
@@ -36,8 +38,9 @@ NonLinearSolver::NonLinearSolver() : m_sigma(std::array<double, 2>
};
}
+template <typename SystemMatrix>
double
-NonLinearSolver::armijo(const Vectord& dx, Vectord& x, const double previousFnorm)
+NonLinearSolver<SystemMatrix>::armijo(const Vectord& dx, Vectord& x, const double previousFnorm)
{
/// Temporaries used in the line search
std::array<double, 3> fnormSqr = { previousFnorm* previousFnorm, 0.0, 0.0 };
@@ -96,8 +99,9 @@ NonLinearSolver::armijo(const Vectord& dx, Vectord& x, const double previousFnor
return currentFnorm;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::parabolicModel(const std::array<double, 3>& fnorm, std::array<double, 3>& lambda)
+NonLinearSolver<SystemMatrix>::parabolicModel(const std::array<double, 3>& fnorm, std::array<double, 3>& lambda)
{
/// Compute the coefficients for the interpolation polynomial:
/// p(lambda) = fnorm[0] + (b*lambda + a*lambda^2)/d1, where
@@ -130,57 +134,69 @@ NonLinearSolver::parabolicModel(const std::array<double, 3>& fnorm, std::array<d
lambda[0] = newLambda;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::setSigma(const std::array<double, 2>& newSigma)
+NonLinearSolver<SystemMatrix>::setSigma(const std::array<double, 2>& newSigma)
{
m_sigma = newSigma;
}
+template <typename SystemMatrix>
const std::array<double, 2>&
-NonLinearSolver::getSigma() const
+NonLinearSolver<SystemMatrix>::getSigma() const
{
return m_sigma;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::setAlpha(const double newAlpha)
+NonLinearSolver<SystemMatrix>::setAlpha(const double newAlpha)
{
m_alpha = newAlpha;
}
+template <typename SystemMatrix>
double
-NonLinearSolver::getAlpha() const
+NonLinearSolver<SystemMatrix>::getAlpha() const
{
return m_alpha;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::setArmijoMax(const size_t newArmijoMax)
+NonLinearSolver<SystemMatrix>::setArmijoMax(const size_t newArmijoMax)
{
m_armijoMax = newArmijoMax;
}
+template <typename SystemMatrix>
size_t
-NonLinearSolver::getArmijoMax() const
+NonLinearSolver<SystemMatrix>::getArmijoMax() const
{
return m_armijoMax;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::setSystem(std::shared_ptr<NonLinearSystem> newSystem)
+NonLinearSolver<SystemMatrix>::setSystem(std::shared_ptr<NonLinearSystem<SystemMatrix>> newSystem)
{
m_nonLinearSystem = newSystem;
}
-std::shared_ptr<NonLinearSystem>
-NonLinearSolver::getSystem() const
+template <typename SystemMatrix>
+std::shared_ptr<NonLinearSystem<SystemMatrix>>
+NonLinearSolver<SystemMatrix>::getSystem() const
{
return m_nonLinearSystem;
}
+template <typename SystemMatrix>
void
-NonLinearSolver::setUpdateIterate(const NonLinearSolver::UpdateIterateType& newUpdateIterate)
+NonLinearSolver<SystemMatrix>::setUpdateIterate(const NonLinearSolver<SystemMatrix>::UpdateIterateType& newUpdateIterate)
{
m_updateIterate = newUpdateIterate;
}
-} // imstk
\ No newline at end of file
+
+template class NonLinearSolver<SparseMatrixd>;
+template class NonLinearSolver<Matrixd>;
+} // imstk
diff --git a/Source/Solvers/imstkNonLinearSolver.h b/Source/Solvers/imstkNonLinearSolver.h
index 05c2ea9604bb521074b48c2b2020c6fd10e254f0..3902fb310a54ca81f1e11733ed556ce8efcc15de 100644
--- a/Source/Solvers/imstkNonLinearSolver.h
+++ b/Source/Solvers/imstkNonLinearSolver.h
@@ -26,7 +26,7 @@
#include <functional>
// imstk includes
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
#include "imstkMath.h"
#include "imstkSolverBase.h"
@@ -35,12 +35,13 @@ namespace imstk
///
/// \brief Base class for non-linear solvers
///
+template <typename SystemMatrix>
class NonLinearSolver : public SolverBase
{
public:
- using JacobianType = std::function<const SparseMatrixd& (const Vectord&)>;
+ using JacobianType = std::function<const SystemMatrix& (const Vectord&)>;
using UpdateIterateType = std::function<void (const Vectord&, Vectord&)>;
- using FunctionType = NonLinearSystem::VectorFunctionType;
+ using FunctionType = NonLinearSystem<SparseMatrixd>::VectorFunctionType;
public:
///
@@ -105,8 +106,8 @@ public:
///
/// \param newSystem Non-linear system replacement.
///
- void setSystem(std::shared_ptr<NonLinearSystem> newSystem);
- std::shared_ptr<NonLinearSystem> getSystem() const;
+ void setSystem(std::shared_ptr<NonLinearSystem<SystemMatrix>> newSystem);
+ std::shared_ptr<NonLinearSystem<SystemMatrix>> getSystem() const;
///
/// \brief Set a customized iterate update function.
@@ -136,7 +137,7 @@ protected:
double m_alpha; ///< Parameter to measure decrease
size_t m_armijoMax; ///< Maximum number of step length reductions
- std::shared_ptr<NonLinearSystem> m_nonLinearSystem; ///< System of non-linear equations
+ std::shared_ptr<NonLinearSystem<SystemMatrix>> m_nonLinearSystem; ///< System of non-linear equations
UpdateIterateType m_updateIterate; ///< Update iteration function
bool m_isSemiImplicit = true; ///> Semi-Implicit solver
};
diff --git a/Source/Solvers/imstkNonlinearSystem.cpp b/Source/Solvers/imstkNonLinearSystem.cpp
similarity index 58%
rename from Source/Solvers/imstkNonlinearSystem.cpp
rename to Source/Solvers/imstkNonLinearSystem.cpp
index 8da916dd00e427f72673dbe9523a2e453e114d26..92c688c54b8a1ac4b4e0409814c7fdfb58b6bc80 100644
--- a/Source/Solvers/imstkNonlinearSystem.cpp
+++ b/Source/Solvers/imstkNonLinearSystem.cpp
@@ -19,37 +19,47 @@
=========================================================================*/
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
+#include "imstkMath.h"
-#include <g3log/g3log.hpp>
+// #include <g3log/g3log.hpp>
namespace imstk
{
-NonLinearSystem::NonLinearSystem(const VectorFunctionType& F, const MatrixFunctionType& dF) : m_F(F), m_dF(dF)
+template <typename Matrix>
+NonLinearSystem<Matrix>::NonLinearSystem(const VectorFunctionType& F, const MatrixFunctionType& dF) : m_F(F), m_dF(dF)
{
}
+template <typename Matrix>
void
-NonLinearSystem::setFunction(const VectorFunctionType& function)
+NonLinearSystem<Matrix>::setFunction(const VectorFunctionType& function)
{
m_F = function;
}
+template <typename Matrix>
void
-NonLinearSystem::setJacobian(const MatrixFunctionType& function)
+NonLinearSystem<Matrix>::setJacobian(const MatrixFunctionType& function)
{
m_dF = function;
}
+template <typename Matrix>
const Vectord&
-NonLinearSystem::evaluateF(const Vectord& x, const bool isSemiImplicit)
+NonLinearSystem<Matrix>::evaluateF(const Vectord& x, const bool isSemiImplicit)
{
return m_F(x, isSemiImplicit);
}
-const SparseMatrixd&
-NonLinearSystem::evaluateJacobian(const Vectord& x)
+template <typename Matrix>
+const Matrix&
+NonLinearSystem<Matrix>::evaluateJacobian(const Vectord& x)
{
return m_dF(x);
}
-} //imstk
\ No newline at end of file
+
+template class NonLinearSystem<SparseMatrixd>;
+template class NonLinearSystem<Matrixd>;
+
+} //imstk
diff --git a/Source/Solvers/imstkNonlinearSystem.h b/Source/Solvers/imstkNonLinearSystem.h
similarity index 95%
rename from Source/Solvers/imstkNonlinearSystem.h
rename to Source/Solvers/imstkNonLinearSystem.h
index 941387064002f836d9160338239a061a439a4019..d631063b9b287bb562a671059b09ea35e92ca1bf 100644
--- a/Source/Solvers/imstkNonlinearSystem.h
+++ b/Source/Solvers/imstkNonLinearSystem.h
@@ -31,11 +31,12 @@ namespace imstk
///
/// \brief Base class for a multi-variable nonlinear system
///
+template <typename Matrix>
class NonLinearSystem
{
public:
using VectorFunctionType = std::function<const Vectord& (const Vectord&, const bool)>;
- using MatrixFunctionType = std::function<const SparseMatrixd& (const Vectord&)>;
+ using MatrixFunctionType = std::function<const Matrix& (const Vectord&)>;
using UpdateFunctionType = std::function<void (const Vectord&, const bool)>;
using UpdatePrevStateFunctionType = std::function<void ()>;
@@ -66,7 +67,7 @@ public:
///
/// \brief Evaluate gradient of the function at a given state
///
- virtual const SparseMatrixd& evaluateJacobian(const Vectord& x);
+ virtual const Matrix& evaluateJacobian(const Vectord& x);
///
/// \brief
diff --git a/Source/Solvers/imstkSOR.h b/Source/Solvers/imstkSOR.h
index f310b48bc70d1dde34f439fdf74cdf76864ade3d..d71a93fd9af8aeea16b06a432ab856f6723d3c0e 100644
--- a/Source/Solvers/imstkSOR.h
+++ b/Source/Solvers/imstkSOR.h
@@ -22,7 +22,7 @@
#pragma once
// iMSTK includes
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
#include "imstkIterativeLinearSolver.h"
// Eigen includes
diff --git a/Source/apiUtilities/imstkAPIUtilities.cpp b/Source/apiUtilities/imstkAPIUtilities.cpp
index ef6930160251d2533a82d66a71c6ee7714b018a8..82bfcb51e895fd2093cfd7cf062fa32cf12ff387 100644
--- a/Source/apiUtilities/imstkAPIUtilities.cpp
+++ b/Source/apiUtilities/imstkAPIUtilities.cpp
@@ -30,7 +30,7 @@
#include "imstkColor.h"
// Solvers
-#include "imstkNonlinearSystem.h"
+#include "imstkNonLinearSystem.h"
// Geometry
#include "imstkPlane.h"
@@ -177,12 +177,12 @@ createAndAddVisualSceneObject(std::shared_ptr<Scene> scene,
///
/// \brief Create a non-linear system using FEM dynamic model
///
-std::shared_ptr<NonLinearSystem>
+std::shared_ptr<NonLinearSystem<SparseMatrixd>>
createNonLinearSystem(std::shared_ptr<FEMDeformableBodyModel> dynaModel)
{
CHECK(dynaModel != nullptr) << "createNonLinearSystem: Dynamic model is not valid!";
- auto nlSystem = std::make_shared<NonLinearSystem>(
+ auto nlSystem = std::make_shared<NonLinearSystem<SparseMatrixd>>(
dynaModel->getFunction(),
dynaModel->getFunctionGradient());
@@ -255,4 +255,4 @@ getMeshGraph(std::shared_ptr<TetrahedralMesh> m)
return gMesh;
}
} //apiutils
-} // imstk
\ No newline at end of file
+} // imstk
diff --git a/Source/apiUtilities/imstkAPIUtilities.h b/Source/apiUtilities/imstkAPIUtilities.h
index a13b953d01e3555e384094305c87e028d98665ef..f39515b79aa620b9c7a23580a178117a2fe608cc 100644
--- a/Source/apiUtilities/imstkAPIUtilities.h
+++ b/Source/apiUtilities/imstkAPIUtilities.h
@@ -33,6 +33,7 @@ class PointSet;
class SceneManager;
class Scene;
class FEMDeformableBodyModel;
+template <typename Matrix>
class NonLinearSystem;
class CollidingObject;
class Graph;
@@ -67,7 +68,7 @@ std::shared_ptr<SceneObject> createAndAddVisualSceneObject(std::shared_ptr<Scene
///
/// \brief Create a non-linear system using FEM dynamic model
///
-std::shared_ptr<NonLinearSystem> createNonLinearSystem(std::shared_ptr<FEMDeformableBodyModel> dynaModel);
+std::shared_ptr<NonLinearSystem<SparseMatrixd>> createNonLinearSystem(std::shared_ptr<FEMDeformableBodyModel> dynaModel);
///
/// \brief Print number of updates for second for a given scene