Commit 84f3390c authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Create the TimeIntegrator module and VegeFEM interfaces.

parent 6d08cada
Pipeline #4780 passed with stage
......@@ -27,11 +27,13 @@
#include <Core/RenderDelegate.h>
SceneObjectDeformable::SceneObjectDeformable() :
OdeSystem(),
renderSecondaryMesh(false),
topologyAltered(false),
pulledVertex(-1),
timestepCounter(0),
subTimestepCounter(0)
subTimestepCounter(0),
integrationScheme(TimeIntegrator::ImplicitEuler)
{
this->setRenderDelegate(
Factory<RenderDelegate>::createConcreteClassForGroup(
......
......@@ -31,10 +31,13 @@
#include "Mesh/SurfaceMesh.h"
#include "Core/SceneObject.h"
#include "Core/ErrorLog.h"
#include "VegaObjectConfig.h"
#include "Simulators/VegaObjectConfig.h"
#include "TimeIntegrators/TimeIntegrator.h"
#include "TimeIntegrators/BackwarEuler.h"
#include "TimeIntegrators/ForwardEuler.h"
/// \brief Base class for any scene object that is defmormable
class SceneObjectDeformable : public SceneObject
class SceneObjectDeformable : public SceneObject, public OdeSystem
{
public:
......@@ -138,6 +141,49 @@ public:
return this->f_ext;
}
void setTimeIntegrator(TimeIntegrator::IntegratorType integrator)
{
this->integrationScheme = integrator;
}
bool init()
{
auto thisPointer = this->safeDownCast<SceneObjectDeformable>();
switch(integrationScheme)
{
case TimeIntegrator::ImplicitEuler:
{
this->odeSolver = std::make_shared<BackwardEuler>(thisPointer);
}
case TimeIntegrator::ExplicitEuler:
{
this->odeSolver = std::make_shared<ForwardEuler>(thisPointer);
}
default:
{
std::cerr << "Invalid time integration scheme." << std::endl;
}
return false;
}
}
void update(double dt)
{
if(!this->odeSolver)
{
std::cerr << "Ode solver needs to be set." << std::endl;
return;
}
this->odeSolver->solve(currentState,newState,dt);
currentState.swap(previousState);
currentState.swap(newState);
// TODO: Check state validity
}
protected:
friend class SceneObjectDeformableRenderDelegate;
......@@ -171,6 +217,22 @@ protected:
std::shared_ptr<SurfaceMesh> primarySurfaceMesh;
std::shared_ptr<SurfaceMesh> secondarySurfaceMesh;
std::shared_ptr<TimeIntegrator> odeSolver; ///> Integration scheme
// Consecutive system states
std::shared_ptr<OdeSystemState> currentState;
std::shared_ptr<OdeSystemState> previousState;
std::shared_ptr<OdeSystemState> newState;
std::shared_ptr<OdeSystemState> initialState;
core::Matrixd M; ///> Mass matrix
core::Matrixd C; ///> Damping matrix
core::Matrixd K; ///> Stiffness matrix
core::Vectord f; ///> Force loads
TimeIntegrator::IntegratorType integrationScheme;
};
#endif // SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
......@@ -127,12 +127,17 @@ void VegaFemSceneObject::initialize()
loadFixedBC();
// make room for deformation and force vectors
// TODO: Clean these variables
u.resize(3 * numNodes);
uvel.resize(3 * numNodes);
uaccel.resize(3 * numNodes);
f_extBase.resize(3 * numNodes);
f_ext.resize(3 * numNodes);
// New state variables
this->initialState->resize(3*numNodes);
*this->currentState = *this->initialState;
loadInitialStates();
loadScriptedExternalForces();
createForceModel();
......
......@@ -44,10 +44,8 @@ void InexactNewton::solve(core::Vectord &x)
return;
}
// Set tolerances and other temporaries
this->f.resize(x.size());
this->nonLinearSystem->eval(x, this->f);
double fnorm = this->f.norm();
// Compute norms, set tolerances and other temporaries
double fnorm = this->nonLinearSystem->eval(x).norm();
double stopTolerance = this->absoluteTolerance + this->relativeTolerance * fnorm;
this->linearSolver->setTolerance(stopTolerance);
core::Vectord dx = x;
......@@ -61,16 +59,17 @@ void InexactNewton::solve(core::Vectord &x)
this->updateJacobian(x);
this->linearSolver->solve(dx);
x -= dx;
this->updateIterate(-dx,x);
double newNorm = this->armijo(dx, x);
if(this->forcingTerm > 0 && fnorm > 0)
if(this->forcingTerm > 0 && newNorm > stopTolerance)
{
double ratio = newNorm / fnorm; // Ratio of succesive residual norms
this->updateForcingTerm(ratio, stopTolerance, fnorm);
// Reset tolerance in the linear solver according to the new forcing term
// to avoid over solving of the system.
if(this->linearSolver)
{
this->linearSolver->setTolerance(this->forcingTerm);
......@@ -114,22 +113,22 @@ std::shared_ptr< InexactNewton::LinearSolverType > InexactNewton::getLinearSolve
void InexactNewton::updateJacobian(const core::Vectord &x)
{
// Evaluate the Jacobian and sets the matrix
if(this->jacobian)
{
this->jacobianMatrix.resize(x.size(), x.size());
this->jacobian(x, this->jacobianMatrix);
if(!this->jacobian)
{// TODO: Print message or log error.
return;
}
if(this->jacobianMatrix.innerSize() > 0)
{
auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(
this->jacobianMatrix, this->f);
this->linearSolver->setSystem(linearSystem);
}
else
{
std::cerr << "Error updating Jacobian" << std::endl;
const core::SparseMatrixd &jacobianMatrix = this->jacobian(x);
if(!(jacobianMatrix.innerSize() > 0))
{// TODO: Print message and/or log error.
return;
}
auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(
jacobianMatrix, this->f);
this->linearSolver->setSystem(linearSystem);
}
//---------------------------------------------------------------------------
......
......@@ -27,65 +27,67 @@ NonLinearSolver::NonLinearSolver():
sigma(std::array<double, 2> {.1, .5}),
alpha(1e-4),
armijoMax(30)
{}
{
this->updateIterate = [](const core::Vectord &dx, core::Vectord &x)
{
x += dx;
};
}
//---------------------------------------------------------------------------
double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x)
{
/// Temporaries used in the line search
std::array<double, 3> lambda = {this->sigma[0] *this->sigma[1], 1.0, 1.0};
std::array<double, 3> fnormSqr = {this->f.squaredNorm(), 0.0, 0.0};
std::array<double, 3> fnorm = {std::sqrt(fnormSqr[0]), 0.0, 0.0};
double previousFnorm = this->nonLinearSystem.getF().norm();
std::array<double, 3> fnormSqr = {previousFnorm*previousFnorm, 0.0, 0.0};
std::array<double, 3> lambda = {this->sigma[0]*this->sigma[1], 1.0, 1.0};
/// Initialize temporaries
this->nonLinearSystem->eval(x, this->f);
fnormSqr[2] = fnormSqr[1] = this->f.squaredNorm();
fnorm[1] = std::sqrt(fnormSqr[1]);
double currentFnorm = this->nonLinearSystem->eval(x).norm();
// Exit if the function norm satisfies the Armijo-Goldstain condition
if(fnorm[1] < (1.0 - this->alpha * lambda[0]) *fnorm[0])
// Exit if the function norm satisfies the Armijo-Goldstein condition
if(currentFnorm < (1.0 - this->alpha * lambda[0])*previousFnorm)
{
return fnorm[1];
return currentFnorm;
}
// Starts armijo line serach loop
// Starts Armijo line search loop
for(size_t i = 0; i < this->armijoMax; ++i)
{
/// Update x and keep books on lambda
x -= lambda[0] * dx;
this->updateIterate(-lambda[0]*dx,x);
lambda[2] = lambda[1];
lambda[1] = lambda[0];
this->nonLinearSystem->eval(x, this->f);
fnorm[1] = this->f.norm();
currentFnorm = this->nonLinearSystem->eval(x).norm();
// Exit if the function norm satisfies the Armijo-Goldstain condition
if(fnorm[1] < (1.0 - this->alpha * lambda[0]) *fnorm[0])
// Exit if the function norm satisfies the Armijo-Goldstein condition
if(currentFnorm < (1.0 - this->alpha * lambda[0])*previousFnorm)
{
return fnorm[1];
return currentFnorm;
}
/// Update function norms
fnormSqr[2] = fnormSqr[1];
fnormSqr[1] = fnorm[1] * fnorm[1];
fnormSqr[1] = currentFnorm * currentFnorm;
/// Apply the three point parabolic model
this->parabolicModel(fnormSqr, lambda);
}
return fnorm[1];
return currentFnorm;
}
//---------------------------------------------------------------------------
void NonLinearSolver::parabolicModel(const std::array< double, int(3) > &fnorm, std::array< double, int(3) > &lambda)
void NonLinearSolver::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
/// d1 = (lambda[1] - lambda[2])*lambda[1]*lambda[2] < 0
/// if a > 0, then we have a concave up curvature and lambda defaults to:
/// lambda = sigma[0]*lambda
double a1 = lambda[2] * (fnorm[1] - fnorm[0]);
double a1 = lambda[2] * (currentFnorm - fnorm[0]);
double a2 = lambda[1] * (fnorm[2] - fnorm[0]);
double a = a1 - a2;
......@@ -112,7 +114,7 @@ void NonLinearSolver::parabolicModel(const std::array< double, int(3) > &fnorm,
}
//---------------------------------------------------------------------------
void NonLinearSolver::setSigma(const std::array< double, int(2) > &newSigma)
void NonLinearSolver::setSigma(const std::array<double,2> &newSigma)
{
this->sigma = newSigma;
}
......
......@@ -38,7 +38,8 @@
class NonLinearSolver
{
public:
using JacobianType = std::function<void(const core::Vectord &, core::SparseMatrixd &)>;
using JacobianType = std::function<const core::SparseMatrixd&(const core::Vectord&)>;
using UpdateIterateType = std::function<void(const core::Vectord&,core::Vectord&)>;
using FunctionType = SystemOfEquations::FunctionType;
public:
......@@ -114,12 +115,23 @@ public:
///
void setSystem(const FunctionType &F);
///
/// \brief Set a customized iterate update function.
///
/// \param newUpdateIterate Function used to update iterates.
///
void setUpdateIterate(const UpdateIterateType &newUpdateIterate)
{
this->updateIterate = newUpdateIterate;
}
protected:
core::Vectord f; ///< Storage for function evaluations
std::array<double, 2> sigma; ///< Safeguarding bounds for the lineseach
double alpha; ///< Parameter to measure decrease
size_t armijoMax; ///< Maximum number of steplength reductions
std::shared_ptr<SystemOfEquations> nonLinearSystem; ///< System of non-linear equations
UpdateIterateType updateIterate; ///< Update iteration function
};
#endif // SM_LINEAR_SOLVER
......@@ -37,8 +37,7 @@
class SystemOfEquations
{
public:
using FunctionType = std::function < core::Vectord &(const core::Vectord &,
core::Vectord &) >;
using FunctionType = std::function<core::Vectord&(const core::Vectord &)>;
public:
///
......@@ -58,16 +57,28 @@ public:
///
/// \brief Evaluate function at specified argument
///
/// \param x
/// \param x Value.
/// \return Function value.
///
inline void eval(const core::Vectord &x, core::Vectord &y)
inline const core::Vectord &evalF(const core::Vectord &x)
{
this->F(x,y);
this->F(x,this->f);
return this->f;
}
///
/// \brief Get the value of the function F
///
/// \return Function value.
///
inline core::Vectord &getF()
{
return this->f;
}
protected:
//!< Function associated with the system of equation to solve.
FunctionType F;
FunctionType F; ///> Function associated with the system of equation to solve.
core::Vectord f; ///> Storage for function values
};
///
......@@ -138,7 +149,7 @@ public:
///
inline core::Vectord &computeResidual(const core::Vectord &x, core::Vectord &r) const
{
r = this->rhs - this->F(x, r);
r = this->rhs - this->F(x);
return r;
}
......
......@@ -23,39 +23,52 @@
#include "BackwarEuler.h"
void BackwardEuler::solve(core::Vectord &x0, double timeStep)
void BackwardEuler::solve(const OdeSystemState &state, OdeSystemState &newState, double timeStep)
{
if(!this->F || !this->DF)
if(!this->system)
{
return;
}
auto G = [ &, this](const core::Vectord & x, core::Vectord & y) -> core::Vectord &
this->computeSystemRHS(state,newState,timeStep);
auto G = [&,this](const core::Vectord &) -> core::Vectord&
{
this->F(x, y);
y = x - (x0 + timeStep * y);
return y;
return this->rhs;
};
core::SparseMatrixd I(x0.size(), x0.size());
I.setIdentity();
auto DG = [ &, this](const core::Vectord & x, core::SparseMatrixd & J) -> void
auto DG = [&,this](const core::Vectord &) -> const core::SparseMatrixd&
{
this->DF(x, J);
J = I - timeStep * J;
this->computeSystemMatrix(newState,timeStep);
return this->systemMatrix;
};
auto updateIterate = [](const core::Vectord &dv, core::Vectord &v)
{
v += dv;
newState.getPositions() = state.getPositions() + timeStep*v;
};
auto NewtonSolver = std::make_shared<InexactNewton>();
auto solution = x0;
newState = state;
NewtonSolver->setSystem(G);
NewtonSolver->setJacobian(DG);
NewtonSolver->solve(solution);
x0 = solution;
NewtonSolver->solve(newState.getVelocities());
}
//---------------------------------------------------------------------------
void BackwardEuler::setJacobian(const NonLinearSolver::JacobianType &newJacobian)
BackwardEuler::SystemMatrixType BackwardEuler::getSystemMatrix()
{
this->DF = newJacobian;
if(!this->system)
{
return;
}
auto DG = [ &, this](const core::Vectord &x) -> void
{
auto &C = this->system->evalDFv(x);
auto &K = this->system->evalDFx(x);
};
}
......@@ -28,37 +28,78 @@
#include "Core/Vector.h"
#include "Solvers/InexactNewton.h"
#include "TimeIntegrators/OdeSystem.h"
///
/// @brief Approximates the solution to the IVP: dv/dt = F(x,dx/dt,t), x(t0) = x0 with
/// @brief Approximates the solution to the IVP: Mdv/dt = f(x,v,t), x(t0) = x0 with
/// a backward Euler scheme.
///
class BackwardEuler : public TimeIntegrator
{
public:
using SystemMatrixType = std::function<const SparseMatrixd&(const Vectord &x)>;
public:
///
/// @brief Constructor/Destructor.
/// @brief Default constructor/destructor.
///
BackwardEuler() = default;
~BackwardEuler() = default;
///
/// \brief Constructor. Takes the system describing the ODE.
///
BackwardEuler(std::shared_ptr<OdeSystem> system) : TimeIntegrator(system)
{}
///
/// @brief Perform one iteration of the Backward Euler method.
///
/// \param x Current iterate.
/// \param timeStep current timeStep.
/// \param state Current iterate.
/// \param newState New computed state.
/// \param timeStep Time step used to discretize the system.
///
void solve(core::Vectord &x0, double timeStep) override;
void solve(const OdeSystemState &state, OdeSystemState &newState, double timeStep) override;
///
/// \brief Compute and store the system matrix.
///
/// \param state Current state
/// \param newState New state
/// \param timeStep Time step used to discretize the ODE.
///
void computeSystemMatrix(const OdeSystemState &state, OdeSystemState &newState, double timeStep, bool computeRHS = true)
{
auto &M = this->system->evalMass(newState);
auto &K = this->system->evalDFv(newState);
auto &C = this->system->evalDFx(newState);
this->systemMatrix = (1.0/dt) * M;
this->systemMatrix += C;
this->systemMatrix += dt * K;
if(computeRHS)
{
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() - state.getPositions() - newState.getVelocities()*timeStep);
this->rhs -= M*(newState.getVelocities()-state.getVelocities())/timeStep;
}
}
///
/// @brief Set the Jacobian of the right hand side jacobian.
/// \brief Compute and store the right hand side of the system.
///
/// \param newJacobian Function for updating the jacobian matrix.
/// \param state Current state
/// \param newState New state
/// \param timeStep Time step used to discretize the ODE.
///
void setJacobian(const InexactNewton::JacobianType &newJacobian);
void computeSystemRHS(const OdeSystemState &state, OdeSystemState &newState, double timeStep)
{
auto &M = this->system->evalMass(newState);
auto &K = this->system->evalDFv(newState);
private:
InexactNewton::JacobianType DF;
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() - state.getPositions() - newState.getVelocities()*timeStep);
this->rhs -= M*(newState.getVelocities()-state.getVelocities())/timeStep;
}
};
#endif // BACKWAREULER_H
......@@ -28,14 +28,14 @@
#include "Core/Vector.h"
///
/// @brief Approximates the solution to the Initial value problem: dx/dt = F(x,t),
/// v(t0) = v0 using a forward Euler scheme; v(t1) = v(t0) + dtF(v(t0),t0)
/// @brief Approximates the solution to the Initial value problem: Mdv/dt = F(x,v,t),
/// v(t0) = v0 using a forward Euler scheme; Mv(t1) = Mv(t0) + dtF(v(t0),t0)
///
class ForwardEuler : public TimeIntegrator
{
public:
///
/// @brief Constructor/Destructor.
/// @brief Default constructor/destructor.
///
ForwardEuler() = default;
~ForwardEuler() = default;
......@@ -44,7 +44,7 @@ public:
/// @brief Perform one iteration of the Forward Euler method.
///
/// \param x Current iterate.
/// \param timeStep current timeStep.
/// \param timeStep Current timeStep.
///
void solve(core::Vectord &x, double timeStep) override;
......
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and 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
//
// 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.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
#include "TimeIntegrators/OdeSystem.h"
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and 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
//
// 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.