Commit be9707f8 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Several format fixes.

parent 38632432
Pipeline #4814 passed with stage
......@@ -34,7 +34,8 @@ BackwardGaussSeidel::BackwardGaussSeidel(
//---------------------------------------------------------------------------
void BackwardGaussSeidel::iterate(core::Vectord &x, bool updateResidual)
{
x = this->linearSystem->getRHSVector() - this->linearSystem->getStrictLowerTriangular() * x;
x = this->linearSystem->getRHSVector() -
this->linearSystem->getStrictLowerTriangular() * x;
this->linearSystem->getUpperTrianglular().solveInPlace(x);
if (updateResidual)
......@@ -61,7 +62,8 @@ void BackwardGaussSeidel::relax(core::Vectord& x)
}
//---------------------------------------------------------------------------
void BackwardGaussSeidel::setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
void BackwardGaussSeidel::
setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
{
LinearSolver::setSystem(newSystem);
}
......@@ -26,9 +26,9 @@
BackwardSOR::BackwardSOR(): weight(.9) {}
//---------------------------------------------------------------------------
BackwardSOR::BackwardSOR(const core::SparseMatrixd &A,
const core::Vectord &rhs,
const double &w): BackwardGaussSeidel(A, rhs), weight(w)
BackwardSOR::
BackwardSOR(const core::SparseMatrixd &A, const core::Vectord &rhs, const double &w)
: BackwardGaussSeidel(A, rhs), weight(w)
{}
//---------------------------------------------------------------------------
......
......@@ -23,8 +23,8 @@
#include "Solvers/ConjugateGradient.h"
ConjugateGradient::ConjugateGradient(const core::SparseMatrixd &A,
const core::Vectord &rhs) : solver(A)
ConjugateGradient::
ConjugateGradient(const core::SparseMatrixd &A, const core::Vectord &rhs) : solver(A)
{
this->linearSystem = std::make_shared<LinearSystem<core::SparseMatrixd>>(A, rhs);
this->maxIterations = rhs.size();
......@@ -37,7 +37,7 @@ ConjugateGradient::ConjugateGradient(const core::SparseMatrixd &A,
//---------------------------------------------------------------------------
void ConjugateGradient::iterate(core::Vectord &, bool)
{
std::cout << "Nothing to do\n";
// Nothing to do
}
//---------------------------------------------------------------------------
......@@ -47,8 +47,6 @@ void ConjugateGradient::solve(core::Vectord &x)
{
x = this->solver.solve(this->linearSystem->getRHSVector());
}
this->linearSystem->computeResidual(x, this->residual);
}
//---------------------------------------------------------------------------
......@@ -79,8 +77,19 @@ void ConjugateGradient::setMaximumIterations(const size_t maxIter)
}
//---------------------------------------------------------------------------
void ConjugateGradient::setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
void ConjugateGradient::
setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
{
LinearSolver<core::SparseMatrixd>::setSystem(newSystem);
this->solver.compute(this->linearSystem->getMatrix());
}
//---------------------------------------------------------------------------
void ConjugateGradient::print()
{
/// TODO: Log this
// std::cout << "#iterations: " << this->solver.iterations() << std::endl;
// std::cout << "estimated error: " << this->solver.error() << std::endl;
// std::cout << "tolerance: " << this->solver.tolerance() << std::endl;
// std::cout << "maxIterations: " << this->solver.maxIterations() << std::endl;
}
......@@ -103,6 +103,11 @@ public:
///
void setSystem(std::shared_ptr<LinearSystemType> newSystem) override;
///
/// \brief Print solver information
///
void print();
private:
///> Pointer to the Eigen's Conjugate gradient solver
Eigen::ConjugateGradient<core::SparseMatrixd> solver;
......
......@@ -25,7 +25,8 @@
#include "Solvers/SystemOfEquations.h"
#include <iostream>
DirectLinearSolver<core::Matrixd>::DirectLinearSolver(const core::Matrixd &matrix, const core::Vectord &b)
DirectLinearSolver<core::Matrixd>::
DirectLinearSolver(const core::Matrixd &matrix, const core::Vectord &b)
{
this->linearSystem = std::make_shared<LinearSystem<core::Matrixd>>(matrix, b);
this->solver.compute(matrix);
......@@ -52,14 +53,16 @@ void DirectLinearSolver<core::Matrixd>::solve(const core::Vectord &rhs, core::Ve
}
//---------------------------------------------------------------------------
void DirectLinearSolver<core::Matrixd>::setSystem(std::shared_ptr<LinearSystem<core::Matrixd>> newSystem)
void DirectLinearSolver<core::Matrixd>::
setSystem(std::shared_ptr<LinearSystem<core::Matrixd>> newSystem)
{
LinearSolver<core::Matrixd>::setSystem(newSystem);
this->solver.compute(this->linearSystem->getMatrix());
}
//---------------------------------------------------------------------------
DirectLinearSolver<core::SparseMatrixd>::DirectLinearSolver(const core::SparseMatrixd &matrix, const core::Vectord &b)
DirectLinearSolver<core::SparseMatrixd>::
DirectLinearSolver(const core::SparseMatrixd &matrix, const core::Vectord &b)
{
this->linearSystem = std::make_shared<LinearSystem<core::SparseMatrixd>>(matrix, b);
this->solver.compute(matrix);
......@@ -80,13 +83,15 @@ void DirectLinearSolver<core::SparseMatrixd>::solve(core::Vectord &x)
}
//---------------------------------------------------------------------------
void DirectLinearSolver<core::SparseMatrixd>::solve(const core::Vectord &rhs, core::Vectord &x)
void DirectLinearSolver<core::SparseMatrixd>::
solve(const core::Vectord &rhs, core::Vectord &x)
{
x = this->solver.solve(rhs);
}
//---------------------------------------------------------------------------
void DirectLinearSolver<core::SparseMatrixd>::setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
void DirectLinearSolver<core::SparseMatrixd>::
setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
{
LinearSolver<core::SparseMatrixd>::setSystem(newSystem);
this->solver.compute(this->linearSystem->getMatrix());
......
......@@ -33,7 +33,9 @@ ForwardGaussSeidel::ForwardGaussSeidel(
//---------------------------------------------------------------------------
void ForwardGaussSeidel::iterate(core::Vectord &x, bool updateResidual)
{
x = this->linearSystem->getRHSVector() - this->linearSystem->getStrictUpperTriangular() * x;
x = this->linearSystem->getRHSVector() -
this->linearSystem->getStrictUpperTriangular() * x;
this->linearSystem->getLowerTriangular().solveInPlace(x);
if (updateResidual)
......@@ -62,7 +64,8 @@ void ForwardGaussSeidel::relax(core::Vectord& x)
}
//---------------------------------------------------------------------------
void ForwardGaussSeidel::setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
void ForwardGaussSeidel::
setSystem(std::shared_ptr<LinearSystem<core::SparseMatrixd>> newSystem)
{
LinearSolver::setSystem(newSystem);
}
......@@ -26,9 +26,9 @@
ForwardSOR::ForwardSOR(): weight(.9) {}
//---------------------------------------------------------------------------
ForwardSOR::ForwardSOR(const core::SparseMatrixd &A,
const core::Vectord &rhs,
const double &w): ForwardGaussSeidel(A, rhs), weight(w)
ForwardSOR::
ForwardSOR(const core::SparseMatrixd &A, const core::Vectord &rhs, const double &w)
: ForwardGaussSeidel(A, rhs), weight(w)
{}
//---------------------------------------------------------------------------
......
......@@ -74,3 +74,9 @@ double IterativeLinearSolver::getError(const core::Vectord &x)
this->linearSystem->computeResidual(x, this->residual);
return this->residual.squaredNorm();
}
//---------------------------------------------------------------------------
void IterativeLinearSolver::print()
{
// Nothing to print
}
......@@ -79,6 +79,11 @@ public:
///
virtual double getError(const core::Vectord &x);
///
/// \brief Print solver information.
///
virtual void print();
protected:
size_t maxIterations; ///> Maximum number of iterations to be performed.
core::Vectord residual; ///> Storage for residual vector.
......
......@@ -24,11 +24,13 @@
#define SM_LINEAR_SOLVER_HPP
template<typename SystemMatrixType>
LinearSolver<SystemMatrixType>::LinearSolver() : linearSystem(nullptr), minTolerance(1.0e-6){}
LinearSolver<SystemMatrixType>::
LinearSolver() : linearSystem(nullptr), minTolerance(1.0e-6){}
//---------------------------------------------------------------------------
template<typename SystemMatrixType>
void LinearSolver<SystemMatrixType>::setSystem(std::shared_ptr<LinearSystem<SystemMatrixType>> newSystem)
void LinearSolver<SystemMatrixType>::
setSystem(std::shared_ptr<LinearSystem<SystemMatrixType>> newSystem)
{
this->linearSystem.reset();
this->linearSystem = newSystem;
......@@ -36,7 +38,8 @@ void LinearSolver<SystemMatrixType>::setSystem(std::shared_ptr<LinearSystem<Syst
//---------------------------------------------------------------------------
template<typename SystemMatrixType>
std::shared_ptr<LinearSystem<SystemMatrixType>> LinearSolver<SystemMatrixType>::getSystem() const
std::shared_ptr<LinearSystem<SystemMatrixType>> LinearSolver<SystemMatrixType>::
getSystem() const
{
return this->linearSystem;
}
......
......@@ -80,9 +80,8 @@ void NewtonMethod::solve(core::Vectord &x)
}
//---------------------------------------------------------------------------
void NewtonMethod::updateForcingTerm(const double ratio,
const double stopTolerance,
const double fnorm)
void NewtonMethod::
updateForcingTerm(const double ratio, const double stopTolerance, const double fnorm)
{
double eta = this->gamma * ratio * ratio;
double forcingTermSqr = this->forcingTerm * this->forcingTerm;
......@@ -93,17 +92,19 @@ void NewtonMethod::updateForcingTerm(const double ratio,
eta = std::max(eta, this->gamma * forcingTermSqr);
}
this->forcingTerm = std::max(std::min(eta, this->etaMax), 0.5 * stopTolerance / fnorm);
this->forcingTerm = std::max(std::min(eta, this->etaMax),
0.5 * stopTolerance / fnorm);
}
//---------------------------------------------------------------------------
void NewtonMethod::setLinearSolver(std::shared_ptr< NewtonMethod::LinearSolverType > newLinearSolver)
void NewtonMethod::
setLinearSolver(std::shared_ptr< NewtonMethod::LinearSolverType > newLinearSolver)
{
this->linearSolver = newLinearSolver;
}
//---------------------------------------------------------------------------
std::shared_ptr< NewtonMethod::LinearSolverType > NewtonMethod::getLinearSolver() const
std::shared_ptr<NewtonMethod::LinearSolverType> NewtonMethod::getLinearSolver() const
{
return this->linearSolver;
}
......
......@@ -35,7 +35,8 @@ NonLinearSolver::NonLinearSolver():
}
//---------------------------------------------------------------------------
double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x, const double previousFnorm)
double NonLinearSolver::
armijo(const core::Vectord &dx, core::Vectord &x, const double previousFnorm)
{
/// Temporaries used in the line search
std::array<double, 3> fnormSqr = {previousFnorm*previousFnorm, 0.0, 0.0};
......@@ -89,15 +90,15 @@ double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x, const
}
//---------------------------------------------------------------------------
void NonLinearSolver::parabolicModel(const std::array<double,3> &fnorm,
std::array<double,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] * (currentFnorm - fnorm[0]);
double a1 = lambda[2] * (fnorm[1] - fnorm[0]);
double a2 = lambda[1] * (fnorm[2] - fnorm[0]);
double a = a1 - a2;
......@@ -130,7 +131,7 @@ void NonLinearSolver::setSigma(const std::array<double,2> &newSigma)
}
//---------------------------------------------------------------------------
const std::array< double, int(2) > &NonLinearSolver::getSigma() const
const std::array<double,2> &NonLinearSolver::getSigma() const
{
return this->sigma;
}
......@@ -177,3 +178,10 @@ void NonLinearSolver::setSystem(const NonLinearSolver::FunctionType &F)
this->nonLinearSystem = std::make_shared<SystemOfEquations>();
this->nonLinearSystem->setFunction(F);
}
//---------------------------------------------------------------------------
void NonLinearSolver::
setUpdateIterate(const NonLinearSolver::UpdateIterateType& newUpdateIterate)
{
this->updateIterate = newUpdateIterate;
}
......@@ -62,6 +62,7 @@ void BackwardEuler::solve(const OdeSystemState &state,
this->newtonSolver.setSystem(G);
this->newtonSolver.setJacobian(DG);
this->newtonSolver.setRelativeTolerance(0);
this->newtonSolver.setLinearSolver(std::make_shared<DirectLinearSolver<core::SparseMatrixd>>());
this->newtonSolver.setLinearSolver(
std::make_shared<DirectLinearSolver<core::SparseMatrixd>>());
this->newtonSolver.solve(newState.getVelocities());
}
......@@ -89,7 +89,8 @@ void OdeSystem::computeImplicitSystemLHS(const OdeSystemState &previousState,
if(computeRHS)
{
auto &f = this->evalF(newState);
this->rhs = M * (newState.getVelocities() - previousState.getVelocities()) / timeStep;
this->rhs = M * (newState.getVelocities() -
previousState.getVelocities()) / timeStep;
this->rhs -= (f + K * (newState.getPositions() - previousState.getPositions() -
newState.getVelocities() * timeStep));
......@@ -128,7 +129,8 @@ void OdeSystem::computeImplicitSystemRHS(const OdeSystemState &state,
auto &f = this->evalF(newState);
this->rhs = M * (newState.getVelocities() - state.getVelocities()) / timeStep;
this->rhs -= (f + K * (newState.getPositions() - state.getPositions() - newState.getVelocities() * timeStep));
this->rhs -= (f + K * (newState.getPositions() - state.getPositions() -
newState.getVelocities() * timeStep));
if(this->Damping)
{
......
......@@ -62,13 +62,15 @@ void OdeSystemState::resize(const size_t size)
}
//---------------------------------------------------------------------------
void OdeSystemState::setBoundaryConditions(const std::vector< std::size_t > &boundaryConditions)
void OdeSystemState::
setBoundaryConditions(const std::vector< std::size_t > &boundaryConditions)
{
this->fixedVertices = boundaryConditions;
}
//---------------------------------------------------------------------------
void OdeSystemState::applyBoundaryConditions(core::SparseMatrixd &M, bool withCompliance) const
void OdeSystemState::
applyBoundaryConditions(core::SparseMatrixd &M, bool withCompliance) const
{
double compliance = withCompliance ? 1.0 : 0.0;
......@@ -96,7 +98,8 @@ void OdeSystemState::applyBoundaryConditions(core::SparseMatrixd &M, bool withCo
}
//---------------------------------------------------------------------------
void OdeSystemState::applyBoundaryConditions(core::Matrixd &M, bool withCompliance) const
void OdeSystemState::
applyBoundaryConditions(core::Matrixd &M, bool withCompliance) const
{
double compliance = withCompliance ? 1.0 : 0.0;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment