BoundaryConditionDiscretizationPenalty.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: BoundaryConditionDiscretizationPenalty.hpp,v 1.23 2007/08/29 08:37:14 yakoubix Exp $


#ifndef BOUNDARY_CONDITION_DISCRETIZATION_PENALTY_HPP
#define BOUNDARY_CONDITION_DISCRETIZATION_PENALTY_HPP

#include <BoundaryConditionFDMDiscretization.hpp>

#include <FiniteElementTraits.hpp>
#include <SurfaceMesh.hpp>

#include <UnAssembledMatrix.hpp>

#include <ErrorHandler.hpp>


template<typename MeshType,
         ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class BoundaryConditionDiscretizationPenalty
  : public BoundaryConditionFDMDiscretization<MeshType,
                                              TypeOfDiscretization>
{
public:
  typedef
  typename MeshType::CellType
  CellType;                     
  typedef
  typename FiniteElementTraits<CellType, TypeOfDiscretization>::Type
  FiniteElement;                
  typedef
  typename FiniteElementTraits<CellType, TypeOfDiscretization>::Transformation
  ConformTransformation;        
private:
  const real_t __epsilon;
public:
  void getDiagonal (BaseVector& Z) const
  {
    Vector<real_t>& z = dynamic_cast<Vector<real_t>&>(Z);

    // Natural Boundary conditions
    BoundaryConditionFDMDiscretization<MeshType,
                                       TypeOfDiscretization>
      ::__StandardGetDiagonalVariationalBorderBilinearOperator(*this,
                                                               z);
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *(this->__bcMeshAssociation);

    for (size_t i=0; i<this->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 SurfaceMesh& surfMesh = *(*ibcMesh).second;

        BoundaryConditionFDMDiscretization<MeshType,
                                           TypeOfDiscretization>
	  ::__meshWrapper(surfMesh,
                          typename BoundaryConditionDiscretizationPenalty
                          ::template __GetDiagonalDirichletBoundaryConditions<Vector<real_t> >(D,
                                                                                               i,
                                                                                               z,
                                                                                               *this));

      }
    }
  }


  void setMatrix (ReferenceCounting<BaseMatrix> givenA,
                  ReferenceCounting<BaseVector> b) const
  {
    switch((*givenA).type()) {
    case BaseMatrix::doubleHashedMatrix: {
      DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(*givenA);

      const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *(this->__bcMeshAssociation);

      // Variational Problem's Natural BC
      BoundaryConditionFDMDiscretization<MeshType,
                                         TypeOfDiscretization>
	::__StandardVariationalBorderBilinearOperator(*this,
                                                      A);

      // Now dealling with Dirichlet Boundary Conditions
      for (size_t i=0; i<this->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 SurfaceMesh& surfMesh = *(*ibcMesh).second;

          BoundaryConditionFDMDiscretization<MeshType,
                                             TypeOfDiscretization>
	    ::__meshWrapper(surfMesh,
                            __SetMatrixDirichletBoundaryConditions<DoubleHashedMatrix>(D, i, __epsilon,
                                                                                       A, *this));
        }
      }

      break;
    }
    case BaseMatrix::unAssembled: {
      UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(*givenA);
      A.setBoundaryConditions(this);
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected matrix type",
                         ErrorHandler::unexpected);
    }
    }
  }

  void setSecondMember (ReferenceCounting<BaseMatrix> givenA,
                        ReferenceCounting<BaseVector> givenB) const
  {
    Vector<real_t>& b = dynamic_cast<Vector<real_t>&>(*givenB);

    BoundaryConditionFDMDiscretization<MeshType,
                                       TypeOfDiscretization>::
      __StandardVariationalBorderLinearOperator(*this,
                                                b);

    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *(this->__bcMeshAssociation);

    // Dirichlet on the mesh border ...
    for (size_t i=0; i<this->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 SurfaceMesh& surfMesh = *(*ibcMesh).second;

        BoundaryConditionFDMDiscretization<MeshType,
                                           TypeOfDiscretization>
	  ::__meshWrapper(surfMesh,__SetSecondMemberDirichlet<Vector<real_t> >(D, i, b,*this));
      }
    }
  }

  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);

    // Natural Boundary Conditions.
    BoundaryConditionFDMDiscretization<MeshType,
                                       TypeOfDiscretization>::
      __StandardVariationalBorderBilinearOperatorTimesX(*this,
                                                        x, z);
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *(this->__bcMeshAssociation);

    for (size_t i=0; i<this->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 SurfaceMesh& surfMesh = *(*ibcMesh).second;

        BoundaryConditionFDMDiscretization<MeshType,
                                           TypeOfDiscretization>
	  ::__meshWrapper(surfMesh,
                          typename BoundaryConditionDiscretizationPenalty
                          ::template __TimesXDirichlet<Vector<real_t> >(D, i,
                                                                        x, z,
                                                                        *this));
      }
    }
  }

  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);

    // Natural Boundary Conditions.
    BoundaryConditionFDMDiscretization<MeshType,
                                       TypeOfDiscretization>::
      __StandardVariationalBorderBilinearOperatorTransposedTimesX(*this,
                                                                  x, z);
    const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *(this->__bcMeshAssociation);

    for (size_t i=0; i<this->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 SurfaceMesh& surfMesh = *(*ibcMesh).second;

        BoundaryConditionFDMDiscretization<MeshType,
                                           TypeOfDiscretization>
	  ::__meshWrapper(surfMesh,
                          typename BoundaryConditionDiscretizationPenalty
                          ::template __TransposedTimesXDirichlet<Vector<real_t> >(D, i,
                                                                                  x, z,
                                                                                  *this));
      }
    }
  }

  BoundaryConditionDiscretizationPenalty(const Problem& problem,
                                         const MeshType& mesh,
                                         const DegreeOfFreedomSet& dof,
                                         const real_t epsilon)
    : BoundaryConditionFDMDiscretization<MeshType,
                                         TypeOfDiscretization>(problem, mesh, dof),
      __epsilon(epsilon)
  {
    ;
  }

  ~BoundaryConditionDiscretizationPenalty()
  {
    ;
  }

private:

  template<typename MatrixType>
  class __SetMatrixDirichletBoundaryConditions
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;
    const real_t __epsilon;

    MatrixType& __A;

    const BoundaryConditionDiscretizationPenalty<MeshType,
                                                 TypeOfDiscretization>& __bc;

  public:

    template <typename SurfaceMeshType>
    void eval(const SurfaceMeshType& surfMesh) const
    {
      typedef typename SurfaceMeshType::CellType BoundaryCellType;

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

      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
        = __bc.__degreeOfFreedomSet.positionsSet(0);

      AutoPointer<ConformTransformation> pT;

      const CellType* lastCell = 0;

      for (typename SurfaceMeshType::const_iterator iface(surfMesh);
           not(iface.end()); ++iface) {
        const typename SurfaceMeshType::CellType& face = *iface;
        const BoundaryConformTransformation boundaryT(face);
        
        const CellType& cell
          = static_cast<const CellType&>(face.mother());
        
        if(lastCell != &cell) {
          pT = new ConformTransformation(cell);
          lastCell = &cell;
        }

        const size_t cellNumber = __bc.mesh().cellNumber(cell);

        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)
            * face.volume();

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

          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 indices[FiniteElement::numberOfDegreesOfFreedom];
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
            indices[l] = __bc.__degreeOfFreedomSet(__equationNumber,
                                                   dofPositions(cellNumber,l));
          }

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I = indices[l];
            for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
              const size_t J = indices[m];
              __A(I,J) += weight * WiWj(l,m) / __epsilon;
            }
          }
        }
      }
    }

    __SetMatrixDirichletBoundaryConditions(const Dirichlet& D,
                                           const size_t& equationNumber,
                                           const real_t& epsilon,
                                           MatrixType& A,
                                           const BoundaryConditionDiscretizationPenalty<MeshType,
                                                                                        TypeOfDiscretization>& bc)
      : __dirichlet(D),
        __equationNumber(equationNumber),
        __epsilon(epsilon),
        __A(A),
        __bc(bc)
    {
      ;
    }

    __SetMatrixDirichletBoundaryConditions(const __SetMatrixDirichletBoundaryConditions& S)
      : __dirichlet(S.__dirichlet),
        __equationNumber(S.__equationNumber),
        __epsilon(S.__epsilon),
        __A(S.__A),
        __bc(S.__bc)
    {
      ;
    }

    ~__SetMatrixDirichletBoundaryConditions()
    {
      ;
    }
  };

  template <typename VectorType>
  class __GetDiagonalDirichletBoundaryConditions
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;

    VectorType& __z;

    const BoundaryConditionDiscretizationPenalty<MeshType,
                                                 TypeOfDiscretization>& __bc;

  public:

    template <typename SurfaceMeshType>
    void eval(const SurfaceMeshType& surfaceMesh) const
    {
      typedef typename SurfaceMeshType::CellType BoundaryCellType;

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

      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
        = __bc.__degreeOfFreedomSet.positionsSet(0);

      AutoPointer<ConformTransformation> pT;

      const CellType* lastCell = 0;

      for (typename SurfaceMeshType::const_iterator iface(surfaceMesh);
           not(iface.end()); ++iface) {
        const typename SurfaceMeshType::CellType& face = *iface;
        const BoundaryConformTransformation boundaryT(face);

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

        if(lastCell != &cell) {
          pT = new ConformTransformation(cell);
          lastCell = &cell;
        }

        const size_t cellNumber = __bc.mesh().cellNumber(cell);

        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)
            * face.volume();

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

          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] = __bc.__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 / __bc.__epsilon;
          }
        }
      }
    }


    __GetDiagonalDirichletBoundaryConditions(const Dirichlet& D,
                                             const size_t equationNumber,
                                             VectorType& z,
                                             const BoundaryConditionDiscretizationPenalty& bc)
      : __dirichlet(D),
        __equationNumber(equationNumber),
        __z(z),
        __bc(bc)      
    {
      ;
    }

    __GetDiagonalDirichletBoundaryConditions(const __GetDiagonalDirichletBoundaryConditions& G)
      : __dirichlet(G.__dirichlet),
        __equationNumber(G.__equationNumber),
        __z(G.__z),
        __bc(G.__bc)
    {
      ;
    }

    ~__GetDiagonalDirichletBoundaryConditions()
    {
      ;
    }
  };

  template <typename VectorType>
  class __TimesXDirichlet
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;

    const VectorType& __x;

    VectorType& __z;

    const BoundaryConditionDiscretizationPenalty<MeshType,
                                                 TypeOfDiscretization>& __bc;  

  public:

    template <typename SurfaceMeshType>
    void eval(const SurfaceMeshType& surfaceMesh) const
    {
      typedef typename SurfaceMeshType::CellType BoundaryCellType;

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

      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
        = __bc.__degreeOfFreedomSet.positionsSet(0);

      AutoPointer<ConformTransformation> pT;

      const CellType* lastCell = 0;

      for (typename SurfaceMeshType::const_iterator iface(surfaceMesh);
           not(iface.end()); ++iface) {
        const typename SurfaceMeshType::CellType& face = *iface;
        const BoundaryConformTransformation boundaryT(face);

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

        if(lastCell != &cell) {
          pT = new ConformTransformation(cell);
          lastCell = &cell;
        }

        const size_t cellNumber = __bc.mesh().cellNumber(cell);

        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)
            * face.volume();

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

          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 indices[FiniteElement::numberOfDegreesOfFreedom];
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
            indices[l] = __bc.__degreeOfFreedomSet(__equationNumber,
                                                   dofPositions(cellNumber, l));
          }
            
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I = indices[l];
            if (not(__bc.__dirichletList[I])) {
              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
                const size_t J = indices[m];
                if (not(__bc.__dirichletList[J])) {
                  __z[J] += WiWj(l,m) * weight * __x[I] / __bc.__epsilon;
                }
              }
            }
          }
        }
      }
    }

    __TimesXDirichlet(const Dirichlet& dirichlet,
                      const size_t equationNumber,
                      const VectorType& x,
                      VectorType& z,
                      const BoundaryConditionDiscretizationPenalty<MeshType,
                                                                   TypeOfDiscretization>& bc)
      : __dirichlet(dirichlet),
        __equationNumber(equationNumber),
        __x(x),
        __z(z),
        __bc(bc)
    {
      ;
    }

    __TimesXDirichlet(const __TimesXDirichlet& T)
      : __dirichlet(T.__dirichlet),
        __equationNumber(T.__equationNumber),
        __x(T.__x),
        __z(T.__z),
        __bc(T.__bc)
    {
      ;
    }

    ~__TimesXDirichlet()
    {
      ;
    }
  };

  template <typename VectorType>
  class __TransposedTimesXDirichlet
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;

    const VectorType& __x;

    VectorType& __z;

    const BoundaryConditionDiscretizationPenalty<MeshType,
                                                 TypeOfDiscretization>& __bc;  

  public:

    template <typename SurfaceMeshType>
    void eval(const SurfaceMeshType& surfaceMesh) const
    {
      typedef typename SurfaceMeshType::CellType BoundaryCellType;

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

      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
        = __bc.__degreeOfFreedomSet.positionsSet(0);

      AutoPointer<ConformTransformation> pT;

      const CellType* lastCell = 0;

      for (typename SurfaceMeshType::const_iterator iface(surfaceMesh);
           not(iface.end()); ++iface) {
        const typename SurfaceMeshType::CellType& face = *iface;
        const BoundaryConformTransformation boundaryT(face);

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

        if(lastCell != &cell) {
          pT = new ConformTransformation(cell);
          lastCell = &cell;
        }

        const size_t cellNumber = __bc.mesh().cellNumber(cell);

        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)
            * face.volume();

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

          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 indices[FiniteElement::numberOfDegreesOfFreedom];
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; ++l) {
            indices[l] = __bc.__degreeOfFreedomSet(__equationNumber,
                                                   dofPositions(cellNumber, l));
          }
            
          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I = indices[l];
            if (not(__bc.__dirichletList[I])) {
              for (size_t m=0; m<FiniteElement::numberOfDegreesOfFreedom; m++) {
                const size_t J = indices[m];
                if (not(__bc.__dirichletList[J])) {
                  __z[I] += WiWj(l,m) * weight * __x[J] / __bc.__epsilon;
                }
              }
            }
          }
        }
      }
    }

    __TransposedTimesXDirichlet(const Dirichlet& dirichlet,
                                const size_t equationNumber,
                                const VectorType& x,
                                VectorType& z,
                                const BoundaryConditionDiscretizationPenalty<MeshType,
                                                                             TypeOfDiscretization>& bc)
      : __dirichlet(dirichlet),
        __equationNumber(equationNumber),
        __x(x),
        __z(z),
        __bc(bc)
    {
      ;
    }

    __TransposedTimesXDirichlet(const __TransposedTimesXDirichlet& T)
      : __dirichlet(T.__dirichlet),
        __equationNumber(T.__equationNumber),
        __x(T.__x),
        __z(T.__z),
        __bc(T.__bc)
    {
      ;
    }

    ~__TransposedTimesXDirichlet()
    {
      ;
    }
  };

  template <typename VectorType>
  class __SetSecondMemberDirichlet
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;

    VectorType& __b;

    const BoundaryConditionDiscretizationPenalty<MeshType,
                                                 TypeOfDiscretization>& __bc;
  
  public:

    template <typename SurfaceMeshType>
    void eval(const SurfaceMeshType& surfaceMesh) const
    {
      typedef typename SurfaceMeshType::CellType BoundaryCellType;

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

      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
        = __bc.__degreeOfFreedomSet.positionsSet(0);

      AutoPointer<ConformTransformation> pT;

      const CellType* lastCell = 0;

      for (typename SurfaceMeshType::const_iterator iface(surfaceMesh);
           not(iface.end()); ++iface) {
        const typename SurfaceMeshType::CellType& face = *iface;
        const BoundaryConformTransformation boundaryT(face);

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

        if(lastCell != &cell) {
          pT = new ConformTransformation(cell);
          lastCell = &cell;
        }

        const ConformTransformation T = *pT;

        const size_t cellNumber = __bc.mesh().cellNumber(cell);

        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)
            * face.volume();

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

          const real_t GValue = __dirichlet.g(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] = dofPositions(cellNumber, l);
          }

          for (size_t l=0; l<FiniteElement::numberOfDegreesOfFreedom; l++) {
            const size_t I = __bc.__degreeOfFreedomSet(__equationNumber,
                                                       indices[l]);
            // Assembling second member
            __b[I] += W[l] * weight * GValue / __bc.__epsilon;
          }
        }
      }
    }

    __SetSecondMemberDirichlet(const Dirichlet& dirichlet,
                               const size_t equationNumber,
                               VectorType& b,
                               const BoundaryConditionDiscretizationPenalty<MeshType,
                                                                            TypeOfDiscretization>& bc)
      : __dirichlet(dirichlet),
        __equationNumber(equationNumber),
        __b(b),
        __bc(bc)
    {
      ;
    }

    __SetSecondMemberDirichlet(const __SetSecondMemberDirichlet& S)
      : __dirichlet(S.__dirichlet),
        __equationNumber(S.__equationNumber),
        __b(S.__b),
        __bc(S.__bc)
    {
      ;
    }

    ~__SetSecondMemberDirichlet()
    {
      ;
    }
  };

};

#endif // BOUNDARY_CONDITION_DISCRETIZATION_PENALTY_HPP

---