BoundaryConditionDiscretizationElimination.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: BoundaryConditionDiscretizationElimination.hpp,v 1.20 2007/08/18 15:57:42 delpinux Exp $

#ifndef BOUNDARY_CONDITION_DISCRETIZATION_ELIMINATION_HPP
#define BOUNDARY_CONDITION_DISCRETIZATION_ELIMINATION_HPP

#include <BoundaryConditionFDMDiscretization.hpp>
#include <UnAssembledMatrix.hpp>

template<typename MeshType,
         ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
class BoundaryConditionDiscretizationElimination
  : public BoundaryConditionFDMDiscretization<MeshType,
                                              TypeOfDiscretization>
{
public:
  typedef typename MeshType::CellType CellType;

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

private:
  class __SetSecondMemberDirichlet
  {
  private:
    const Dirichlet& __dirichlet;

    const size_t __equationNumber;

    const BoundaryConditionDiscretizationElimination<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 BoundaryQuadratureType& referenceBoundaryQuadrature
        = BoundaryQuadratureType::instance();

      fferr(1) << __FILE__ << ":" << __LINE__
               << "\n########################################################\n"
               << "Implementation not finished! use of dofPosition required\n"
               << "########################################################\n";

      for (typename SurfaceMeshType::const_iterator icell(surfMesh);
           not(icell.end()); ++icell) {
        const BoundaryCellType& cell = *icell;

        const BoundaryConformTransformation boundaryT(cell);

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

          // Index of the vertex of the mesh the closer to q.
#warning Use of Index is not necessary
          const Index& iv = __bc.mesh().vertexIndex(q);
          const Vertex& V = __bc.mesh().vertex(iv);

          const real_t GValue = __dirichlet.g(V);

          const size_t I = __bc.__degreeOfFreedomSet(__equationNumber,
                                                     __bc.mesh().vertexNumber(V));
          __bc.__dirichletValues[I] = GValue;
          __bc.__dirichletList[I] = true;
        }
      }
    }

    __SetSecondMemberDirichlet(const Dirichlet &D ,
                               const size_t equationNumber,
                               const BoundaryConditionDiscretizationElimination<MeshType,
                                                                                TypeOfDiscretization>& bc)
      : __dirichlet(D),
        __equationNumber(equationNumber),
        __bc(bc)
    {
      ;
    }

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

    ~__SetSecondMemberDirichlet()
    {
      ;
    }
  };


  template <typename MatrixType,
            typename VectorType>
  static void
  __DirichletBorderLinearOperator(const BoundaryConditionDiscretizationElimination<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 SurfaceMesh& surfMesh = *ibcMesh->second;

        BoundaryConditionDiscretizationElimination<MeshType,
                                                   TypeOfDiscretization>
	  ::__meshWrapper(surfMesh,
                          typename BoundaryConditionDiscretizationElimination<MeshType,
                                                                              TypeOfDiscretization>
                          ::__SetSecondMemberDirichlet(D, i, bc));
      }
    }

    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];
      }
    }
  }

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);
    // Dirichlet on the mesh border ...
    BoundaryConditionFDMDiscretization<MeshType,TypeOfDiscretization>
      ::__StandardGetDiagonalDirichletBorderBilinearOperator(*this, z);
  }

  void setMatrix (ReferenceCounting<BaseMatrix> givenA,
                  ReferenceCounting<BaseVector> b) const
  {
    switch(givenA->type()) {
    case BaseMatrix::doubleHashedMatrix: {

      DoubleHashedMatrix& A = dynamic_cast<DoubleHashedMatrix&>(*givenA);

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

      // Dirichlet on the mesh border ...
      BoundaryConditionFDMDiscretization<MeshType,TypeOfDiscretization>
	::__StandardDirichletBorderBilinearOperator(*this, A);
      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);
    BaseMatrix& A = *givenA;

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

    BoundaryConditionDiscretizationElimination<MeshType,TypeOfDiscretization>::
      __DirichletBorderLinearOperator(*this, A, b);
  }

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

    // Dirichlet on the mesh border ...
    BoundaryConditionFDMDiscretization<MeshType,TypeOfDiscretization>::
      __StandardDirichletBorderBilinearOperatorTimesX(*this,
                                                      x, z);
  }

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

    // Dirichlet on the mesh border ...
    BoundaryConditionFDMDiscretization<MeshType,TypeOfDiscretization>::
      __StandardDirichletBorderBilinearOperatorTransposedTimesX(*this,
                                                                x, z);
  }

  BoundaryConditionDiscretizationElimination(const Problem& problem,
                                             const MeshType& mesh,
                                             const DegreeOfFreedomSet& dof)
    : BoundaryConditionFDMDiscretization<MeshType,
                                         TypeOfDiscretization>(problem, mesh, dof)
  {
    ;
  }

  ~BoundaryConditionDiscretizationElimination()
  {
    ;
  }
};

#endif // BOUNDARY_CONDITION_DISCRETIZATION_ELIMINATION_HPP

---