An Object Finite Element Library

Matrix.h

Go to the documentation of this file.

/*==============================================================================

                                    O  F  E  L  I

                           Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2008 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free
   Software Foundation; Version 2 of the License.

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
   details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================

                         Definition of abstract class 'Matrix'

  ==============================================================================*/


#ifndef __MATRIX_H
#define __MATRIX_H


#include <valarray>
#include <iostream>
using std::valarray;
using std::ostream;

#include "Mesh.h"
#include "GraphOfMatrix.h"
#include "Vect.h"

namespace OFELI {

template<class T_> class DMatrix;

class Mesh;
class Element;
class Side;

#ifndef DOXYGEN_SHOULD_SKIP_THIS
typedef std::pair<size_t,size_t> RC;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifdef USE_MTL
#include "mtl/mtl.h"
typedef mtl::matrix<double,mtl::rectangle<>,mtl::compressed<int,mtl::external,mtl::index_from_one>,mtl::row_major>::type MTL_Mat;
#endif

enum MatrixType {
   SKYLINE        = 0x00800000,   
   SPARSE         = 0x01000000,   
   DIAGONAL       = 0x02000000,   
   TRIDIAGONAL    = 0x04000000,   
   SYMMETRIC      = 0x08000000,   
   UNSYMMETRIC    = 0x10000000,   
   IDENTITY       = 0x12000000    
};

enum Iteration {
   DIRECT_SOLVER    = 0,          
   CG_SOLVER        = 1,          
   CGS_SOLVER       = 2,          
   BICG_SOLVER      = 3,          
   BICG_STAB_SOLVER = 4,          
   GMRES_SOLVER     = 5,          
   QMR_SOLVER       = 6           
};

enum Preconditioner {
   IDENT_PREC       = 0,          
   DIAG_PREC        = 1,          
   ILU_PREC         = 2           
};



template<class T_>
class Matrix
{

 public:

    Matrix();
 
    Matrix(const Matrix<T_> &m);

    virtual ~Matrix() { }

    size_t getNbRows() const { return _nb_rows; }

    size_t getNbColumns() const { return _nb_cols; }

    void setPenal(double p) { _penal = p; }

    T_ getDiag(size_t k) const { return _diag[k-1]; }

    size_t getSize() const { return _size; }

    virtual void MultAdd(const Vect<T_> &x, Vect<T_> &y) const = 0;

    virtual void MultAdd(T_ a, const Vect<T_> &v, Vect<T_> &w) const = 0;

    virtual void Mult(const Vect<T_> &x, Vect<T_> &y) const = 0;

    virtual void TMult(const Vect<T_> &v, Vect<T_> &w) const = 0;

    void setDiagonal(const class Mesh &mesh);

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    virtual void setMesh(const class Mesh &mesh, size_t dof=0, size_t type=0)
    {
       _init_set_mesh(mesh,dof,type);
    }

    void _init_set_mesh(const class Mesh &mesh, size_t dof=0, size_t type=0);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    virtual void setMesh(size_t dof, const class Mesh &mesh, int code=0)=0;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    virtual void setMesh(size_t dof, size_t nb_eq, const class Mesh &mesh)=0;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

    void Assembly(const Element *el, T_ *a);

    void Assembly(const Element *el, const DMatrix<T_> &a);

    void Assembly(const Side *sd, T_ *a);

    void Assembly(const Side *sd, const DMatrix<T_> &a);

    void Prescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag=0);

    void Prescribe(const class Mesh &mesh, Vect<T_> &b, int flag=0);

    void Prescribe(size_t dof, const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag=0);

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    void Prescribe1(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag=0);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    void Prescribe(const class Mesh &mesh);
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

    void PrescribeSide(const class Mesh &mesh);

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    void Constraint(const class Mesh &mesh) { Prescribe(mesh); }
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

    virtual void Set(size_t i, size_t j, const T_ &x) = 0;

    virtual void Add(size_t i, size_t j, const T_ &x) = 0;

    virtual int Factor() = 0;

    virtual int Solve(Vect<T_> &b) = 0;

    int Solve(const Vect<T_> &b, Vect<T_> &x);

    int FactorAndSolve(Vect<T_> &b);

    unsigned long getLength() const { return _length; }

    bool isDiagonal() const { return _is_diagonal; }

    bool isFactorized() const { return _fact; }

    virtual size_t getColInd(size_t i) const;

    virtual size_t getRowPtr(size_t i) const;

    virtual T_ & operator()(size_t i, size_t j) = 0;

    virtual const T_ operator()(size_t i) const = 0;

    virtual T_ operator()(size_t i, size_t j) const = 0;

    T_ & operator[](size_t k) { return _a[k]; }

    T_ operator[](size_t k) const { return _a[k]; }

    Matrix & operator=(Matrix<T_> &m);

    Matrix & operator+=(const Matrix<T_> &m);

    Matrix & operator-=(const Matrix<T_> &m);

    Matrix & operator=(const T_ &x);

    Matrix & operator*=(const T_ &x);

    Matrix & operator+=(const T_ &x);

    Matrix & operator-=(const T_ &x);

    virtual T_ getEntry(size_t i, size_t j) const = 0;

    void setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax);

    void getPrintView(size_t &rmin, size_t &rmax, size_t &cmin, size_t &cmax) const;

#ifndef DOXYGEN_SHOULD_SKIP_THIS
// Pointer to an MTL Matrix
#ifdef USE_MTL
   MTL_Mat *A_;
#endif
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

 protected:
#ifndef DOXYGEN_SHOULD_SKIP_THIS
   size_t                  _rmin, _rmax, _cmin, _cmax;
   size_t                  _dof_type, _nb_rows, _nb_cols, _size, _dof;
   unsigned long           _length;
   T_                      _zero, _temp;
   valarray<T_>            _a, _aU, _diag;
   valarray<unsigned long> _row_ptr;
   valarray<size_t>        _col_ind, _ch;
   double                  _penal;
   bool                    _fact, _set_nodes, _set_elements, _set_sides, _is_diagonal;
#endif /* DOXYGEN_SHOULD_SKIP_THIS */
};


//                         I M P L E M E N T A T I O N                       //


template<class T_>
Matrix<T_>::Matrix()
{
   _penal = 1.e20;
   _zero = 0;
   _rmin = _cmin = 1;
   _rmax = _cmax = 10000;
}


template<class T_>
Matrix<T_>::Matrix(const Matrix<T_> &m)
{
   _zero = 0;
   _size = m._size;
   _nb_rows = m._nb_rows;
   _nb_cols = m._nb_cols;
   _length = m._length;
   _penal = m._penal;
   _ch.resize(_size);
   _diag.resize(_size);
   _ch = m._ch;
   _diag = m._diag;
   _rmin = m._rmin;
   _cmin = m._cmin;
   _rmax = m._rmax;
   _cmax = m._cmax;
}


template<class T_>
void Matrix<T_>::setDiagonal(const class Mesh &mesh)
{
   _init_set_mesh(mesh);
   _size = mesh.getNbEq();
   _a.resize(_size);
   _ch.resize(_size);
   _ch[0] = 0;
   for (size_t i=1; i<_size; i++)
      _ch[i] = i+1;
   Clear(_a);
   _fact = false;
   _dof = 0;
   _length = _nb_rows = _nb_cols = _size;
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void Matrix<T_>::_init_set_mesh(const class Mesh &mesh, size_t dof, size_t type)
{
   type = 0;
   _zero = T_(0);
   _dof_type = 0;
   if (mesh.NodesAreDOF())
      _dof_type = NODE_DOF;
   else if (mesh.SidesAreDOF())
      _dof_type = SIDE_DOF;
   else if (mesh.ElementsAreDOF())
      _dof_type = ELEMENT_DOF;
   _dof = dof;
   _size = _nb_rows = _nb_cols = mesh.getNbEq();
   if (_dof_type==NODE_DOF)
      if (_dof)
         _size = _nb_rows = _nb_cols = mesh.getNbNodes();
      else
         _size = _nb_rows = _nb_cols = mesh.getNbEq();
   else if (_dof_type==SIDE_DOF)
      if (_dof)
         _size = _nb_rows = _nb_cols = mesh.getNbSides();
      else
         _size = _nb_rows = _nb_cols = mesh.getNbEq();
   else if (_dof_type==ELEMENT_DOF)
      _size = _nb_rows = _nb_cols = mesh.getNbElements();
   else;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


template<class T_>
void Matrix<T_>::Assembly(const Element *el, T_ *a)
{
   size_t kk = 0;
   for (size_t in=1; in<=el->getNbNodes(); ++in) {
      Node *nd1 = el->getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); ++k) {
         for (size_t jn=1; jn<=el->getNbNodes(); ++jn) {
            Node *nd2 = el->getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); ++l) {
               if (nd1->getDOF(k)!=0 && nd2->getDOF(l)!=0)
                  Add(nd1->getDOF(k),nd2->getDOF(l),a[kk]);
               kk++;
            }
         }
      }
   }
}


template<class T_>
void Matrix<T_>::Assembly(const Element *el, const DMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t in=1; in<=el->getNbNodes(); ++in) {
      Node *nd1 = el->getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); ++k) {
         size_t j = 1;
         for (size_t jn=1; jn<=el->getNbNodes(); ++jn) {
            Node *nd2 = el->getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); ++l) {
               if (nd1->getDOF(k)!=0 && nd2->getDOF(l)!=0)
                  Add(nd1->getDOF(k),nd2->getDOF(l),a(i,j));
               j++;
            }
         }
         i++;
      }
   }
}


template<class T_>
void Matrix<T_>::Assembly(const Side *sd, T_ *a)
{
   size_t kk = 0;
   for (size_t in=1; in<=sd->getNbNodes(); ++in) {
      Node *nd1 = sd->getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); ++k) {
         for (size_t jn=1; jn<=sd->getNbNodes(); ++jn) {
            Node *nd2 = sd->getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); ++l) {
               if (nd1->getDOF(k)!=0 && nd2->getDOF(l)!=0)
                  Add(nd1->getDOF(k),nd2->getDOF(l),a[kk]);
               kk++;
            }
         }
      }
   }
}


template<class T_>
void Matrix<T_>::Assembly(const Side *sd, const DMatrix<T_> &a)
{
   size_t i = 1;
   for (size_t in=1; in<=sd->getNbNodes(); ++in) {
      Node *nd1 = sd->getPtrNode(in);
      for (size_t k=1; k<=nd1->getNbDOF(); ++k) {
         size_t j = 1;
         for (size_t jn=1; jn<=sd->getNbNodes(); ++jn) {
            Node *nd2 = sd->getPtrNode(jn);
            for (size_t l=1; l<=nd2->getNbDOF(); ++l) {
               if (nd1->getDOF(k)!=0 && nd2->getDOF(l)!=0)
                  Add(nd1->getDOF(k),nd2->getDOF(l),a(i,j));
               j++;
            }
         }
         i++;
      }
   }
}


template<class T_>
void Matrix<T_>::Prescribe(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag)
{
   const Node *nd;
   for (mesh.topNode(); (nd=mesh.getNode());) {
      for (size_t i=1; i<=nd->getNbDOF(); ++i) {
         if (nd->getCode(i)>0) {
            size_t k = nd->getDOF(i) - 1;
            if (flag==0) {
               _diag[k] = getEntry(k+1,k+1)*_penal;
               Set(k+1,k+1,_diag[k]);
            }
            b[k] = u[k] * _diag[k];
         }
      }
   }
}


template<class T_>
void Matrix<T_>::Prescribe(const class Mesh &mesh, Vect<T_> &b, int flag)
{
   const Node  *nd;
   for (mesh.topNode(); (nd=mesh.getNode());)
      for (size_t j=1; j<=nd->getNbDOF(); ++j)
         if (nd->getCode(j)>0) {
            size_t k = nd->getDOF(j) - 1;
            if (!flag) {
               _diag[k] = getEntry(k+1,k+1)*_penal;
               Set(k+1,k+1,_diag[k]);
            }
            b[k] = 0;
         }
}


template<class T_>
void Matrix<T_>::Prescribe(size_t dof, const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag)
{
   const Node *nd;
   for (mesh.topNode(); (nd=mesh.getNode());) {
      if (nd->getCode(dof)>0) {
         size_t k = nd->getLabel() - 1;
         if (!flag) {
            _diag[k] = getEntry(k+1,k+1)*_penal;
            Set(k+1,k+1,_diag[k]);
         }
         b[k] = u[k] * _diag[k];
      }
   }
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void Matrix<T_>::Prescribe1(const class Mesh &mesh, Vect<T_> &b, const Vect<T_> &u, int flag)
{
       const Node *nd;
       for (mesh.topNode(); (nd=mesh.getNode());) {
          for (size_t i=1; i<=nd->getNbDOF(); i++)
             if (nd->getCode(i)>0) {
                size_t k = nd->getDOF(i);
                if (!flag)
                   Add(k,k,_penal);
                b(k) = u(k)*_penal;
             }
       }
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
void Matrix<T_>::Prescribe(const class Mesh &mesh)
{
   const Node *nd;
   for (mesh.topNode(); (nd=mesh.getNode());) {
      for (size_t i=1; i<=nd->getNbDOF(); i++)
         if (nd->getCode(i)>0) {
            size_t k = nd->getDOF(i);
            Set(k,k,getEntry(k,k)*_penal);
         }
   }
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


template<class T_>
void Matrix<T_>::PrescribeSide(const class Mesh &mesh)
{
   const Side *sd;
   for (mesh.topSide(); (sd=mesh.getSide());) {
      for (size_t i=1; i<=sd->getNbDOF(); i++)
         if (sd->getCode(i)>0) {
            size_t k = sd->getDOF(i);
            Set(k,k,getEntry(k,k)*_penal);
         }
   }
}


template<class T_>
int Matrix<T_>::Solve(const Vect<T_> &b, Vect<T_> &x)
{
   x = b;
   return Solve(x);
}


template<class T_>
int Matrix<T_>::FactorAndSolve(Vect<T_> &b)
{
   Factor();
   return Solve(b);
}


template<class T_>
size_t Matrix<T_>::getColInd(size_t i) const
{
   i = 0;
   return 0;
}


template<class T_>
size_t Matrix<T_>::getRowPtr(size_t i) const
{
   i = 0;
   return 0;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator=(Matrix<T_> &m)
{
   _a = m._a;
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator+=(const Matrix<T_> &m)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] += m._a[i];
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator-=(const Matrix<T_> &m)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] -= m._a[i];
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator=(const T_ &x)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] = x;
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator*=(const T_ &x)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] *= x;
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator+=(const T_ &x)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] += x;
   return *this;
}


template<class T_>
Matrix<T_> & Matrix<T_>::operator-=(const T_ &x)
{
   for (size_t i=0; i<_length; ++i)
      _a[i] -= x;
   return *this;
}

    
template<class T_>
void Matrix<T_>::setPrintView(size_t rmin, size_t rmax, size_t cmin, size_t cmax)
{
   _rmin = rmin;
   _rmax = rmax;
   if (_rmin==0)
      _rmin = 1;
   if (_rmax>_size)
      _rmax = _size;
   _cmin = cmin;
   _cmax = cmax;
   if (_cmin==0)
      _cmin = 1;
   if (_cmax>_size)
      _cmax = _size;
}


template<class T_>
void Matrix<T_>::getPrintView(size_t &rmin, size_t &rmax, size_t &cmin, size_t &cmax) const
{
   rmin = _rmin;
   rmax = _rmax;
   cmin = _cmin;
   cmax = _cmax;
   if (rmin <= 0)
      rmin = 1;
   if (cmin <= 0)
      cmin = 1;
   if (rmax >= _nb_rows || rmax < rmin)
      rmax = _nb_rows;
   if (cmax >= _nb_cols || cmax < cmin)
      cmax = _nb_cols;
}


} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani