hp1D/matrices.hh

Go to the documentation of this file.

#ifndef matrices1D_hh
#define matrices1D_hh

#include "basics/typedefs.hh"
#include "basics/vectorsMatricesForward.hh"
#include "basics/outputOperator.hh"
#include "integration/quadrature.hh"
#include "integration/karniadakis.hh"
#include "space/element.hh"

namespace hp1D {

  using concepts::Real;

  // **************************************************** StiffElementMatrix **

  class StiffElementMatrix : 
    public concepts::ElementMatrixBase<Real>, public concepts::OutputOperator {
  public:
    StiffElementMatrix(const uint m = 0) : concepts::ElementMatrixBase<Real>()
    { resize(m); }
      
    void resize(uint m);

    template<typename F>
    concepts::ElementMatrix<F> extract(uint m, uint n, const F factor) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    inline Real& A_(const uint i, const uint j) {
      conceptsAssert(i < this->m_, concepts::Assertion());
      conceptsAssert(j < this->n_, concepts::Assertion());
      return this->data_[i * this->n_ + j];
    }
  };

  template<typename F>
  concepts::ElementMatrix<F>
  StiffElementMatrix::extract(uint m, uint n, const F factor) const
  {
    conceptsAssert(m > 1         && n > 1,         concepts::Assertion());
    conceptsAssert(m <= this->m_ && n <= this->n_, concepts::Assertion());

    concepts::ElementMatrix<F> em(m,n); em.zeros();

    em(0,0) = em(1,1) =  factor;
    em(1,0) = em(0,1) = -factor;
    for(uint j = 2; j < std::min(m, n); ++j)
      em(j,j) = (*this)(j,j)*factor;

    return em;
  }

  // **************************************************** MassElementMatrix **

  class MassElementMatrix : 
    public concepts::ElementMatrixBase<Real>, public concepts::OutputOperator {
  public:
    MassElementMatrix(const uint m = 0) : concepts::ElementMatrixBase<Real>()
    { resize(m); }
      
    void resize(uint m);

    template<typename F>
    concepts::ElementMatrix<F> extract(uint m, uint n, const F factor) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    inline Real& M(const uint i, const uint j)
    {
      conceptsAssert(i < this->m_, concepts::Assertion());
      conceptsAssert(j < this->n_, concepts::Assertion());
      return this->data_[i * this->n_ + j];
    }
  };

  template<typename F>
  concepts::ElementMatrix<F>
  MassElementMatrix::extract(uint m, uint n, const F factor) const
  {
    conceptsAssert(m > 1         && n > 1,         concepts::Assertion());
    conceptsAssert(m <= this->m_ && n <= this->n_, concepts::Assertion());

    concepts::ElementMatrix<F> em(m,n); em.zeros();

    em(0,0) = em(1,1) = (*this)(0,0)*factor;
    em(1,0) = em(0,1) = (*this)(1,0)*factor;
    if (m > 2) {
      if (n > 2) {
        em(0,2) = em(1,2) = em(2,1) = em(2,0) = (*this)(0,2)*factor;
        em(2,2) = (*this)(2,2)*factor;
      } else {
        em(2,1) = em(2,0) = (*this)(2,0)*factor;
      }
    } else if (n > 2) {
      em(0,2) = em(1,2) = (*this)(0,2)*factor;
    }
    if (m > 3)
      em(3,0)  = (*this)(3,0)*factor;
    if (n > 3)
      em(0,3)  = (*this)(0,3)*factor;

    for(uint j = 3; j < std::min(m, n); ++j)
      em(j,j)   = (*this)(j,j)  *factor;
    for(uint j = 3; j < std::min(m+2, n); ++j)
      em(j-2,j) = (*this)(j-2,j)*factor;
    for(uint j = 3; j < std::min(m, n+2); ++j)
      em(j,j-2) = (*this)(j,j-2)*factor;

    return em;
  }

  // **************************************************** AdvectionElementMatrix **

  class AdvectionElementMatrix : 
    public concepts::ElementMatrixBase<Real>, public concepts::OutputOperator {
  public:
    AdvectionElementMatrix(const uint m = 0) 
      : concepts::ElementMatrixBase<Real>() { resize(m); }
      
    void resize(uint m);

    template<typename F>
    concepts::ElementMatrix<F> extract(uint m, uint n, const F factor) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    inline Real& C(const uint i, const uint j) {
      conceptsAssert(i < this->m_, concepts::Assertion());
      conceptsAssert(j < this->n_, concepts::Assertion());
      return this->data_[i * this->n_ + j];
    }
  };

  template<typename F>
  concepts::ElementMatrix<F>
  AdvectionElementMatrix::extract(uint m, uint n, const F factor) const
  {
    conceptsAssert(m > 1         && n > 1,         concepts::Assertion());
    conceptsAssert(m <= this->m_ && n <= this->n_, concepts::Assertion());

    concepts::ElementMatrix<F> em(m,n); em.zeros();

    em(0,0) = em(0,1) = (*this)(0,0)*factor;
    em(1,0) = em(1,1) = (*this)(1,0)*factor;
    if (n > 2)
      em(0,2) =  (*this)(0,2)*factor;
    if (m > 2)
      em(2,0) = -(*this)(0,2)*factor;
    for(uint j = 2; j < std::min(m+1, n); ++j)
      em(j-1,j) = (*this)(j-1,j)*factor;
    for(uint j = 2; j < std::min(m, n+1); ++j)
      em(j,j-1) = -(*this)(j-1,j)*factor;

    return em;
  }

  static StiffElementMatrix      stiffElemMatrix;
  static MassElementMatrix       massElemMatrix;
  static AdvectionElementMatrix  advElemMatrix;

} // namespace hp1D

#endif  // matrices1D_hh