ODESystem.h

Go to the documentation of this file.

 
#pragma once
 
#include "EquationSystem.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/IndexValueVector.h"
#include "Types.h"
 
namespace NumLib
{

template <ODESystemTag ODETag, NonlinearSolverTag NLTag>
class ODESystem;

template <>
class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                NonlinearSolverTag::Picard> : public EquationSystem
{
public:
    static const ODESystemTag ODETag =
        ODESystemTag::FirstOrderImplicitQuasilinear;

    virtual void preAssemble(const double t, double const dt,
                             GlobalVector const& x) = 0;

    virtual std::vector<GlobalIndexType>
    getIndicesOfResiduumWithoutInitialCompensation() const
    {
        return {};
    }
 
    virtual void setReleaseNodalForces(GlobalVector const* /*r_neq*/,
                                       int const /*process_id*/)
    {
    }

    virtual void assemble(const double t, double const dt,
                          std::vector<GlobalVector*> const& x,
                          std::vector<GlobalVector*> const& x_prev,
                          int const process_id, GlobalMatrix& M,
                          GlobalMatrix& K, GlobalVector& b) = 0;
 
    using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;

    virtual std::vector<NumLib::IndexValueVector<Index>> const*
    getKnownSolutions(double const /*t*/, GlobalVector const& /*x*/,
                      int const /*process_id*/) const
    {
        return nullptr;  // by default there are no known solutions
    }
 
    virtual void updateConstraints(GlobalVector& /*lower*/,
                                   GlobalVector& /*upper*/,
                                   int const /*process_id*/){};

    virtual bool shouldLinearSolverComputeOnlyUponTimestepChange() const
    {
        return false;
    }
};

template <>
class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                NonlinearSolverTag::Newton>
    : public ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                       NonlinearSolverTag::Picard>
{
public:
    void preAssemble(const double t, double const dt,
                     GlobalVector const& x) override = 0;

    virtual void assembleWithJacobian(const double t, double const dt,
                                      std::vector<GlobalVector*> const& x,
                                      std::vector<GlobalVector*> const& x_prev,
                                      int const process_id, GlobalVector& b,
                                      GlobalMatrix& Jac) = 0;
};

}  // namespace NumLib