Linear.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Sandals project is distributed under the BSD 2-Clause License.                            *
 *                                                                                               *
 * Davide Stocco                                                               Enrico Bertolazzi *
 * University of Trento                                                     University of Trento *
 * e-mail: davide.stocco@unitn.it                             e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#pragma once
 
#ifndef SANDALS_LINEAR_SYSTEM_HXX
#define SANDALS_LINEAR_SYSTEM_HXX
 
namespace Sandals {
 
  /*\
   |   _     _
   |  | |   (_)_ __   ___  __ _ _ __
   |  | |   | | '_ \ / _ \/ _` | '__|
   |  | |___| | | | |  __/ (_| | |
   |  |_____|_|_| |_|\___|\__,_|_|
   |
  \*/

  template <Integer N, Integer M = 0>
  class Linear : public Explicit<N, M>
  {
  public:
    using Pointer  = std::shared_ptr<SemiExplicit<N, M>>; 
    using VectorF  = typename Explicit<N, M>::VectorF;    
    using MatrixJF = typename Explicit<N, M>::MatrixJF;   
    using MatrixE  = typename Explicit<N, M>::MatrixJF;   
    using MatrixA  = typename Explicit<N, M>::MatrixJF;   
    using VectorB  = typename Explicit<N, M>::VectorF;    
    using Type     = typename Explicit<N, M>::Type;       
 
  private:
    Eigen::FullPivLU<MatrixE> m_lu; 
 
  public:
    Linear() : Explicit<N, M>(Type::LINEAR, "(missing name)") {}

    Linear(std::string t_name) : Explicit<N, M>(Type::LINEAR, t_name) {}

    VectorF F(VectorF const &x, VectorF const &x_dot, Real t) const override
    {
      return this->E(t)*x_dot - this->A(t)*x - this->b(t);
    }

    MatrixJF JF_x(VectorF const &/*x*/, VectorF const &/*x_dot*/, Real t) const override {return -this->A(t);}

    MatrixJF JF_x_dot(VectorF const &/*x*/, VectorF const &/*x_dot*/, Real t) const override {return this->E(t);}

    VectorF f(VectorF const &x, Real t) const override
    {
      this->m_lu.compute(this->E(t));
      SANDALS_ASSERT(this->m_lu.rank() == N, "Sandals:Linear::f(...): singular mass matrix E(t) detected.");
      return this->m_lu.solve(this->A(t)*x + this->b(t));
    }

    MatrixJF Jf_x(VectorF const &/*x*/, Real t) const override
    {
      this->m_lu.compute(this->E(t));
      SANDALS_ASSERT(this->m_lu.rank() == N, "Sandals:Linear::Jf_x(...): singular mass matrix E(t) detected.");
      return this->m_lu.solve(this->A(t));
    }

    virtual MatrixE E(Real t) const = 0;

    virtual MatrixA A(Real t) const = 0;

    virtual VectorB b(Real t) const = 0;
 
  }; // class Linear
 
} // namespace Sandals
 
#endif // SANDALS_LINEAR_SYSTEM_HXX