BaseFEMDiscretization.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: BaseFEMDiscretization.hpp,v 1.14 2008/02/13 21:20:26 delpinux Exp $

#ifndef BASE_FEM_DISCRETIZATION_HPP
#define BASE_FEM_DISCRETIZATION_HPP

#include <FiniteElementTraits.hpp>
#include <Discretization.hpp>

#include <DiscretizedOperators.hpp>
#include <DegreeOfFreedomSet.hpp>

#include <PDE.hpp>
#include <PDEProblem.hpp>
#include <MassOperator.hpp>
#include <FirstOrderOperator.hpp>
#include <DivMuGrad.hpp>
#include <SecondOrderOperator.hpp>

#include <VariationalProblem.hpp>

#include <VariationalOperatorFV.hpp>
#include <VariationalOperatorFdxGV.hpp>
#include <VariationalOperatorFdxV.hpp>
#include <VariationalOperatorFgradGgradV.hpp>

#include <ConformTransformation.hpp>

#include <ErrorHandler.hpp>

template <typename GivenMeshType,
          ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class BaseFEMDiscretization
  : public Discretization
{
protected:
  typedef GivenMeshType MeshType;

  typedef typename MeshType::CellType CellType;

  typedef typename FiniteElementTraits<CellType,
                                       TypeOfDiscretization>::Type FiniteElement;

  typedef typename FiniteElement::ElementaryMatrix ElementaryMatrixType;

  typedef typename FiniteElement::ElementaryVector ElementaryVectorType;

  typedef
  typename FiniteElementTraits<CellType,
                               TypeOfDiscretization>::Transformation
  ConformTransformation;

  typedef
  typename FiniteElementTraits<CellType,
                               TypeOfDiscretization>::JacobianTransformation
  JacobianTransformation;

  const MeshType& __mesh;

  mutable ElementaryMatrixSet <ElementaryMatrixType> __eSet;

  mutable DiscretizedOperators<ElementaryMatrixType> __discretizedOperators;

  const DegreeOfFreedomSet& __degreeOfFreedomSet;

  void
  generatesElementaryVector(ElementaryVectorType& eVector,
                            const JacobianTransformation& J,
                            const TinyVector<FiniteElement::numberOfQuadraturePoints,real_t>& f) const
  {
    FiniteElement::instance().integrateWj(eVector,J,f);
  }

  void
  generatesElementaryMatrix(ElementaryMatrixSet<ElementaryMatrixType>& eSet,
                            const JacobianTransformation& J) const
  {
    if (eSet.isMassOperator()) {
      generatesElementaryMatrix(PDEOperator::massop,
                                J, eSet.massOperator());
    }
    if (eSet.isFirstOrderOperator()) {
      for (size_t i=0; i<3; ++i) {
        if (eSet.isFirstOrderUdxV(i)) {
          generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
                                    J,eSet.firstOrderOperatorUdxV(i),i);
        }
        if (eSet.isFirstOrderDxUV(i)) {
          generatesElementaryMatrix(PDEOperator::firstorderop, J,
                                    eSet.firstOrderOperatorDxUV(i),i);
        }
      }
    }
    if (eSet.isSecondOrderOperator()) {
      for (size_t i=0; i<3; ++i)
        for (size_t j=0; j<3; ++j) {
          if (eSet.isSecondOrderOperator(i,j)) {
            generatesElementaryMatrix(PDEOperator::secondorderop, J,
                                      eSet.secondOrderOperator(i,j),i,j);
          }
        }
    }
    if (eSet.isDivMuGrad()) {
      generatesElementaryMatrix(PDEOperator::divmugrad, J, eSet.divMuGrad());
    }
  }

  void
  generatesElementaryMatrix(const PDEOperator::Type operatorType,
                            const JacobianTransformation& J,
                            ElementaryMatrixType& matelem,
                            const size_t i = 0, const size_t j = 0) const
  {
    matelem = 0;
    switch(operatorType) {
    case PDEOperator::firstorderop: {
      FiniteElement::instance().integrateDWjWi(matelem,i,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::firstorderopTransposed: {
      FiniteElement::instance().integrateWjDWi(matelem,i,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::divmugrad: {
      FiniteElement::instance().integrateDWjDWi(matelem,0,0,J);
      FiniteElement::instance().integrateDWjDWi(matelem,1,1,J);
      FiniteElement::instance().integrateDWjDWi(matelem,2,2,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::secondorderop: {
      FiniteElement::instance().integrateDWjDWi(matelem,i,j,J);
      matelem *= J.jacobianDet();
      break;
    }
    case PDEOperator::massop: {
      FiniteElement::instance().integrateWjWi(matelem,J);
      matelem *= J.jacobianDet();
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unable to built elementary matrix for this operator",
                         ErrorHandler::unexpected);
    }
    }
  }

  BaseFEMDiscretization(const Problem& p,
                        const MeshType& mesh,
                        const DiscretizationType& discretisationType,
                        BaseMatrix& a,
                        BaseVector& bb,
                        const DegreeOfFreedomSet& dof)
    : Discretization(discretisationType, p, a, bb),
      __mesh(mesh),
      __eSet(problem()),
      __discretizedOperators(__eSet,problem()),
      __degreeOfFreedomSet(dof)
  {
    ;
  }

  virtual ~BaseFEMDiscretization()
  {
    ;
  }
};

#endif // BASE_FEM_DISCRETIZATION_HPP

---