BoundaryConditionCommonFEMDiscretization.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: BoundaryConditionCommonFEMDiscretization.hpp,v 1.30 2007/11/02 20:36:13 delpinux Exp $


#ifndef BOUNDARY_CONDITION_COMMON_FEM_DISCRETIZATION_HPP
#define BOUNDARY_CONDITION_COMMON_FEM_DISCRETIZATION_HPP

#include <FiniteElementTraits.hpp>

#include <BoundaryConditionDiscretization.hpp>

#include <SurfaceMeshGenerator.hpp>

#include <SurfaceMeshOfQuadrangles.hpp>
#include <SurfaceMeshOfTriangles.hpp>

#include <MatrixManagement.hpp>

#include <AutoPointer.hpp>
#include <ConformTransformation.hpp>

#include <BoundaryMeshAssociation.hpp>
#include <BoundaryConditionSurfaceMeshAssociation.hpp>

#include <Stringify.hpp>
#include <ErrorHandler.hpp>

template <typename MeshType,
          ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class BoundaryConditionCommonFEMDiscretization
  : public BoundaryConditionDiscretization
{
public:
  typedef
  typename MeshType::CellType
  CellType;                     
  typedef
  typename FiniteElementTraits<CellType, TypeOfDiscretization>::Type
  FiniteElement;                
  typedef
  typename FiniteElement::QuadratureType
  QuadratureType;               
  typedef
  typename FiniteElementTraits<CellType, TypeOfDiscretization>::Transformation
  ConformTransformation;        
protected:
  const MeshType& __mesh;       
  ReferenceCounting<BoundaryConditionSurfaceMeshAssociation>
  __bcMeshAssociation;          
  void __checkBoundaryMeshAssociation(const BoundaryMeshAssociation& bma)
  {
    for (BoundaryMeshAssociation::const_iterator i = bma.begin();
         i != bma.end(); ++i) {
      if (i->second == 0) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "The mesh of boundary "
                           +stringify(*(i->first))+" was not generated",
                           ErrorHandler::unexpected);
      } else {
        const SurfaceMesh& s = *(i->second);
        if (not(s.isAssociatedTo(__mesh))) {
          throw ErrorHandler(__FILE__,__LINE__,
                             "The mesh of boundary "
                             +stringify(*(i->first))
                             +" is not related to the volume mesh"
                             +"(check your mesh and/or report the error)",
                             ErrorHandler::normal);
        }
      }
    }
  }

  void __associatesDefinedMeshToBoundaryConditions(const BoundaryMeshAssociation& bma)
  {
    __bcMeshAssociation = new BoundaryConditionSurfaceMeshAssociation(__problem,bma);
  }

public:

  virtual void associatesMeshesToBoundaryConditions()
  {
    BoundaryMeshAssociation bma(__problem, __mesh);
    this->__checkBoundaryMeshAssociation(bma);
    this->__associatesDefinedMeshToBoundaryConditions(bma);
  }

  const MeshType& mesh() const
  {
    return __mesh;
  }

  MeshType& mesh()
  {
    return __mesh;
  }

  BoundaryConditionCommonFEMDiscretization(const Problem& problem,
                                           const MeshType& aMesh,
                                           const DegreeOfFreedomSet& dof)
    : BoundaryConditionDiscretization(problem,dof),
      __mesh(aMesh)
  {
    ;
  }

  virtual ~BoundaryConditionCommonFEMDiscretization()
  {
    ;
  }

protected:

  template <typename MatrixType,
            typename BoundaryMeshType>
  void
  __setVariationalBoundaryConditionsAlphaUV(const ScalarFunctionBase& alpha,
                                            const size_t& equationNumber,
                                            const size_t& unknownNumber,
                                            MatrixType& A,
                                            const BoundaryMeshType& surfaceMesh) const
  {
    typedef typename BoundaryMeshType::CellType BoundaryCellType;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::Transformation
      BoundaryConformTransformation;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::JacobianTransformation
      BoundaryConformTransformationJacobian;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::Type
      BoundaryFiniteElement;

    typedef
      typename BoundaryFiniteElement::QuadratureType
      BoundaryQuadratureType;

    const FiniteElement& finiteElement
      = FiniteElement::instance();

    const BoundaryQuadratureType& referenceBoundaryQuadrature
      = BoundaryQuadratureType::instance();

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    AutoPointer<ConformTransformation> pT;

    const CellType* currentCell = 0;

    for (typename BoundaryMeshType::const_iterator icell(surfaceMesh);
         not(icell.end()); ++icell) {
      const BoundaryCellType& cell = *icell;
      const BoundaryConformTransformation boundaryT(cell);
      const BoundaryConformTransformationJacobian boundaryJ(boundaryT);

      const CellType& K
        = static_cast<const CellType&>(cell.mother());

      if(currentCell != &K) {
        pT = new ConformTransformation(K);
        currentCell = &K;
      }
      const size_t cellNumber = __mesh.cellNumber(K);

      const ConformTransformation& T = *pT;

      for (size_t k=0; k<BoundaryQuadratureType::numberOfQuadraturePoints;
           k++) {
        // computes local quadrature vertex
        TinyVector<3, real_t> q;
        boundaryT.value(referenceBoundaryQuadrature[k], q);

        const real_t weight
          = referenceBoundaryQuadrature.weight(k)
          * boundaryJ.jacobianDet();

        TinyVector<3> coordinates;
        T.invertT(q, coordinates);

        const real_t AlphaValue = alpha(q);

        typename FiniteElement::ElementaryVector W;
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++)
          W[l] = finiteElement.W(l,coordinates);

        typename FiniteElement::ElementaryMatrix WiWj;

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++)
          for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++)
            WiWj(l,m) = W[l] * W[m];

        size_t indicesI[FiniteElement::numberOfDegreesOfFreedom];
        size_t indicesJ[FiniteElement::numberOfDegreesOfFreedom];
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          const size_t variableDOFNumber = dofPositions(cellNumber,l);
          indicesI[l] = __degreeOfFreedomSet(unknownNumber,
                                             variableDOFNumber);
          indicesJ[l] = __degreeOfFreedomSet(equationNumber,
                                             variableDOFNumber);
        }

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I = indicesI[l];
          for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
            const size_t J = indicesJ[m];
            A(I,J) += weight * AlphaValue * WiWj(l,m);
          }
        }
      }
    }
  }  

  template <typename BoundaryMeshType>
  void __setSecondMemberDirichlet(const Dirichlet& dirichlet,
                                  const size_t& equationNumber,
                                  const BoundaryMeshType& surfMesh) const
  {
    typedef FiniteElementTraits<typename BoundaryMeshType::CellType,
                                TypeOfDiscretization>
      FEType;

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);

    typedef typename FEType::Transformation ConformTransformation;


    for (typename BoundaryMeshType::const_iterator iBorderCell(surfMesh);
         not(iBorderCell.end()); ++iBorderCell) {
      const typename BoundaryMeshType::CellType& borderCell = *iBorderCell;

      const size_t cellNumber
        = __mesh.cellNumber(static_cast<const CellType&>(borderCell.mother()));
      const size_t& faceNumber = borderCell.motherCellFaceNumber();

      for (size_t i=0; i<FiniteElement::numberOfFaceLivingDegreesOfFreedom; ++i) {
        const size_t& dofNumber = dofPositions(cellNumber,
                                               FiniteElement::facesDOF[faceNumber][i]);
        const size_t I = __degreeOfFreedomSet(equationNumber,
                                              dofPositions(cellNumber,
                                                           FiniteElement::facesDOF[faceNumber][i]));
        if (not(__dirichletList[I])) {
          const real_t GValue = dirichlet.g(dofPositions.vertex(dofNumber));
          __dirichletValues[I] = GValue;
          __dirichletList[I] = true;
        }
      }
    }
  }

  template <typename VectorType,
            typename BoundaryMeshType>
  void
  __getDiagonalNaturalBoundaryConditions(const ScalarFunctionBase& alpha,
                                         const size_t& equationNumber,
                                         VectorType& z,
                                         const BoundaryMeshType& surfaceMesh) const
  {
    typedef typename BoundaryMeshType::CellType BoundaryCellType;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::Transformation
      BoundaryConformTransformation;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::JacobianTransformation
      BoundaryConformTransformationJacobian;

    typedef
      typename FiniteElementTraits<BoundaryCellType,TypeOfDiscretization>::Type
      BoundaryFiniteElement;

    typedef typename BoundaryFiniteElement::QuadratureType BoundaryQuadratureType;

    const FiniteElement& finiteElement
      = FiniteElement::instance();

    const BoundaryQuadratureType& referenceBoundaryQuadrature
      = BoundaryQuadratureType::instance();
  
    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    AutoPointer<ConformTransformation> pT;

    const CellType* currentCell = 0;

    for (typename BoundaryMeshType::const_iterator icell(surfaceMesh);
         not(icell.end()); ++icell) {
      const BoundaryCellType& cell = *icell;
      const BoundaryConformTransformation boundaryT(cell);
      const BoundaryConformTransformationJacobian boundaryJ(boundaryT);

      const CellType& K
        = static_cast<const CellType&>(cell.mother());

      if(currentCell != &K) {
        pT = new ConformTransformation(K);
        currentCell = &K;
      }
      const size_t cellNumber = __mesh.cellNumber(K);
      const ConformTransformation& T = *pT;

      for (size_t k=0; k<BoundaryQuadratureType::numberOfQuadraturePoints;
           k++) {
        // computes local quadrature vertex
        TinyVector<3, real_t> q;
        boundaryT.value(referenceBoundaryQuadrature[k], q);

        const real_t weight
          = referenceBoundaryQuadrature.weight(k)
          * boundaryJ.jacobianDet();

        TinyVector<3, real_t> coordinates;
        T.invertT(q, coordinates);

        const real_t AlphaValue = alpha(q);

        typename FiniteElement::ElementaryVector W;

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) 
          W[l] = finiteElement.W(l,coordinates);

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

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t& I = indices[l];
          z[I] += W[l] * W[l] * weight * AlphaValue;
        }
      }
    }
  }

  template <typename VectorType,
            typename BoundaryMeshType>
  void
  __setVariationalBoundaryConditionsFV(const ScalarFunctionBase& g,
                                       const size_t& equationNumber,
                                       VectorType& b,
                                       const BoundaryMeshType& boundaryMesh) const
  {
    typedef typename BoundaryMeshType::CellType BoundaryCellType;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::Transformation
      BoundaryConformTransformation;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::JacobianTransformation
      BoundaryConformTransformationJacobian;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
      TypeOfDiscretization>::Type
      BoundaryFiniteElement;

    typedef
      typename BoundaryFiniteElement::QuadratureType
      BoundaryQuadratureType;

    const FiniteElement& finiteElement
      = FiniteElement::instance();

    const BoundaryQuadratureType& referenceBoundaryQuadrature
      = BoundaryQuadratureType::instance();

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    AutoPointer<ConformTransformation> pT;

    const CellType* currentCell = 0;

    for (typename BoundaryMeshType::const_iterator icell(boundaryMesh);
         not(icell.end()); ++icell) {
      const BoundaryCellType& cell = *icell;
      const BoundaryConformTransformation boundaryT(cell);
      const BoundaryConformTransformationJacobian boundaryJ(boundaryT);

      const CellType& K
        = static_cast<const CellType&>(cell.mother());

      if (currentCell != &K) {
        pT = new ConformTransformation(K);
        currentCell = &K;
      }

      const ConformTransformation& T = *pT;

      for (size_t k=0; k<BoundaryQuadratureType::numberOfQuadraturePoints;
           k++) {
        // computes local quadrature vertex
        TinyVector<3, real_t> q;
        boundaryT.value(referenceBoundaryQuadrature[k], q);

        const real_t weight
          = referenceBoundaryQuadrature.weight(k)
          * boundaryJ.jacobianDet();

        TinyVector<3> coordinates;
        T.invertT(q, coordinates);

        const real_t GValue = g(q);

        typename FiniteElement::ElementaryVector W;
        size_t indices[FiniteElement::numberOfDegreesOfFreedom];

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          W[l] = finiteElement.W(l,coordinates);
          indices[l] = dofPositions(__mesh.cellNumber(K),l);
        }

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I = __degreeOfFreedomSet(equationNumber,
                                                indices[l]);
          b[I] += weight * GValue * W[l];
        }
      }
    }
  }

  template <typename VectorType,
            typename BoundaryMeshType>
  void
  __variationalBoundaryConditionsAlphaUVTimesX(const ScalarFunctionBase& alpha,
                                               const size_t& equationNumber,
                                               const size_t& unknownNumber,
                                               const VectorType& x,
                                               VectorType& z,
                                               const BoundaryMeshType& boundaryMesh) const
  {
    typedef typename BoundaryMeshType::CellType BoundaryCellType;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                   TypeOfDiscretization>::Transformation
      BoundaryConformTransformation;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                  TypeOfDiscretization>::JacobianTransformation
      BoundaryConformTransformationJacobian;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                   TypeOfDiscretization>::Type
      BoundaryFiniteElement;

    typedef
      typename BoundaryFiniteElement::QuadratureType
      BoundaryQuadratureType;

    const FiniteElement& finiteElement
      = FiniteElement::instance();

    const BoundaryQuadratureType& referenceBoundaryQuadrature
      = BoundaryQuadratureType::instance();

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    AutoPointer<ConformTransformation> pT;

    const CellType* currentCell = 0;

    for (typename BoundaryMeshType::const_iterator icell(boundaryMesh);
         not(icell.end()); ++icell) {
      const BoundaryCellType& cell = *icell;
      const BoundaryConformTransformation boundaryT(cell);
      const BoundaryConformTransformationJacobian boundaryJ(boundaryT);

      const CellType& K = static_cast<const CellType&>(cell.mother());

      if(currentCell != &K) {
        pT = new ConformTransformation(K);
        currentCell = &K;
      }

      const size_t cellNumber = __mesh.cellNumber(K);

      const ConformTransformation& T = *pT;

      for (size_t k=0; k<BoundaryQuadratureType::numberOfQuadraturePoints;
           k++) {
        // computes local quadrature vertex
        TinyVector<3, real_t> q;
        boundaryT.value(referenceBoundaryQuadrature[k], q);

        const real_t weight
          = referenceBoundaryQuadrature.weight(k)
          * boundaryJ.jacobianDet();

        TinyVector<3> coordinates;
        T.invertT(q, coordinates);

        const real_t AlphaValue = alpha(q);

        typename FiniteElement::ElementaryVector W;
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) 
          W[l] = finiteElement.W(l,coordinates);

        typename FiniteElement::ElementaryMatrix WiWj;

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++)
          for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++)
            WiWj(l,m) = W[l] * W[m];

        size_t indicesI[FiniteElement::numberOfDegreesOfFreedom];
        size_t indicesJ[FiniteElement::numberOfDegreesOfFreedom];
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          const size_t variableDOFNumber = dofPositions(cellNumber,l);
          indicesI[l] = __degreeOfFreedomSet(unknownNumber,
                                             variableDOFNumber);
          indicesJ[l] = __degreeOfFreedomSet(equationNumber,
                                             variableDOFNumber);
        }
            
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I = indicesI[l];
          if (not(__dirichletList[I])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t J = indicesJ[m];
              if (not(__dirichletList[J])) {
                z[I] += WiWj(l,m) * weight * AlphaValue * x[J];
              }
            }
          }
        }
      }
    }
  }

  template <typename VectorType,
            typename BoundaryMeshType>
  void
  __variationalBoundaryConditionsAlphaUVTransposedTimesX(const ScalarFunctionBase& alpha,
                                                         const size_t& equationNumber,
                                                         const size_t& unknownNumber,
                                                         const VectorType& x,
                                                         VectorType& z,
                                                         const BoundaryMeshType& boundaryMesh) const
  {
    typedef typename BoundaryMeshType::CellType BoundaryCellType;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                   TypeOfDiscretization>::Transformation
      BoundaryConformTransformation;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                   TypeOfDiscretization>::JacobianTransformation
      BoundaryConformTransformationJacobian;

    typedef
      typename FiniteElementTraits<BoundaryCellType,
                                   TypeOfDiscretization>::Type
      BoundaryFiniteElement;

    typedef
      typename BoundaryFiniteElement::QuadratureType
      BoundaryQuadratureType;

    const FiniteElement& finiteElement
      = FiniteElement::instance();

    const BoundaryQuadratureType& referenceBoundaryQuadrature
      = BoundaryQuadratureType::instance();

    const ScalarDegreeOfFreedomPositionsSet& dofPositions
      = this->__degreeOfFreedomSet.positionsSet(0);
  
    AutoPointer<ConformTransformation> pT;

    const CellType* currentCell = 0;

    for (typename BoundaryMeshType::const_iterator icell(boundaryMesh);
         not(icell.end()); ++icell) {
      const BoundaryCellType& cell = *icell;
      const BoundaryConformTransformation boundaryT(cell);
      const BoundaryConformTransformationJacobian boundaryJ(boundaryT);

      const CellType& K = static_cast<const CellType&>(cell.mother());

      if(currentCell != &K) {
        pT = new ConformTransformation(K);
        currentCell = &K;
      }

      const size_t cellNumber = __mesh.cellNumber(K);

      const ConformTransformation& T = *pT;

      for (size_t k=0; k<BoundaryQuadratureType::numberOfQuadraturePoints;
           k++) {
        // computes local quadrature vertex
        TinyVector<3, real_t> q;
        boundaryT.value(referenceBoundaryQuadrature[k], q);

        const real_t weight
          = referenceBoundaryQuadrature.weight(k)
          * boundaryJ.jacobianDet();

        TinyVector<3> coordinates;
        T.invertT(q, coordinates);

        const real_t AlphaValue = alpha(q);

        typename FiniteElement::ElementaryVector W;
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) 
          W[l] = finiteElement.W(l,coordinates);

        typename FiniteElement::ElementaryMatrix WiWj;

        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++)
          for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++)
            WiWj(l,m) = W[l] * W[m];

        size_t indicesI[FiniteElement::numberOfDegreesOfFreedom];
        size_t indicesJ[FiniteElement::numberOfDegreesOfFreedom];
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
          const size_t variableDOFNumber = dofPositions(cellNumber,l);
          indicesI[l] = __degreeOfFreedomSet(unknownNumber,
                                             variableDOFNumber);
          indicesJ[l] = __degreeOfFreedomSet(equationNumber,
                                             variableDOFNumber);
        }
            
        for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
          const size_t I = indicesI[l];
          if (not(__dirichletList[I])) {
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t J = indicesJ[m];
              if (not(__dirichletList[J])) {
                z[J] += WiWj(l,m) * weight * AlphaValue * x[I];
              }
            }
          }
        }
      }
    }
  }

  template <typename FunctionClass>
  static void __meshWrapper(const SurfaceMesh& surfmesh, const FunctionClass& F)
  {
    switch (surfmesh.type()) {
    case Mesh::surfaceMeshTriangles: {
      const SurfaceMeshOfTriangles& surfTria
        = dynamic_cast<const SurfaceMeshOfTriangles&>(surfmesh);
      F.eval(surfTria);
      break;
    }
    case Mesh::surfaceMeshQuadrangles: {
      const SurfaceMeshOfQuadrangles& surfQuad
        = dynamic_cast<const SurfaceMeshOfQuadrangles&>(surfmesh);
      F.eval(surfQuad);
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unknown surface mesh type",
                         ErrorHandler::unexpected);
    }
    }
  }



  template <typename VectorType>
  static void
  __StandardGetDiagonalVariationalBorderBilinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                         TypeOfDiscretization>& __bc,
                                                         VectorType& z)
  {
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bc.__bcMeshAssociation;

    for (BoundaryConditionSurfaceMeshAssociation
           ::BilinearBorderOperatorMeshAssociation
           ::const_iterator i = bcMeshAssociation.beginOfBilinear();
         i != bcMeshAssociation.endOfBilinear(); ++i) {

      switch (i->first->type()) {
      case VariationalBilinearBorderOperator::alphaUV: {
        const VariationalBorderOperatorAlphaUV& V
          = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);
        const Mesh& boundaryMesh = (*i->second);
        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          __bc.__getDiagonalNaturalBoundaryConditions(V.alpha(),
                                                      V.testFunctionNumber(), z,
                                                      dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          __bc.__getDiagonalNaturalBoundaryConditions(V.alpha(),
                                                      V.testFunctionNumber(), z,
                                                      dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardGetDiagonalDirichletBorderBilinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                                      TypeOfDiscretization>& bc,
                                                       VectorType& z)
  {
    for (size_t i=0; i<bc.__dirichletList.size(); ++i) {
      if (bc.__dirichletList[i]) {
        z[i] = 1;
      }
    }
  }

  template <typename MatrixType>
  static void
  __StandardVariationalBorderBilinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                             TypeOfDiscretization>& __bc,
                                              MatrixType& A)
  {
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bc.__bcMeshAssociation;

    for (BoundaryConditionSurfaceMeshAssociation
           ::BilinearBorderOperatorMeshAssociation
           ::const_iterator i = bcMeshAssociation.beginOfBilinear();
         i != bcMeshAssociation.endOfBilinear(); ++i) {

      switch (i->first->type()) {
      case VariationalBilinearBorderOperator::alphaUV: {
        const VariationalBorderOperatorAlphaUV& V
          = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);
        
        const Mesh& boundaryMesh = *i->second;

        switch (boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          __bc.__setVariationalBoundaryConditionsAlphaUV<MatrixType> (V.alpha(),
                                                                      V.testFunctionNumber(),
                                                                      V.unknownNumber(),
                                                                      A,
                                                                      dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          __bc.__setVariationalBoundaryConditionsAlphaUV<MatrixType> (V.alpha(),
                                                                      V.testFunctionNumber(),
                                                                      V.unknownNumber(),
                                                                      A,
                                                                      dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
    }
  }

  template <typename MatrixType>
  static void
  __StandardDirichletBorderBilinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                           TypeOfDiscretization>& __bc,
                                            MatrixType& A)

  {
    for (size_t i=0; i<__bc.__dirichletList.size(); ++i) {
      if (__bc.__dirichletList[i]) {
        for (DoubleHashedMatrix::iterator lineBrowser = A.beginOfLine(i);
             lineBrowser != A.endOfLine(i); ++lineBrowser) {
#warning should use a syntax like "M.line(i) = 0;"
          lineBrowser.second() = 0;
        }
        // makes the column i to be zero. Here the elimination modifies the second member.
        for (DoubleHashedMatrix::iterator columnBrowser = A.beginOfColumn(i);
             columnBrowser != A.endOfColumn(i); ++columnBrowser) {
          columnBrowser.second() = 0;
        }
        A(i,i) = 1.; // sets the diagonal to 1
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardVariationalBorderLinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                           TypeOfDiscretization>& __bc,
                                            VectorType& b)
  {
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bc.__bcMeshAssociation;

    for (BoundaryConditionSurfaceMeshAssociation
           ::LinearBorderOperatorMeshAssociation
           ::const_iterator i = bcMeshAssociation.beginOfLinear();
         i != bcMeshAssociation.endOfLinear(); ++i) {

      switch (i->first->type()) {
      case VariationalLinearBorderOperator::FV: {
        const VariationalBorderOperatorFV& V
          = dynamic_cast<const VariationalBorderOperatorFV&>(*i->first);

        const Mesh& boundaryMesh = (*i->second);
        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          __bc.__setVariationalBoundaryConditionsFV(V.f(), V.testFunctionNumber(), b,
                                                    dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          __bc.__setVariationalBoundaryConditionsFV(V.f(), V.testFunctionNumber(), b,
                                                    dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
    }
  }

  template <typename MatrixType,
            typename VectorType>
  static void
  __StandardDirichletBorderLinearOperator(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                         TypeOfDiscretization>& bc,
                                          MatrixType& A,
                                          VectorType& b)
  {
    bc.__dirichletValues.resize(bc.__degreeOfFreedomSet.size());
    bc.__dirichletValues = 0;
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *bc.__bcMeshAssociation;

    // Dirichlet on the mesh border ...
    for (size_t i=0; i<bc.problem().numberOfUnknown(); ++i) {
      for (BoundaryConditionSurfaceMeshAssociation
             ::DirichletMeshAssociation
             ::const_iterator ibcMesh = bcMeshAssociation.bc(i).begin();
           ibcMesh != bcMeshAssociation.bc(i).end();
           ++ibcMesh) {
        const Dirichlet& D = *ibcMesh->first;
        const Mesh& boundaryMesh = *ibcMesh->second;
        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          bc.__setSecondMemberDirichlet(D,i,
                                        dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          bc.__setSecondMemberDirichlet(D,i,
                                        dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
      }
    }

    Vector<bool> temp(bc.__dirichletList);
    bc.__dirichletList=false;
    Vector<real_t> y (bc.__dirichletValues.size());
    A.timesX(bc.__dirichletValues,y);
    b -= y;
    bc.__dirichletList=temp;

    for (size_t i=0; i<b.size(); ++i) {
      if (bc.__dirichletList[i]) {
        b[i] = bc.__dirichletValues[i];
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardVariationalBorderBilinearOperatorTimesX(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                                   TypeOfDiscretization>& __bc,
                                                    const VectorType& x,
                                                    VectorType& z)
  {
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bc.__bcMeshAssociation;
    for (BoundaryConditionSurfaceMeshAssociation
           ::BilinearBorderOperatorMeshAssociation
           ::const_iterator i = bcMeshAssociation.beginOfBilinear();
         i != bcMeshAssociation.endOfBilinear(); ++i) {

      switch (i->first->type()) {
      case VariationalBilinearBorderOperator::alphaUV: {

        const VariationalBorderOperatorAlphaUV& V
          = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);

        const Mesh& boundaryMesh = (*i->second);

        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles : {
          __bc.__variationalBoundaryConditionsAlphaUVTimesX(V.alpha(),
                                                            V.testFunctionNumber(),
                                                            V.unknownNumber(),
                                                            x,
                                                            z,
                                                            dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles : {
          __bc.__variationalBoundaryConditionsAlphaUVTimesX(V.alpha(),
                                                            V.testFunctionNumber(),
                                                            V.unknownNumber(),
                                                            x,
                                                            z,
                                                            dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardVariationalBorderBilinearOperatorTransposedTimesX(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                                             TypeOfDiscretization>& __bc,
                                                              const VectorType& x,
                                                              VectorType& z)
  {
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bc.__bcMeshAssociation;
    for (BoundaryConditionSurfaceMeshAssociation
           ::BilinearBorderOperatorMeshAssociation
           ::const_iterator i = bcMeshAssociation.beginOfBilinear();
         i != bcMeshAssociation.endOfBilinear(); ++i) {

      switch (i->first->type()) {
      case VariationalBilinearBorderOperator::alphaUV: {

        const VariationalBorderOperatorAlphaUV& V
          = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);

        const Mesh& boundaryMesh = *i->second;

        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          __bc.__variationalBoundaryConditionsAlphaUVTransposedTimesX(V.alpha(),
                                                                      V.testFunctionNumber(),
                                                                      V.unknownNumber(),
                                                                      x, z,
                                                                      dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          __bc.__variationalBoundaryConditionsAlphaUVTransposedTimesX(V.alpha(),
                                                                      V.testFunctionNumber(),
                                                                      V.unknownNumber(),
                                                                      x, z,
                                                                      dynamic_cast<const SurfaceMeshOfQuadrangles&>(boundaryMesh));
          break;
        }
        default: {        
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardDirichletBorderBilinearOperatorTimesX(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                                 TypeOfDiscretization>& bc,
                                                  const VectorType& x,
                                                  VectorType& z)
  {
    for (size_t i=0; i<bc.__dirichletList.size(); ++i) {
      if (bc.__dirichletList[i]) {
        z[i] = x[i];
      }
    }
  }

  template <typename VectorType>
  static void
  __StandardDirichletBorderBilinearOperatorTransposedTimesX(const BoundaryConditionCommonFEMDiscretization<MeshType,
                                                                                                           TypeOfDiscretization>& bc,
                                                            const VectorType& x,
                                                            VectorType& z)
  {
    for (size_t i=0; i<bc.__dirichletList.size(); ++i) {
      if (bc.__dirichletList[i]) {
        z[i] = x[i];
      }
    }
  }
};

#endif // BOUNDARY_CONDITION_COMMON_FEM_DISCRETIZATION_HPP

---