Commit 467ea338 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Fix compilation errors.

parent 472c037f
Pipeline #4827 passed with stage
......@@ -158,8 +158,7 @@ public:
void generateWeigths(std::shared_ptr<SurfaceMesh> surfaceMesh,
double radius = 1.0,
const bool saveToDisk = false,
const std::string &filename = "mesh.interp"
);
const std::string &filename = "mesh.interp");
///
/// \brief Apply a tranlation to all the vertices, including attached meshes.
......
......@@ -108,7 +108,7 @@ public:
return;
}
this->odeSolver->solve(this->currentState,this->newState,dt);
this->odeSolver->solve(*this->currentState,*this->newState,dt);
this->currentState.swap(this->previousState);
this->currentState.swap(this->newState);
......@@ -151,7 +151,7 @@ protected:
core::SparseMatrixd M; ///> Mass matrix
core::SparseMatrixd C; ///> Raleigh Damping matrix
core::SparseMatrixd D; ///> Laplacian Damping matrix
core::SparseMatrixd D; ///> Raleigh Damping matrix
core::SparseMatrixd K; ///> Stiffness matrix
core::Vectord f; ///> Accumulative forces vector
......
......@@ -27,13 +27,11 @@
#include <Core/RenderDelegate.h>
SceneObjectDeformable::SceneObjectDeformable() :
OdeSystem(),
renderSecondaryMesh(false),
topologyAltered(false),
pulledVertex(-1),
timestepCounter(0),
subTimestepCounter(0),
integrationScheme(TimeIntegrator::ImplicitEuler)
subTimestepCounter(0)
{
this->setRenderDelegate(
Factory<RenderDelegate>::createConcreteClassForGroup(
......@@ -41,6 +39,10 @@ SceneObjectDeformable::SceneObjectDeformable() :
RenderDelegate::RendererType::VTK));
}
SceneObjectDeformable::~SceneObjectDeformable()
{
}
void SceneObjectDeformable::applyContactForces()
{
for(const auto &cf : this->getContactForces())
......
......@@ -31,23 +31,22 @@
#include "Mesh/SurfaceMesh.h"
#include "Core/SceneObject.h"
#include "Core/ErrorLog.h"
#include "Simulators/VegaObjectConfig.h"
#include "TimeIntegrators/TimeIntegrator.h"
#include "TimeIntegrators/BackwarEuler.h"
#include "TimeIntegrators/ForwardEuler.h"
#include "VegaObjectConfig.h"
/// \brief Base class for any scene object that is defmormable
class SceneObjectDeformable : public SceneObject, public OdeSystem
class SceneObjectDeformable : public SceneObject
{
public:
///
/// \brief Constructor
///
SceneObjectDeformable();
/// \brief Destructor
virtual ~SceneObjectDeformable() = default;
~SceneObjectDeformable();
/// \brief rest the object to inital configuration
/// and reset initial states
virtual void resetToInitialState() = 0;
/// \brief Load specified meshes
virtual void loadVolumeMesh() = 0;
......@@ -139,77 +138,6 @@ 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;
}
}
///
/// \brief Update states
///
void update(double dt)
{
if(!this->odeSolver)
{
std::cerr << "Ode solver needs to be set." << std::endl;
return;
}
this->odeSolver->solve(this->currentState,this->newState,dt);
this->currentState.swap(this->previousState);
this->currentState.swap(this->newState);
// TODO: Check state validity
}
///
/// \brief Reset the current state to the initial state
///
virtual void resetToInitialState()
{
*this->currentState = *this->initialState;
*this->previousState = *this->initialState;
}
///
/// \brief Return the current state.
///
std::shared_ptr<OdeSystemState> getCurrentState()
{
this->currentState;
}
///
/// \brief Return the previous state.
///
std::shared_ptr<OdeSystemState> getPreviousState()
{
this->previousState;
}
protected:
friend class SceneObjectDeformableRenderDelegate;
......@@ -243,20 +171,6 @@ 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;
core::SparseMatrixd M; ///> Mass matrix
core::SparseMatrixd C; ///> Damping matrix
core::SparseMatrixd K; ///> Stiffness matrix
core::Vectord f; ///> Force loads
TimeIntegrator::IntegratorType integrationScheme;
};
#endif // SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
......@@ -134,10 +134,6 @@ void VegaFemSceneObject::initialize()
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();
......
......@@ -34,7 +34,8 @@ InexactNewton::InexactNewton():
relativeTolerance(1e-6),
gamma(0.9),
etaMax(0.9),
maxIterations(50) {}
maxIterations(50),
useArmijo(true) {}
//---------------------------------------------------------------------------
void InexactNewton::solve(core::Vectord &x)
......@@ -61,11 +62,15 @@ void InexactNewton::solve(core::Vectord &x)
this->linearSolver->solve(dx);
this->updateIterate(-dx,x);
double newNorm = this->armijo(dx, x);
double newNorm = fnorm;
if(this->useArmijo)
{
newNorm = this->armijo(dx, x);
}
if(this->forcingTerm > 0 && newNorm > stopTolerance)
{
double ratio = newNorm / fnorm; // Ratio of succesive residual norms
double ratio = newNorm / fnorm; // Ratio of successive residual norms
this->updateForcingTerm(ratio, stopTolerance, fnorm);
// Reset tolerance in the linear solver according to the new forcing term
......@@ -88,7 +93,7 @@ void InexactNewton::updateForcingTerm(const double ratio,
double eta = this->gamma * ratio * ratio;
double forcingTermSqr = this->forcingTerm * this->forcingTerm;
// Saveguard to prevent the forcing term to become too small for far away iterates
// Save guard to prevent the forcing term to become too small for far away iterates
if(this->gamma * forcingTermSqr > .1)
{
eta = std::max(eta, this->gamma * forcingTermSqr);
......@@ -143,24 +148,6 @@ const NonLinearSolver::JacobianType &InexactNewton::getJacobian() const
return this->jacobian;
}
//---------------------------------------------------------------------------
void InexactNewton::setJacobianMatrix(const core::SparseMatrixd &newJacobianMatrix)
{
this->jacobianMatrix = newJacobianMatrix;
}
//---------------------------------------------------------------------------
const core::SparseMatrixd &InexactNewton::getJacobianMatrix() const
{
return this->jacobianMatrix;
}
//---------------------------------------------------------------------------
core::SparseMatrixd &InexactNewton::getJacobianMatrix()
{
return this->jacobianMatrix;
}
//---------------------------------------------------------------------------
void InexactNewton::setAbsoluteTolerance(const double aTolerance)
{
......
......@@ -28,10 +28,10 @@
#include "Solvers/LinearSolver.h"
///
/// @brief Inexact Newtom method. This version of the newton method is based on the work
/// Tim Kelly and others at NC Statu University. The algorithm is globally convergent
/// @brief Inexact Newton method. This version of the newton method is based on the work
/// Tim Kelly and others at NC State University. The algorithm is globally convergent
/// in the sense that for any initial iterate the iteration either converge to a
/// root of F or fails. Global convergence is achieved usign a line search subprocess
/// root of F or fails. Global convergence is achieved using a line search sub-process
/// and the Armijo rule.
///
class InexactNewton : public NonLinearSolver
......@@ -201,16 +201,33 @@ public:
return this->maxIterations;
}
///
/// \brief Set the useArmijo flag. If useArmijo is true a line search is performed
/// using the Armijo-Goldstein condition.
///
void setUseArmijo(const bool value)
{
this->useArmijo = value;
}
///
/// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
///
bool getUseArmijo() const
{
return this->useArmijo;
}
private:
std::shared_ptr<LinearSolverType> linearSolver;
JacobianType jacobian;
core::SparseMatrixd jacobianMatrix;
double forcingTerm;
double absoluteTolerance;
double relativeTolerance;
double gamma;
double etaMax;
size_t maxIterations;
std::shared_ptr<LinearSolverType> linearSolver; ///> Linear solver to use. Default: Conjugate gradient.
JacobianType jacobian; ///> Jacobian matrix function.
double forcingTerm; ///> Method's forcing term (Default: 0.9).
double absoluteTolerance; ///> Tolerance for the method (Default: 1.0e-3).
double relativeTolerance; ///> Relative (to the rhs) tolerance (Default: 1.0e-6).
double gamma; ///> Internal parameter used to update the forcing term (Default: 0.9).
double etaMax; ///> Maximum tolerance for the linear solver (Default: 0.9).
size_t maxIterations; ///> Maximum number of nonlinear iterations (Default: 40).
bool useArmijo; ///> True if Armijo liner search is desired (Default: true).
};
#endif // INEXACT_NEWTON_H
......@@ -38,7 +38,7 @@ NonLinearSolver::NonLinearSolver():
double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x)
{
/// Temporaries used in the line search
double previousFnorm = this->nonLinearSystem.getF().norm();
double previousFnorm = this->nonLinearSystem->getFunctionValue().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};
......
......@@ -76,7 +76,7 @@ public:
std::array<double, 3> &lambda);
///
/// \brief Set/Get Sigma. Safeguard parameter for the the linesearch method.
/// \brief Set/Get Sigma. Safeguard parameter for the the line search method.
///
/// \param newSigma New sigma parameter.
///
......@@ -84,7 +84,7 @@ public:
const std::array<double, 2> &getSigma() const;
///
/// \brief Set/Get Alpha. Parameter to measure sufficient decrease in the linerseach.
/// \brief Set/Get Alpha. Parameter to measure sufficient decrease in the line search.
///
/// \param newAlpha New alpha parameter.
///
......@@ -118,18 +118,15 @@ public:
///
/// \brief Set a customized iterate update function.
///
/// \param newUpdateIterate Function used to update iterates.
/// \param newUpdateIterate Function used to update iterates. Default: x+=dx.
///
void setUpdateIterate(const UpdateIterateType &newUpdateIterate)
{
this->updateIterate = newUpdateIterate;
}
void setUpdateIterate(const UpdateIterateType &newUpdateIterate);
protected:
core::Vectord f; ///< Storage for function evaluations
std::array<double, 2> sigma; ///< Safeguarding bounds for the lineseach
std::array<double, 2> sigma; ///< Safeguarding bounds for the line search
double alpha; ///< Parameter to measure decrease
size_t armijoMax; ///< Maximum number of steplength reductions
size_t armijoMax; ///< Maximum number of step length reductions
std::shared_ptr<SystemOfEquations> nonLinearSystem; ///< System of non-linear equations
UpdateIterateType updateIterate; ///< Update iteration function
};
......
......@@ -30,7 +30,7 @@
#include "Core/Matrix.h"
///
/// \class systemOfEquations
/// \class SystemOfEquations
///
/// \brief Base class for system of equations
///
......@@ -60,10 +60,9 @@ public:
/// \param x Value.
/// \return Function value.
///
inline const core::Vectord &evalF(const core::Vectord &x)
inline const core::Vectord &eval(const core::Vectord &x)
{
this->F(x,this->f);
return this->f;
return this->F(x);
}
///
......@@ -71,7 +70,7 @@ public:
///
/// \return Function value.
///
inline core::Vectord &getF()
inline core::Vectord &getFunctionValue()
{
return this->f;
}
......@@ -109,9 +108,10 @@ public:
A(matrix),
rhs(b)
{
this->F = [this](const core::Vectord & x, core::Vectord & y) -> core::Vectord &
this->F = [this](const core::Vectord & x) -> core::Vectord &
{
return y = this->A * x;
this->f = this->A * x;
return this->f;
};
}
......
......@@ -31,44 +31,31 @@ void BackwardEuler::solve(const OdeSystemState &state, OdeSystemState &newState,
}
this->computeSystemRHS(state,newState,timeStep);
auto G = [&,this](const core::Vectord &) -> core::Vectord&
auto G = [this](const core::Vectord &) -> core::Vectord&
{
return this->rhs;
};
auto DG = [&,this](const core::Vectord &) -> const core::SparseMatrixd&
const auto &cState = state;
auto &cNewState = newState;
auto DG = [&,this](const core::Vectord &)
{
this->computeSystemMatrix(newState,timeStep);
this->computeSystemMatrix(cState,cNewState,timeStep);
return this->systemMatrix;
};
auto updateIterate = [](const core::Vectord &dv, core::Vectord &v)
auto updateIterate = [&](const core::Vectord &dv, core::Vectord &v)
{
v += dv;
newState.getPositions() = state.getPositions() + timeStep*v;
cNewState.getPositions() = cState.getPositions() + timeStep*v;
};
auto NewtonSolver = std::make_shared<InexactNewton>();
newState = state;
NewtonSolver->setUpdateIterate(updateIterate);
NewtonSolver->setSystem(G);
NewtonSolver->setJacobian(DG);
NewtonSolver->solve(newState.getVelocities());
}
//---------------------------------------------------------------------------
BackwardEuler::SystemMatrixType BackwardEuler::getSystemMatrix()
{
if(!this->system)
{
return;
}
auto DG = [ &, this](const core::Vectord &x) -> void
{
auto &C = this->system->evalDFv(x);
auto &K = this->system->evalDFx(x);
};
}
......@@ -37,7 +37,7 @@
class BackwardEuler : public TimeIntegrator
{
public:
using SystemMatrixType = std::function<const SparseMatrixd&(const Vectord &x)>;
using SystemMatrixType = OdeSystem::MatrixFunctionType;
public:
///
......@@ -68,7 +68,10 @@ public:
/// \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)
void computeSystemMatrix(const OdeSystemState &state,
OdeSystemState &newState,
const double timeStep,
bool computeRHS = true)
{
auto &M = this->system->evalMass(newState);
auto &K = this->system->evalDFv(newState);
......@@ -81,8 +84,10 @@ public:
if(computeRHS)
{
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() - state.getPositions() - newState.getVelocities()*timeStep);
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() -
state.getPositions() - newState.getVelocities()*timeStep);
this->rhs -= M*(newState.getVelocities()-state.getVelocities())/timeStep;
state.applyBoundaryConditions(this->rhs);
}
}
......@@ -93,13 +98,17 @@ public:
/// \param newState New state
/// \param timeStep Time step used to discretize the ODE.
///
void computeSystemRHS(const OdeSystemState &state, OdeSystemState &newState, double timeStep)
void computeSystemRHS(const OdeSystemState &state,
OdeSystemState &newState,
double timeStep)
{
auto &M = this->system->evalMass(newState);
auto &K = this->system->evalDFv(newState);
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() - state.getPositions() - newState.getVelocities()*timeStep);
this->rhs = this->system->evalF(newState) + K*(newState.getPositions() -
state.getPositions() - newState.getVelocities()*timeStep);
this->rhs -= M*(newState.getVelocities()-state.getVelocities())/timeStep;
state.applyBoundaryConditions(this->rhs);
}
};
......
......@@ -23,8 +23,19 @@
#include "ForwardEuler.h"
void ForwardEuler::solve(core::Vectord &x, double timeStep)
void ForwardEuler::solve(const OdeSystemState &state, OdeSystemState &newState,double timeStep)
{
core::Vectord y(x.size());
x += timeStep * this->F(x, y);
this->computeSystemMatrix(state,newState,timeStep);
auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(
this->systemMatrix, this->rhs);
this->linearSolver->setSystem(linearSystem);
newState.getPositions() = state.getPositions() + timeStep * state.getVelocities();
this->solution = state.getVelocities();
this->linearSolver->solve(solution);
newState.getVelocities() += state.getVelocities();
}
......@@ -26,6 +26,7 @@
#include "TimeIntegrators/TimeIntegrator.h"
#include "Core/Vector.h"
#include "Solvers/ConjugateGradient.h"
///
/// @brief Approximates the solution to the Initial value problem: Mdv/dt = F(x,v,t),
......@@ -33,11 +34,15 @@
///
class ForwardEuler : public TimeIntegrator
{
public:
using LinearSolverType = LinearSolver<core::SparseMatrixd>;
public:
///
/// @brief Default constructor/destructor.
///
ForwardEuler() = default;
ForwardEuler():
linearSolver(std::make_shared<ConjugateGradient>()){}
~ForwardEuler() = default;
///
......@@ -46,8 +51,32 @@ public:
/// \param x Current iterate.
/// \param timeStep Current timeStep.
///
void solve(core::Vectord &x, 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 &, double timeStep, bool computeRHS = true)
{
auto &M = this->system->evalMass(state);
this->systemMatrix = (1.0/timeStep) * M;
state.applyBoundaryConditions(this->systemMatrix);
if(computeRHS)
{
this->rhs = this->system->evalF(state);
state.applyBoundaryConditions(this->rhs);
}
}
private:
std::shared_ptr<LinearSolverType> linearSolver; ///> Linear solver to use. Default: Conjugate gradient.
core::Vectord solution; ///> Solution to the linear solve.
};
#endif // FORWARDEULER_H
......@@ -36,8 +36,7 @@
class OdeSystem
{
public:
using MatrixFunctionType = std::function<const SparseMatrixd&(const core::Vectord &,
const core::Vectord &)>;
using MatrixFunctionType = std::function<const core::SparseMatrixd&(const OdeSystemState &s)>;
using FunctionType = std::function<const core::Vectord &(const OdeSystemState &s)>;
public:
......@@ -137,7 +136,7 @@ public:
///
/// \param newState Current positions and velocities.
///