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

#ifndef FEM_DISCRETIZATION_HPP
#define FEM_DISCRETIZATION_HPP

#include <ElementaryMatrixSet.hpp>

#include <BaseFEMDiscretization.hpp>

#include <Mesh.hpp>
#include <Structured3DMesh.hpp>

#include <Timer.hpp>

#include <DoubleHashedMatrix.hpp>
#include <UnAssembledMatrix.hpp>

#include <ErrorHandler.hpp>

#include <FEMFunctionBase.hpp>

#include <SolverInformationCenter.hpp>

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

  typedef typename MeshType::CellType CellType;

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

  typedef typename FiniteElement::QuadratureType QuadratureType;

  typedef typename FiniteElement::ElementaryMatrix ElementaryMatrixType;

  typedef typename FiniteElement::ElementaryVector ElementaryVectorType;

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

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

public:
  void assembleMatrix()
  {
    const ScalarDegreeOfFreedomPositionsSet& dofPositions = this->__degreeOfFreedomSet.positionsSet(0);
  
    switch ((this->__A).type()) {
    case BaseMatrix::doubleHashedMatrix: {
      Timer t;
      t.start();

      ffout(2) << "- assembling matrix\n";

      DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(this->__A);

      for(typename MeshType::const_iterator icell(this->__mesh);
          not(icell.end()); ++icell) {
        const CellType& K = *icell;

        const size_t cellNumber = icell.number();
        TinyVector<FiniteElement::numberOfDegreesOfFreedom,
                   size_t> index;

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          index[l] = dofPositions(cellNumber,l);
        }

        ConformTransformation T(K);
        JacobianTransformation J(T);
        generatesElementaryMatrix(this->__eSet,J);

        TinyVector<3> massCenter;
        T.value(FiniteElement::massCenter(), massCenter);

        for(typename DiscretizedOperators<ElementaryMatrixType>
              ::iterator iOperator = this->__discretizedOperators.begin();
            iOperator != this->__discretizedOperators.end(); ++iOperator) {
          const real_t coef = ((*iOperator).second)(massCenter);

          const ElementaryMatrixType& elementaryMatrix
            = (*(*iOperator).first);

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I // Index of Ith value
              = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t J // Index of Jth value
                = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
              A(I,J) += coef*elementaryMatrix(l,m);
            }
          }
        }
      }

      ffout(2) << "- assembling matrix: done";
      t.stop();
      ffout(3) << " [cost: " << t << ']';
      ffout(2) << '\n';
      break;
    }
    case BaseMatrix::unAssembled: {
      UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(this->__A);
      A.setDiscretization(this);
      break;
    }

    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected matrix type",
                       ErrorHandler::unexpected);
    }
    }
  }

  void timesX(const BaseVector& X, BaseVector& Z) const
  {
    const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
    z=0;

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    for(typename MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;

      ConformTransformation T(K);
      JacobianTransformation J(T);
      generatesElementaryMatrix(this->__eSet,J);

      const size_t cellNumber = icell.number();

      TinyVector<FiniteElement::numberOfDegreesOfFreedom,size_t> index;
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        const real_t coef = ((*iOperator).second)(massCenter);
        const ElementaryMatrixType& elementaryMatrix
          = (*(*iOperator).first);

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I1 // Index of Ith value
            = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
          if (not((*this->__dirichletList)[I1])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t I2 // Index of Jth value
                = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
              if (not((*this->__dirichletList)[I2])) {
                z[I1] += coef*elementaryMatrix(l,m)*x[I2];
              }
            }
          }
        }
      }
    }
  }

  void transposedTimesX(const BaseVector& X, BaseVector& Z) const
  {
    const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
    z=0;

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    for(typename MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;

      ConformTransformation T(K);
      JacobianTransformation J(T);
      generatesElementaryMatrix(this->__eSet,J);

      const size_t cellNumber = icell.number();

      TinyVector<FiniteElement::numberOfDegreesOfFreedom,size_t> index;
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        const real_t coef = ((*iOperator).second)(massCenter);
        const ElementaryMatrixType& elementaryMatrix
          = (*(*iOperator).first);

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I1 // Index of Ith value
            = __degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
          if (not((*this->__dirichletList)[I1])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t I2 // Index of Jth value
                = __degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
              if (not((*this->__dirichletList)[I2])) {
                z[I2] += coef*elementaryMatrix(l,m)*x[I1];
              }
            }
          }
        }
      }
    }
  }

  void getDiagonal(BaseVector& Z) const
  {
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    for(typename MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;
      const size_t cellNumber = icell.number();

      ConformTransformation T(K);
      JacobianTransformation J(T);
      generatesElementaryMatrix(this->__eSet,J);

      size_t index[FiniteElement::numberOfDegreesOfFreedom];
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        // Only diagonal opertors have to be considerate
        if (((*iOperator).second).i() == ((*iOperator).second).j()) {
          const real_t coef = ((*iOperator).second)(massCenter);

          const ElementaryMatrixType& elementaryMatrix
            = (*(*iOperator).first);

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I // Index of Ith value
              = this->__degreeOfFreedomSet(((*iOperator).second).i(),
                                           index[l]);
            z[I] += coef*elementaryMatrix(l,l);
          }
        }
      }
    }
  }

  void assembleSecondMember()
  {
    Timer t;
    t.start();

    ffout(2) << "- assembling second member\n";
    ElementaryVectorType eVect;
    Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));
    b = 0;

    const QuadratureType& referenceQuadrature = QuadratureType::instance();

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    switch(this->problem().type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem
        = static_cast<const PDESystem&>(this->problem());

      for(typename MeshType::const_iterator icell(this->__mesh);
          not(icell.end()); ++icell) {
        const CellType& K = *icell;
        eVect = 0;
        ConformTransformation T(K);
        JacobianTransformation J(T);

        TinyVector<QuadratureType::numberOfQuadraturePoints,
                   TinyVector<3> > quadrature;
        for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
          T.value(referenceQuadrature[l], quadrature[l]);
        }

        const size_t cellNumber = icell.number();

        TinyVector<FiniteElement::numberOfDegreesOfFreedom, int> index;
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          index[l] = dofPositions(cellNumber,l);
        }

        for (size_t i=0; i<this->problem().numberOfUnknown(); ++i) {
          TinyVector<FiniteElement::numberOfQuadraturePoints, real_t> f;
          const ScalarFunctionBase& F = pdeSystem.secondMember(i);

          for (size_t l=0; l<FiniteElement::numberOfQuadraturePoints; ++l)
            f[l] = F(quadrature[l]);

          generatesElementaryVector(eVect,J,f);

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
            const size_t dof = __degreeOfFreedomSet(i,index[l]);
            b[dof] += eVect[l]*J.jacobianDet();
          }
        }
      }
      break;
    }
    case Problem::variational: {
      const VariationalProblem& variationalProblem
        =  dynamic_cast<const VariationalProblem&>(this->problem());

      if (variationalProblem.hasJumpOrMean()) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "cannot discretize problems with jumps or means in FEM !",
                           ErrorHandler::normal);
      }

      for (VariationalProblem::linearOperatorConst_iterator i
             = variationalProblem.beginLinearOperator();
           i != variationalProblem.endLinearOperator(); ++i) {
        switch ((*(*i)).type()) {
        case VariationalLinearOperator::FV: {
          const VariationalOperatorFV& fv
            = dynamic_cast<const VariationalOperatorFV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          for(typename MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;
            eVect = 0;
            ConformTransformation T(K);
            JacobianTransformation J(T);

            TinyVector<QuadratureType::numberOfQuadraturePoints,
                       TinyVector<3> > quadrature;
            for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
              T.value(referenceQuadrature[l], quadrature[l]);
            }

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            TinyVector<FiniteElement::numberOfQuadraturePoints,real_t> f;

            for (size_t l=0; l<FiniteElement::numberOfQuadraturePoints; ++l)
              f[l] = F(quadrature[l]);

            generatesElementaryVector(eVect,J,f);

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof = __degreeOfFreedomSet(testFunctionNumber,index[l]);
              b[dof] += eVect[l]*J.jacobianDet();
            }
          }
          break;
        }
        case VariationalLinearOperator::FdxGV: {
          const VariationalOperatorFdxGV& fv
            = dynamic_cast<const VariationalOperatorFdxGV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const ScalarFunctionBase& G = fv.g();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          Vector<real_t> gValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            gValues[i] = G(dofPositions.vertex(i));
          }


          for (typename MeshType::const_iterator icell(this->__mesh);
               not(icell.end()); ++icell) {
            const CellType& K = *icell;

            ConformTransformation T(K);
            JacobianTransformation J(T);

            ElementaryMatrixType edxUV = 0;
            generatesElementaryMatrix(PDEOperator::firstorderop,
                                      J, edxUV, fv.number());

            TinyVector<3> massCenter;
            T.value(FiniteElement::massCenter(), massCenter);

            const real_t fValue = F(massCenter);

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = __degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += edxUV(l,m)*gValues[index[m]]*fValue;
              }
            }
          }
          break;
        }
        case VariationalLinearOperator::FdxV: {
          const VariationalOperatorFdxV& fv
            = dynamic_cast<const VariationalOperatorFdxV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          Vector<real_t> fValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            fValues[i] = F(dofPositions.vertex(i));
          }

          for(typename MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;

            ConformTransformation T(K);
            JacobianTransformation J(T);

            ElementaryMatrixType eUdxV = 0;
            generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
                                      J, eUdxV, fv.number());

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += eUdxV(l,m)*fValues[index[m]];
              }
            }
          }
          break;
        }
        case VariationalLinearOperator::FgradGgradV: {
          const VariationalOperatorFgradGgradV& fv
            = dynamic_cast<const VariationalOperatorFgradGgradV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const ScalarFunctionBase& G = fv.g();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          Vector<real_t> gValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            gValues[i] = G(dofPositions.vertex(i));
          }

          for(typename MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;

            ConformTransformation T(K);
            JacobianTransformation J(T);

            ElementaryMatrixType e = 0;
            generatesElementaryMatrix(PDEOperator::divmugrad,
                                      J, e);

            TinyVector<3> massCenter;
            T.value(FiniteElement::massCenter(), massCenter);

            const real_t fValue = F(massCenter);

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = __degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += e(l,m)*gValues[index[m]]*fValue;
              }
            }
          }
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unknown variational operator",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unknown problem type",
                         ErrorHandler::unexpected);
    }
    }
    ffout(2) << "- assembling second member: done";
    t.stop();
    ffout(3) << " [cost: " << t << ']';
    ffout(2) << '\n';
  }

public:

  FEMDiscretization(const Problem& p,
                    const MeshType& m,
                    const DiscretizationType& discretizationType,
                    BaseMatrix& a,
                    BaseVector& bb,
                    const DegreeOfFreedomSet& dof)
    : BaseFEMDiscretization<MeshType,
                            TypeOfDiscretization>(p, m, discretizationType, a, bb, dof)
  {
    SolverInformationCenter::instance().pushMesh(&m);
    SolverInformationCenter::instance().pushDiscretizationType(&(this->__discretizationType));
  }

  ~FEMDiscretization()
  {
    SolverInformationCenter::instance().pop();
  }
};


template <>
template <ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class FEMDiscretization<Structured3DMesh,TypeOfDiscretization>
  : public BaseFEMDiscretization<Structured3DMesh,
                                 TypeOfDiscretization>
{
public:
  typedef Structured3DMesh MeshType;

  typedef MeshType::CellType CellType;

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

  typedef typename FiniteElement::QuadratureType QuadratureType;

  typedef typename FiniteElement::ElementaryMatrix ElementaryMatrixType;

  typedef typename FiniteElement::ElementaryVector ElementaryVectorType;

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

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

public:
  void assembleMatrix()
  {
    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
    
    switch (this->__A.type()) {
    case BaseMatrix::doubleHashedMatrix: {
      Timer t;
      t.start();

      ffout(2) << "- assembling matrix\n";

      DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(this->__A);

      MeshType::const_iterator icell0(this->__mesh);
      ConformTransformation T0(*icell0);
      JacobianTransformation J(T0);
      generatesElementaryMatrix(this->__eSet,J);

      for(MeshType::const_iterator icell(this->__mesh);
          not(icell.end()); ++icell) {
        const CellType& K = *icell;
        ConformTransformation T(K);
        TinyVector<3> massCenter;
        T.value(FiniteElement::massCenter(), massCenter);

        TinyVector<FiniteElement::numberOfDegreesOfFreedom,
          size_t> index;
        const size_t cellNumber = icell.number();

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          index[l] = dofPositions(cellNumber,l);
        }
        for(typename DiscretizedOperators<ElementaryMatrixType>
              ::iterator iOperator = this->__discretizedOperators.begin();
            iOperator != this->__discretizedOperators.end(); ++iOperator) {

          const real_t coef = ((*iOperator).second)(massCenter);

          const ElementaryMatrixType& elementaryMatrix
            = (*(*iOperator).first);

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I1 // Index of Ith value
              = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t I2 // Index of Jth value
                = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);

              A(I1,I2) += coef*elementaryMatrix(l,m);
            }
          }
        }
      }
      ffout(2) << "- assembling matrix: done";
      t.stop();
      ffout(3) << " [cost: " << t << ']';
      ffout(2) << '\n';

      break;
    }
    case BaseMatrix::unAssembled: {
      UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(this->__A);
      A.setDiscretization(this);
      break;
    }

    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected matrix type",
                         ErrorHandler::unexpected);
    }
    }
  }

  void timesX(const BaseVector& X, BaseVector& Z) const
  {
    const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);

    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
    z=0;

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    MeshType::const_iterator icell(this->__mesh);
    ConformTransformation T0(*icell);
    JacobianTransformation J(T0);
    generatesElementaryMatrix(this->__eSet,J);

    for(MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;
      ConformTransformation T(K);

      const size_t cellNumber = icell.number();

      TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        const real_t coef = ((*iOperator).second)(massCenter);

        const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);

#warning can optimize this loops by precomputing (i and) j
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          // Number of the degree of freedom
          const size_t I1 = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
          if (not((*this->__dirichletList)[I1])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              // Number of the degree of freedom
              const size_t I2 = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
              if (not((*this->__dirichletList)[I2])) {
                z[I1] += coef*elementaryMatrix(l,m)*x[I2];
              }
            }
          }
        }
      }
    }
  }

  void transposedTimesX(const BaseVector& X, BaseVector& Z) const
  {
    const Vector<real_t>& x = dynamic_cast<const Vector<real_t>&>(X);

    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);
    z=0;

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    MeshType::const_iterator icell(this->__mesh);
    ConformTransformation T0(*icell);
    JacobianTransformation J(T0);
    generatesElementaryMatrix(this->__eSet,J);

    for(MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;
      ConformTransformation T(K);

      const size_t cellNumber = icell.number();

      TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        const real_t coef = ((*iOperator).second)(massCenter);

        const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);

#warning can optimize this loops by precomputing (i and) j
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          // Number of the degree of freedom
          const size_t I1 = this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
          if (not((*this->__dirichletList)[I1])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              // Number of the degree of freedom
              const size_t I2 = this->__degreeOfFreedomSet(((*iOperator).second).j(), index[m]);
              if (not((*this->__dirichletList)[I2])) {
                z[I2] += coef*elementaryMatrix(l,m)*x[I1];
              }
            }
          }
        }
      }
    }
  }

  void getDiagonal(BaseVector& Z) const
  {
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    MeshType::const_iterator icell0(this->__mesh);
    ConformTransformation T0(*icell0);
    JacobianTransformation J(T0);
    generatesElementaryMatrix(this->__eSet,J);

    for(MeshType::const_iterator icell(this->__mesh);
        not(icell.end()); ++icell) {
      const CellType& K = *icell;
      ConformTransformation T(K);

      const size_t cellNumber = icell.number();

      TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
      for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
        index[l] = dofPositions(cellNumber,l);
      }

      TinyVector<3> massCenter;
      T.value(FiniteElement::massCenter(), massCenter);

      for(typename DiscretizedOperators<ElementaryMatrixType>
            ::iterator iOperator = this->__discretizedOperators.begin();
          iOperator != this->__discretizedOperators.end(); ++iOperator) {
        // Only diagonal opertors have to be considerate
        if (((*iOperator).second).i() == ((*iOperator).second).j()) {
          const real_t coef = ((*iOperator).second)(massCenter);
        
          const ElementaryMatrixType& elementaryMatrix = (*(*iOperator).first);
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I = // Index of Ith value
              this->__degreeOfFreedomSet(((*iOperator).second).i(), index[l]);
            z[I] += coef*elementaryMatrix(l,l);
          }
        }
      }
    }
  }

  void assembleSecondMember()
  {
    Timer t;
    t.start();

    ffout(2) << "- assembling second member\n";

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    ElementaryVectorType eVect;
    Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));
    b = 0;

    const QuadratureType& referenceQuadrature = QuadratureType::instance();

    switch (this->problem().type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem
        = static_cast<const PDESystem&>(this->problem());

      for(MeshType::const_iterator icell(this->__mesh);
            not(icell.end()); ++icell) {
        const CellType& K = *icell;
        eVect = 0;
        ConformTransformation T(K);
        JacobianTransformation J(T);
        TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
        TinyVector<QuadratureType::numberOfQuadraturePoints,
                   TinyVector<3> > quadrature;

        for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
          T.value(referenceQuadrature[l], quadrature[l]);
        }

        const size_t cellNumber = icell.number();
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          index[l] = this->__degreeOfFreedomSet.positionsSet(0)(cellNumber,l);
        }

        for (size_t i=0; i<pdeSystem.numberOfUnknown(); ++i) {
          const ScalarFunctionBase& F = pdeSystem.secondMember(i);

          TinyVector<QuadratureType::numberOfQuadraturePoints> f;

          for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
            f[l] = F(quadrature[l]);
          }

          generatesElementaryVector(eVect,J,f);

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
            const size_t dof = this->__degreeOfFreedomSet(i,index[l]);
            b[dof] += eVect[l]*J.jacobianDet();
          }
        }
      }
      break;
    }
    case Problem::variational: {
      const VariationalProblem& variationalProblem
        =  dynamic_cast<const VariationalProblem&>(this->problem());

      if (variationalProblem.hasJumpOrMean()) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "cannot discretize problems with jumps or means in FEM !",
                           ErrorHandler::normal);
      }

      for (VariationalProblem::linearOperatorConst_iterator i
             = variationalProblem.beginLinearOperator();
           i != variationalProblem.endLinearOperator(); ++i) {
        switch ((*(*i)).type()) {
        case VariationalLinearOperator::FV: {
          const VariationalOperatorFV& fv
            = dynamic_cast<const VariationalOperatorFV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          for(MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;
            eVect = 0;
            ConformTransformation T(K);
            JacobianTransformation J(T);

            // computes local quadrature vertex
            TinyVector<QuadratureType::numberOfQuadraturePoints,
                       TinyVector<3, real_t> > quadrature;
            for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
              T.value(referenceQuadrature[l], quadrature[l]);
            }

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            ElementaryVectorType f;

            for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
              f[l] = F(quadrature[l]);
            }

            generatesElementaryVector(eVect,J,f);

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);
              b[dof] += eVect[l]*J.jacobianDet();
            }
          }
          break;
        }
        case VariationalLinearOperator::FdxGV: {
          const VariationalOperatorFdxGV& fv
            = dynamic_cast<const VariationalOperatorFdxGV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const ScalarFunctionBase& G = fv.g();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          Vector<real_t> gValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            gValues[i] = G(dofPositions.vertex(i));
          }

          MeshType::const_iterator icell0(this->__mesh);
          ConformTransformation T0(*icell0);
          JacobianTransformation J(T0);

          ElementaryMatrixType edxUV = 0;
          generatesElementaryMatrix(PDEOperator::firstorderop,
                                    J, edxUV, fv.number());

          for(MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;

            ConformTransformation T(K);

            TinyVector<3> massCenter;
            T.value(FiniteElement::massCenter(), massCenter);
            const real_t fValue = F(massCenter);

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += edxUV(l,m)*gValues[index[m]]*fValue;
              }
            }
          }
          break;
        }
        case VariationalLinearOperator::FdxV: {
          const VariationalOperatorFdxV& fv
            = dynamic_cast<const VariationalOperatorFdxV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const size_t testFunctionNumber = fv.testFunctionNumber();

          Vector<real_t> fValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            fValues[i] = F(dofPositions.vertex(i));
          }

          MeshType::const_iterator icell0(this->__mesh);
          ConformTransformation T0(*icell0);
          JacobianTransformation J(T0);

          ElementaryMatrixType eUdxV = 0;
          generatesElementaryMatrix(PDEOperator::firstorderopTransposed,
                                    J, eUdxV, fv.number());

          for(MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;

            ConformTransformation T(K);

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = this->__degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += eUdxV(l,m)*fValues[index[m]];
              }
            }
          }
          break;
        }
        case VariationalLinearOperator::FgradGgradV: {
          const VariationalOperatorFgradGgradV& fv
            = dynamic_cast<const VariationalOperatorFgradGgradV&>(*(*i));

          const ScalarFunctionBase& F = fv.f();
          const ScalarFunctionBase& G = fv.g();

          const size_t testFunctionNumber = fv.testFunctionNumber();

          if (G.type() == ScalarFunctionBase::femfunction) {
            const FEMFunctionBase& gFEM = dynamic_cast<const FEMFunctionBase&>(G);
            if (gFEM.baseMesh() != & this->__mesh) {
              if (gFEM.baseMesh()->type() == Mesh::cartesianHexahedraMesh) {
                const MeshType& quadratureMesh
                  = dynamic_cast<const MeshType&>(*gFEM.baseMesh());
                for (MeshType::const_iterator iQuadratureCell(quadratureMesh);
                     not(iQuadratureCell.end()); ++iQuadratureCell) {
                  const CellType& KQuadrature = *iQuadratureCell;
                  ConformTransformation TQuadrature(KQuadrature);
                  
                  for (size_t l=0; l<QuadratureType::numberOfQuadraturePoints; ++l) {
                    TinyVector<3, real_t> quadrature;
                    TQuadrature.value(referenceQuadrature[l], quadrature);

                    MeshType::const_iterator icell = this->__mesh.find(quadrature);
                    if (icell.end()) continue;

                    const CellType& K = *icell;
                    if (not K.isFictitious()) {
                      const TinyVector<3,real_t> gradGValue = gFEM.gradient(quadrature);
                      const real_t fValue = F(quadrature);

                      ConformTransformation T(K);
                      JacobianTransformation J(T);

                      TinyVector<3,real_t> Xhat;
                      T.invertT(quadrature,Xhat);

                      const size_t cellNumber = icell.number();

                      TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
                      for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                        index[m] = dofPositions(cellNumber,m);
                      }

                      for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                        TinyVector<3,real_t> gradV;

                        gradV[0] = FiniteElement::instance().dxW(m,Xhat)
                          * J.invJacobian(0,0);
                        gradV[1] = FiniteElement::instance().dyW(m,Xhat)
                          * J.invJacobian(1,1);
                        gradV[2] = FiniteElement::instance().dzW(m,Xhat)
                          * J.invJacobian(2,2);

                        const size_t dof
                          = this->__degreeOfFreedomSet(testFunctionNumber,index[m]);

                        b[dof] += KQuadrature.volume()*QuadratureType::instance().weight(l)*fValue*(gradGValue*gradV);
                      }
                    }
                  }
                }
                break;
              }
            }
          }

          Vector<real_t> gValues (dofPositions.number());

          for (size_t i=0; i<dofPositions.number(); ++i) {
            gValues[i] = G(dofPositions.vertex(i));
          }

          MeshType::const_iterator icell0(this->__mesh);
          ConformTransformation T0(*icell0);
          JacobianTransformation J(T0);

          ElementaryMatrixType e = 0;
          generatesElementaryMatrix(PDEOperator::divmugrad,
                                    J, e);

          for(MeshType::const_iterator icell(this->__mesh);
              not(icell.end()); ++icell) {
            const CellType& K = *icell;
            ConformTransformation T(K);

            TinyVector<3> massCenter;
            T.value(FiniteElement::massCenter(), massCenter);

            const real_t fValue = F(massCenter);

            const size_t cellNumber = icell.number();

            TinyVector<FiniteElement::numberOfDegreesOfFreedom,int> index;
            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              index[l] = dofPositions(cellNumber,l);
            }

            for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
              const size_t dof
                = __degreeOfFreedomSet(testFunctionNumber,index[l]);

              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; ++m) {
                b[dof] += e(l,m)*gValues[index[m]]*fValue;
              }
            }
          }
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unknown variational operator type",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unknown problem type",
                         ErrorHandler::unexpected);
    }
    }
    ffout(2) << "- assembling second member: done";
    t.stop();
    ffout(3) << " [cost: " << t << ']';
    ffout(2) << '\n';
  }

public:
  FEMDiscretization(const Problem& p,
                    const Structured3DMesh& m,
                    const DiscretizationType& discretisationType,
                    BaseMatrix& a,
                    BaseVector& bb,
                    const DegreeOfFreedomSet& dof)
    : BaseFEMDiscretization<Structured3DMesh,
                            TypeOfDiscretization>(p, m, discretisationType, a, bb, dof)
  {
    SolverInformationCenter::instance().pushMesh(&m);
    SolverInformationCenter::instance().pushDiscretizationType(&(this->__discretizationType));
  }

  ~FEMDiscretization()
  {
    SolverInformationCenter::instance().pop();
  }
};


#endif // FEM_DISCRETIZATION_HPP

---