getmatrix.h

Go to the documentation of this file.

#ifndef GETMATRIX_FILE
#define GETMATRIX_FILE
 
#include "geom.h"
#include <cassert>
#include <vector>
// #####################################################################
class Matrix
{
public:
    Matrix(int nrows, int ncols);
    //Matrix();
 
    Matrix(Matrix const &) = default;
    virtual ~Matrix();
 
    bool isSquare() const
    {
        return _nrows == _ncols;
    }
 
    int Nrows() const
    {
        return _nrows;
    }
 
    int Ncols() const
    {
        return _ncols;
    }
 
    virtual void Debug() const = 0;
 
    virtual void GetDiag(std::vector<double> &d) const = 0;
 
    std::vector<double> const &GetDiag() const
    {
        if ( 0 == _dd.size() )        // GH: better?   Nrows()>static_cast<int>(_dd.size())
        {
            _dd.resize(Nrows());
            this->GetDiag(_dd);
        }
        assert( Nrows() == static_cast<int>(_dd.size()) );
        return _dd;
    }
 
    virtual void Mult(std::vector<double> &w, std::vector<double> const &u) const = 0;
 
    virtual void Defect(
        std::vector<double> &w,
        std::vector<double> const &f, std::vector<double> const &u) const = 0;
 
    virtual int fetch(int row, int col) const = 0;
 
protected:
    int _nrows;              
    int _ncols;              
    mutable std::vector<double> _dd; 
};
 
// #####################################################################
class BisectInterpolation;  // class forward declaration
class CRS_Matrix: public Matrix
{
public:
    CRS_Matrix();
 
    CRS_Matrix(CRS_Matrix const &) = default;
 
 
    virtual ~CRS_Matrix() override;
    void GetDiag(std::vector<double> &d) const override;
    //* Extracts the diagonal elements of the sparse matrix.
    //*
    //* @return  d  vector of diagonal elements
    //*/
    //std::vector<double> const & GetDiag() const override;
 
    void Mult(std::vector<double> &w, std::vector<double> const &u) const override;
 
    void Defect(std::vector<double> &w,
                std::vector<double> const &f, std::vector<double> const &u) const override;
 
    void Debug() const override;
 
    int fetch(int row, int col) const override;
 
    bool Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const;
 
    friend void DefectRestrict(CRS_Matrix const &SK, BisectInterpolation const &P,
                               std::vector<double> &fc, std::vector<double> &ff, std::vector<double> &uf);
 
protected:
    //int _nrows;              //!< number of rows in matrix
    //int _ncols;              //!< number of columns in matrix
    int _nnz;                
    std::vector<int> _id;    
    std::vector<int> _ik;    
    std::vector<double> _sk; 
};
 
 
class FEM_Matrix: public CRS_Matrix
{
public:
    explicit FEM_Matrix(Mesh const &mesh);
 
    FEM_Matrix(FEM_Matrix const &) = default;
 
    ~FEM_Matrix() override;
 
    void Derive_Matrix_Pattern()
    {
        //Derive_Matrix_Pattern_slow();
        Derive_Matrix_Pattern_fast();
    }
    void Derive_Matrix_Pattern_fast();
    void Derive_Matrix_Pattern_slow();
 
 
    void CalculateLaplace(std::vector<double> &f);
 
    void ApplyDirichletBC(std::vector<double> const &u, std::vector<double> &f);
 
    //* Extracts the diagonal elements of the sparse matrix.
    //*
    //* @param[in,out]  d  (prellocated) vector of diagonal elements
    //*/
    //void GetDiag(std::vector<double> &d) const;   // override in MPI parallel
 
 
    void AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], std::vector<double> &f);
 
 
private:
    Mesh const &_mesh;       
 
};
 
 
//* Prolongation matrix in CRS format (compressed row storage; also named CSR),
//* see an <a href="https://en.wikipedia.org/wiki/Sparse_matrix">introduction</a>.
//*
//* The prolongation is applied for each node from the coarse mesh to the fine mesh and
//* is derived only geometrically (no operator weighted prolongation).
//*/
//class Prolongation: public CRS_Matrix
//{
//public:
//* Intializes the CRS matrix structure from the given discetization in @p mesh.
//*
//* The sparse matrix pattern is generated but the values are 0.
//*
//* @param[in] cmesh coarse mesh
//* @param[in] fmesh fine mesh
//*
//* @warning A reference to the discretizations @p fmesh  @p cmesh are stored inside this class.
//*          Therefore, changing these meshes outside requires also
//*          to call method @p Derive_Matrix_Pattern explicitely.
//*
//* @see Derive_Matrix_Pattern
//*/
//Prolongation(Mesh const & cmesh, Mesh const & fmesh);
 
//* Destructor.
//*/
//~Prolongation() override
//{}
 
//* Generates the sparse matrix pattern and overwrites the existing pattern.
//*
//* The sparse matrix pattern is generated but the values are 0.
//*/
//void Derive_Matrix_Pattern() override;
 
//private:
//Mesh const & _cmesh;      //!< reference to coarse discretization
//Mesh const & _fmesh;      //!< reference to fine discretization
//};
 
// *********************************************************************
 
 
class BisectInterpolation: public Matrix
{
public:
    explicit BisectInterpolation(std::vector<int> const &fathers);
    BisectInterpolation();
 
    BisectInterpolation(BisectInterpolation const &) = default;
 
    ~BisectInterpolation() override;
 
    void GetDiag(std::vector<double> &d) const override;
    //* Extracts the diagonal elements of the sparse matrix.
    //*
    //* @return  d  vector of diagonal elements
    //*/
    //std::vector<double> const & GetDiag() const override;
 
    void Mult(std::vector<double> &wf, std::vector<double> const &uc) const override;
 
    void MultT(std::vector<double> const &wf, std::vector<double> &uc) const;
 
    void Defect(std::vector<double> &w,
                std::vector<double> const &f, std::vector<double> const &u) const override;
 
    void Debug() const override;
 
    int fetch(int row, int col) const override;
 
    friend void DefectRestrict(CRS_Matrix const &SK, BisectInterpolation const &P,
                               std::vector<double> &fc, std::vector<double> &ff, std::vector<double> &uf);
 
protected:
    std::vector<int> _iv;     
    std::vector<double> _vv;  
};
 
class BisectIntDirichlet: public BisectInterpolation
{
public:
    BisectIntDirichlet()
        : BisectInterpolation()
    {}
 
    BisectIntDirichlet(std::vector<int> const &fathers, std::vector<int> const &idxc_dir);
 
    BisectIntDirichlet(BisectIntDirichlet const &) = default;
 
 
    ~BisectIntDirichlet() override;
};
 
 
 
// *********************************************************************
 
void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]);
 
void AddElem(int const ial[3], double const ske[3][3], double const fe[3],
             int const id[], int const ik[], double sk[], double f[]);
 
 
 
#endif