LagrangianFiniteElement.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  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; either version 2, or (at your option)
//  any later version.

//  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.  

//  $Id: LagrangianFiniteElement.hpp,v 1.5 2007/05/24 23:49:13 delpinux Exp $

#ifndef LAGRANGIAN_FINITE_ELEMENT_HPP
#define LAGRANGIAN_FINITE_ELEMENT_HPP

template <size_t numberOfDegreesOfFreedom,
          typename FiniteElementType,
          typename GivenQuadratureType>
class LagrangianFiniteElement
{
public:
  typedef GivenQuadratureType
  QuadratureType;               
  enum {
    numberOfQuadraturePoints = QuadratureType::numberOfQuadraturePoints
  };

  typedef
  TinyVector<numberOfDegreesOfFreedom>
  ElementaryVector;             
  typedef
  TinyMatrix<numberOfDegreesOfFreedom,
             numberOfDegreesOfFreedom>
  ElementaryMatrix;             
private:  
  FiniteElementType& self()
  {
    return static_cast<FiniteElementType&>(*this);
  }

protected:
  real_t __W  (const size_t& i, const size_t& j)
  {
    return self().W(i,self().integrationVertices()[j]);
  }

  real_t __dxW(const size_t& i, const size_t& j)
  {
    return self().dxW(i,self().integrationVertices()[j]);
  }

  real_t __dyW(const size_t& i, const size_t& j)
  {
    return self().dyW(i,self().integrationVertices()[j]);
  }

  real_t __dzW(const size_t& i, const size_t& j)
  {
    return self().dzW(i,self().integrationVertices()[j]);
  }

  
  TinyMatrix<numberOfDegreesOfFreedom,numberOfQuadraturePoints>
  __w;                          
  TinyMatrix<numberOfDegreesOfFreedom,numberOfQuadraturePoints>
  __dxw;                        
  TinyMatrix<numberOfDegreesOfFreedom,numberOfQuadraturePoints>
  __dyw;                        
  TinyMatrix<numberOfDegreesOfFreedom,numberOfQuadraturePoints>
  __dzw;                        
public:
  inline const real_t& W(const size_t& i, const size_t& j) const
  {
    return __w(i,j);
  }

  inline const real_t& dxW(const size_t& i, const size_t& j) const
  {
    return __dxw(i,j);
  }

  inline const real_t& dyW(const size_t& i, const size_t& j) const
  {
    return __dyw(i,j);
  }

  inline const real_t& dzW(const size_t& i, const size_t& j) const
  {
    return __dzw(i,j);
  }

public:
  template <typename ConformTransformation>
  inline void integrateWjWi(ElementaryMatrix& matElem,
                            const ConformTransformation& T) const
  {
    ElementaryMatrix tmp = 0;

    for (size_t k=0; k<numberOfQuadraturePoints; ++k) { // Loop on integration vertices
      for (size_t j=0; j<numberOfDegreesOfFreedom; ++j) {
        for (size_t i=0; i<=j; ++i) {
          tmp(i,j)
            += W(i,k) * W(j,k) * QuadratureType::instance().weight(k);
        }
      }
    }

    // for this operator, matElem is symetric.
    for (size_t j=0; j<numberOfDegreesOfFreedom; ++j)
      for (size_t i=j+1; i<numberOfDegreesOfFreedom; ++i)
        tmp(i,j) = tmp(j,i);

    matElem += tmp;
  }

  template <typename ConformTransformation>
  inline void integrateDWjWi(ElementaryMatrix& matElem,
                             const size_t& n,
                             const ConformTransformation& T) const
  {
    ElementaryMatrix tmp = 0;

    for (size_t k=0; k<numberOfQuadraturePoints; ++k) { // Loop on integration vertices
      for (size_t j=0; j<numberOfDegreesOfFreedom; ++j) {
        const real_t fj
          = dxW(j,k)*T.invJacobian(0,n)
          + dyW(j,k)*T.invJacobian(1,n)
          + dzW(j,k)*T.invJacobian(2,n);
        for (size_t i=0; i<numberOfDegreesOfFreedom; ++i) {
          tmp(i,j)
            += fj * W(i,k)  * QuadratureType::instance().weight(k);
        }
      }
    }
    matElem += tmp;
  }

  template <typename ConformTransformation>
  inline void integrateWjDWi(ElementaryMatrix& matElem,
                             const size_t& n,
                             const ConformTransformation& T) const
  {
    ElementaryMatrix tmp = 0;

    for (size_t k=0; k<numberOfQuadraturePoints; ++k) { // Loop on integration vertices
      for (size_t i=0; i<numberOfDegreesOfFreedom; ++i) {
        const real_t fi
          = dxW(i,k)*T.invJacobian(0,n)
          + dyW(i,k)*T.invJacobian(1,n)
          + dzW(i,k)*T.invJacobian(2,n);
        for (size_t j=0; j<numberOfDegreesOfFreedom; ++j) {
          tmp(i,j)
            += fi * W(j,k) * QuadratureType::instance().weight(k);
        }
      }
    }
    matElem += tmp;
  }

  template <typename ConformTransformation>
  inline void integrateDWjDWi(ElementaryMatrix& matElem,
                              const size_t& n,
                              const size_t& m,
                              const ConformTransformation& T) const
  {
    ElementaryMatrix tmp = 0;

    for (size_t k=0; k<numberOfQuadraturePoints; ++k) { // Loop on integration vertices
      for (size_t j=0; j<numberOfDegreesOfFreedom; ++j) {
        const real_t fj
          = dxW(j,k)*T.invJacobian(0,n)
          + dyW(j,k)*T.invJacobian(1,n)
          + dzW(j,k)*T.invJacobian(2,n);
        for (size_t i=0; i<numberOfDegreesOfFreedom; ++i) {
          tmp(i,j)
            += fj
            * (  dxW(i,k)*T.invJacobian(0,m)
               + dyW(i,k)*T.invJacobian(1,m)
               + dzW(i,k)*T.invJacobian(2,m) )
            * QuadratureType::instance().weight(k);
        }
      }
    }

    matElem += tmp;
  }

  template <typename ConformTransformation>
  inline void integrateWj(ElementaryVector& vectElem,
                          const ConformTransformation& T,
                          const TinyVector<numberOfQuadraturePoints, real_t>& f) const
  {
    vectElem = 0;

    for (size_t k=0; k<numberOfQuadraturePoints; ++k)
      for (size_t j=0; j<numberOfDegreesOfFreedom; ++j) {
        vectElem[j]
          += W(j,k)
          *  f[k]
          *  QuadratureType::instance().weight(k);
      }
  }

  LagrangianFiniteElement()
  {
    try {
      for (size_t i=0; i<numberOfDegreesOfFreedom; ++i) {
        for(size_t j=0; j<numberOfQuadraturePoints; ++j) {
          __w(i,j)   =   __W(i,j);
          __dxw(i,j) = __dxW(i,j);
          __dyw(i,j) = __dyW(i,j);
          __dzw(i,j) = __dzW(i,j);
        }
      }
    }
    catch(ErrorHandler e) {
      e.writeErrorMessage();    
    }
    catch(...) {
      fferr(0) << "error: Unknown exception caught!\n";
      fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
    }
  }

  virtual ~LagrangianFiniteElement()
  {
    ;
  }
};

#endif // LAGRANGIAN_FINITE_ELEMENT_HPP

---