BoundaryConditionDiscretizationSpectralNonConform.cpp

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: BoundaryConditionDiscretizationSpectralNonConform.cpp,v 1.3 2008/03/02 22:40:35 delpinux Exp $

#include <BoundaryConditionDiscretizationSpectralNonConform.hpp>
#include <ErrorHandler.hpp>

#include <VariationalProblem.hpp>
#include <DegreeOfFreedomSet.hpp>
#include <VariationalBorderOperatorAlphaUV.hpp>
#include <VariationalBorderOperatorFV.hpp>
#include <VariationalLinearOperator.hpp>
#include <SpectralConformTransformation.hpp>

#include <BaseMatrix.hpp>
#include <BaseVector.hpp>

#include <Discretization.hpp>
#include <ScalarDiscretizationTypeSpectral.hpp>

#include <BoundaryReferences.hpp>

#include <UnAssembledMatrix.hpp>

#include <DoubleHashedMatrix.hpp>

#include <ConformTransformation.hpp>

//typedef  QuadratureFormulaP1Triangle3D() quadratureType;

void BoundaryConditionDiscretizationSpectralNonConform::
getDiagonal(BaseVector& X) const
{
  Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(X);
  
  const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bcMeshAssociation;
  
  for (BoundaryConditionSurfaceMeshAssociation
         ::BilinearBorderOperatorMeshAssociation
         ::const_iterator i = bcMeshAssociation.beginOfBilinear();
       i != bcMeshAssociation.endOfBilinear(); ++i) {
    
    switch (i->first->type()) {
    case VariationalBilinearBorderOperator::alphaUV: {
      const VariationalBorderOperatorAlphaUV& operatorAlphaUV
        = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);
      if (operatorAlphaUV.unknownNumber() == operatorAlphaUV.testFunctionNumber()) {

        const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();   

        size_t unknownNumber = operatorAlphaUV.unknownNumber();
        const ScalarDiscretizationTypeSpectral& discretization
          = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(__discretizationType[unknownNumber]);
      
        TinyVector<3, ConstReferenceCounting<Interval> >  interval;
        TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > spectralTransform;
        TinyVector<3, ConstReferenceCounting<LegendreBasis> > basis;

        for(size_t n=0; n<3; ++n){
          interval[n]= new Interval(discretization.a()[n],discretization.b()[n]);
          spectralTransform[n] = new SpectralConformTransformation(*interval[n]);
          basis[n]= new LegendreBasis(discretization.degrees()[n]);
        }
        
        const Mesh& boundaryMesh = (*i->second);
        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          typedef  QuadratureFormulaP1Triangle3D QuadratureType;

          const QuadratureType& quadrature = QuadratureType::instance();
          const TinyVector<QuadratureType::numberOfQuadraturePoints, TinyVector<3,real_t> >& 
            nodesRef =  quadrature.vertices() ;

        for (SurfaceMeshOfTriangles::const_iterator  itriangle(dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
             not(itriangle.end()); ++ itriangle) {
          const Triangle& T = * itriangle;

          TinyVector<3,Vector<real_t> > baseValues; 
          for(size_t j = 0; j < 3; j++){
            baseValues[j].resize(discretization.degrees()[j] + 1);
          }
          
          ConformTransformationP1Triangle::ConformTransformationP1Triangle transform(T);

          for (size_t q=0; q< QuadratureType::numberOfQuadraturePoints; ++q) {
        
            //Transformation of the points of Gauss in Triangle by ConformTransformation and spectralConformTransformation
            
            TinyVector<3,real_t> nodes;
            transform.value(nodesRef[q],nodes);
            const real_t xi= (spectralTransform[0])->inverse(nodes[0]);
            const real_t yi= (spectralTransform[1])->inverse(nodes[1]);
            const real_t zi= (spectralTransform[2])->inverse(nodes[2]);
            (basis[0])->getValues(xi,baseValues[0]);
            (basis[1])->getValues(yi,baseValues[1]);
            (basis[2])->getValues(zi,baseValues[2]);

            const real_t weight = quadrature.weight(q) ;  
            const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];

            for (size_t i=0; i< (basis[0])->dimension(); ++i){
              const real_t q1 = alpha(nodes)* weight * T.volume() * baseValues[0][i]* baseValues[0][i];
              for (size_t j=0; j<( basis[1])->dimension(); ++j){
                const real_t q2 = q1 * baseValues[1][j]*baseValues[1][j];
                for (size_t k=0; k< (basis[2])->dimension(); ++k){
                  const real_t q3 = q2 * baseValues[2][k] * baseValues[2][k];
                  
                  size_t dofNumberU  = i*dofShapeU[1]*dofShapeU[2] + j*dofShapeU[2] + k;

                  v[__degreeOfFreedomSet(unknownNumber,dofNumberU)] += q3;
                }
              }
            }
          }
        }
//        throw ErrorHandler(__FILE__,__LINE__,
//                           "BC not implemented in triangle meshes",
//                           ErrorHandler::unexpected);
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "non conforming boundary condition is not implemented in quadrangle meshes",
                             ErrorHandler::normal);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented",
                         ErrorHandler::unexpected);
    }
    }
  }
}


void BoundaryConditionDiscretizationSpectralNonConform::
setMatrix(ReferenceCounting<BaseMatrix> givenA,
          ReferenceCounting<BaseVector> b) const
{
  switch((*givenA).type()) {
  case BaseMatrix::doubleHashedMatrix: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);       
    break;
  }
  case BaseMatrix::unAssembled: {
    UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(*givenA);
    A.setBoundaryConditions(this);
    break;
  }

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


void BoundaryConditionDiscretizationSpectralNonConform::
setSecondMember(ReferenceCounting<BaseMatrix> givenA,
                ReferenceCounting<BaseVector> givenb) const
{
  ffout(2)<<"- assembling second membre-boundary\n";
  Vector<real_t>& b= (static_cast<Vector<real_t>&>(*givenb));

  const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bcMeshAssociation;

  for (BoundaryConditionSurfaceMeshAssociation
         ::LinearBorderOperatorMeshAssociation
         ::const_iterator i = bcMeshAssociation.beginOfLinear();
       i != bcMeshAssociation.endOfLinear(); ++i) {
    
    switch (i->first->type()) {
    case VariationalLinearBorderOperator::FV: {
      const VariationalBorderOperatorFV& operatorFV
        = dynamic_cast<const VariationalBorderOperatorFV&>(*i->first);
      const ScalarFunctionBase& f = operatorFV.f();
      
      size_t testFunctionNumber = operatorFV.testFunctionNumber();
      
      const ScalarDiscretizationTypeSpectral& discretization
        = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(__discretizationType[testFunctionNumber]);
      
      TinyVector<3, ConstReferenceCounting<Interval> >  interval;
      TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > spectralTransform;
      TinyVector<3, ConstReferenceCounting<LegendreBasis> > basis;
      
      for(size_t n=0; n<3; ++n){
        interval[n]= new Interval(discretization.a()[n],discretization.b()[n]);
        spectralTransform[n] = new SpectralConformTransformation(*interval[n]);
        basis[n]= new LegendreBasis(discretization.degrees()[n]);
      }
      const Mesh& boundaryMesh = (*i->second);  
      switch(boundaryMesh.type()) {
      case Mesh::surfaceMeshTriangles: {
        
        typedef  QuadratureFormulaP1Triangle3D QuadratureType;

        const QuadratureType& quadrature = QuadratureType::instance();
          
        const TinyVector<QuadratureType::numberOfQuadraturePoints, TinyVector<3,real_t> >& 
          nodesRef =  quadrature.vertices() ;
          
        for (SurfaceMeshOfTriangles::const_iterator  itriangle(dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
             not(itriangle.end()); ++ itriangle) {
          const Triangle& T = *itriangle;
          
          TinyVector<3,Vector<real_t> > baseValues;
          for(size_t j = 0; j < 3; j++){
            baseValues[j].resize(discretization.degrees()[j] + 1);
          }
          const TinyVector<3,size_t>& dofShapeV =__dofShape[operatorFV.testFunctionNumber()];      
          ConformTransformationP1Triangle::ConformTransformationP1Triangle transform(T);
          
          for (size_t q=0; q< QuadratureType::numberOfQuadraturePoints; ++q) {

            //Transformation of the points of Gauss in Triangle by ConformTransformation and spectralConformTransformation
            TinyVector<3,real_t> nodes;
            
            transform.value(nodesRef[q],nodes);
            const real_t weight = quadrature.weight(q);
            const real_t xi= (spectralTransform[0])->inverse(nodes[0]);
            const real_t yi= (spectralTransform[1])->inverse(nodes[1]);
            const real_t zi= (spectralTransform[2])->inverse(nodes[2]);
            
            (basis[0])->getValues(xi,baseValues[0]);
            (basis[1])->getValues(yi,baseValues[1]);
            (basis[2])->getValues(zi,baseValues[2]);
            
            for (size_t m=0; m< basis[0]->dimension(); ++m){
              const real_t b_m =T.volume()* weight * f(nodes)* baseValues[0][m] ;
              for (size_t n=0; n< basis[1]->dimension(); ++n){
                const real_t b_mn = b_m * baseValues[1][n];
                for (size_t p=0; p< basis[2]->dimension(); ++p){
                  size_t dofNumberV = m*dofShapeV[1]*dofShapeV[2] + n*dofShapeV[2] + p;
                  b[__degreeOfFreedomSet(testFunctionNumber,dofNumberV)]
                    += b_mn* baseValues[2][p];
                }
              }
            }
          }
        }
        
        break;
      }
      case Mesh::surfaceMeshQuadrangles: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "non conforming boundary condition is not implemented in quadrangle meshes",
                           ErrorHandler::normal);
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "not implemented",
                           ErrorHandler::unexpected);
      }
      }
        break;
      }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented",
                         ErrorHandler::unexpected);
    }
    }
    }
}


void BoundaryConditionDiscretizationSpectralNonConform::
timesX(const BaseVector& X, BaseVector& Z) const
{
  const Vector<real_t>& u = dynamic_cast<const Vector<real_t>&>(X);
  Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(Z);

  const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bcMeshAssociation;
  
  for (BoundaryConditionSurfaceMeshAssociation
         ::BilinearBorderOperatorMeshAssociation
         ::const_iterator i = bcMeshAssociation.beginOfBilinear();
       i != bcMeshAssociation.endOfBilinear(); ++i) {
    
    switch (i->first->type()) {
    case VariationalBilinearBorderOperator::alphaUV: {
      const VariationalBorderOperatorAlphaUV& operatorAlphaUV
        = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);
      const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();   
      
      size_t unknownNumber = operatorAlphaUV.unknownNumber();   
      size_t testFunctionNumber = operatorAlphaUV.testFunctionNumber();
      
      const ScalarDiscretizationTypeSpectral& discretization
        = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(__discretizationType[unknownNumber]);
      
      TinyVector<3, ConstReferenceCounting<Interval> >  interval;
      TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > spectralTransform;
      TinyVector<3, ConstReferenceCounting<LegendreBasis> > basis;
      
      for(size_t n=0;n<3;++n){
        interval[n]= new Interval(discretization.a()[n],discretization.b()[n]);
        spectralTransform[n] = new SpectralConformTransformation(*interval[n]);
        basis[n]= new LegendreBasis(discretization.degrees()[n]);
      }      
      
      const Mesh& boundaryMesh = (*i->second);
      switch(boundaryMesh.type()) {
      case Mesh::surfaceMeshTriangles: {
        typedef  QuadratureFormulaP1Triangle3D QuadratureType;
        
        const QuadratureType& quadrature = QuadratureType::instance();
        const TinyVector<QuadratureType::numberOfQuadraturePoints, TinyVector<3,real_t> >& 
          nodesRef =  quadrature.vertices() ;
          
        for (SurfaceMeshOfTriangles::const_iterator  itriangle(dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
             not(itriangle.end()); ++ itriangle) {
          const Triangle& T = * itriangle;

          TinyVector<3,Vector<real_t> > baseValues; 
          for(size_t j = 0; j < 3; j++){
            baseValues[j].resize(discretization.degrees()[j] + 1);
          }
          
          ConformTransformationP1Triangle::ConformTransformationP1Triangle transform(T);

          for (size_t q=0; q< QuadratureType::numberOfQuadraturePoints; ++q) {
        
            //Transformation of the points of Gauss in Triangle by ConformTransformation and spectralConformTransformation
            
            TinyVector<3,real_t> nodes;
            transform.value(nodesRef[q],nodes);
            const real_t xi= (spectralTransform[0])->inverse(nodes[0]);
            const real_t yi= (spectralTransform[1])->inverse(nodes[1]);
            const real_t zi= (spectralTransform[2])->inverse(nodes[2]);
            (basis[0])->getValues(xi,baseValues[0]);
            (basis[1])->getValues(yi,baseValues[1]);
            (basis[2])->getValues(zi,baseValues[2]);

            const real_t weight = quadrature.weight(q) ;  
            const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];

            real_t u_q = 0;

            for (size_t i=0; i< (basis[0])->dimension(); ++i){
              const real_t q1 = alpha(nodes)* weight * T.volume() * baseValues[0][i];
              for (size_t j=0; j<( basis[1])->dimension(); ++j){
                const real_t q2 = q1* baseValues[1][j];
                for (size_t k=0; k< (basis[2])->dimension(); ++k){
                  const real_t q3 = baseValues[2][k] *q2;
                  size_t dofNumberU  = i*dofShapeU[1]*dofShapeU[2] + j*dofShapeU[2] + k;
                  u_q += q3 * u[__degreeOfFreedomSet(unknownNumber,dofNumberU)];
                }
              }
            }
            const TinyVector<3,size_t>& dofShapeV = __dofShape[testFunctionNumber];        
            
            for (size_t m=0; m< basis[0]->dimension(); ++m){
              const real_t v_qm = u_q* baseValues[0][m];
              for (size_t n=0; n< basis[1]->dimension(); ++n){
                const real_t v_mn = v_qm * baseValues[1][n];
                for (size_t p=0; p< basis[2] ->dimension(); ++p){
                  size_t dofNumberV = m*dofShapeV[1]*dofShapeV[2] + n*dofShapeV[2] + p;
                  v[__degreeOfFreedomSet(testFunctionNumber,dofNumberV)]
                    += baseValues[2][p]* v_mn;
                }
              }
            }
          }
        }
        break;
      }
        case Mesh::surfaceMeshQuadrangles: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "non conforming boundary condition is not implemented in quadrangle meshes",
                             ErrorHandler::normal);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
      
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented",
                         ErrorHandler::unexpected);
    }
    }
  }
}

void BoundaryConditionDiscretizationSpectralNonConform::
transposedTimesX(const BaseVector& X, BaseVector& Z) const
{
  const Vector<real_t>& u = dynamic_cast<const Vector<real_t>&>(X);
  Vector<real_t>& v = dynamic_cast<Vector<real_t>&>(Z);

  const BoundaryConditionSurfaceMeshAssociation& bcMeshAssociation = *__bcMeshAssociation;
  
  for (BoundaryConditionSurfaceMeshAssociation
         ::BilinearBorderOperatorMeshAssociation
         ::const_iterator i = bcMeshAssociation.beginOfBilinear();
       i != bcMeshAssociation.endOfBilinear(); ++i) {
    
    switch (i->first->type()) {
    case VariationalBilinearBorderOperator::alphaUV: {
      const VariationalBorderOperatorAlphaUV& operatorAlphaUV
        = dynamic_cast<const VariationalBorderOperatorAlphaUV&>(*i->first);
      if (operatorAlphaUV.unknownNumber() == operatorAlphaUV.testFunctionNumber()) {
        
        const ScalarFunctionBase& alpha = operatorAlphaUV.alpha();   
        
        size_t unknownNumber = operatorAlphaUV.unknownNumber();   
        size_t testFunctionNumber = operatorAlphaUV.testFunctionNumber();
        
        const ScalarDiscretizationTypeSpectral& discretization
          = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(__discretizationType[unknownNumber]);
        
        TinyVector<3, ConstReferenceCounting<Interval> >  interval;
        TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > spectralTransform;
        TinyVector<3, ConstReferenceCounting<LegendreBasis> > basis;
        
        for(size_t n=0;n<3;++n){
          interval[n]= new Interval(discretization.a()[n],discretization.b()[n]);
          spectralTransform[n] = new SpectralConformTransformation(*interval[n]);
          basis[n]= new LegendreBasis(discretization.degrees()[n]);
        }      
        
        const Mesh& boundaryMesh = (*i->second);
        switch(boundaryMesh.type()) {
        case Mesh::surfaceMeshTriangles: {
          typedef  QuadratureFormulaP1Triangle3D QuadratureType;
          
          const QuadratureType& quadrature = QuadratureType::instance();
          const TinyVector<QuadratureType::numberOfQuadraturePoints, TinyVector<3,real_t> >& 
            nodesRef =  quadrature.vertices() ;
          
          for (SurfaceMeshOfTriangles::const_iterator  itriangle(dynamic_cast<const SurfaceMeshOfTriangles&>(boundaryMesh));
               not(itriangle.end()); ++ itriangle) {
            const Triangle& T = * itriangle;
            TinyVector<3,Vector<real_t> > baseValues; 
            for(size_t j = 0; j < 3; j++){
              baseValues[j].resize(discretization.degrees()[j] + 1);
            }
            
            ConformTransformationP1Triangle::ConformTransformationP1Triangle transform(T);
            
            for (size_t q=0; q< QuadratureType::numberOfQuadraturePoints; ++q) {
              
              //Transformation of the points of Gauss in Triangle by ConformTransformation and spectralConformTransformation
              
              TinyVector<3,real_t> nodes;
              transform.value(nodesRef[q],nodes);
              const real_t xi= (spectralTransform[0])->inverse(nodes[0]);
              const real_t yi= (spectralTransform[1])->inverse(nodes[1]);
              const real_t zi= (spectralTransform[2])->inverse(nodes[2]);
              (basis[0])->getValues(xi,baseValues[0]);
              (basis[1])->getValues(yi,baseValues[1]);
              (basis[2])->getValues(zi,baseValues[2]);
              
              const real_t weight = quadrature.weight(q);
              const TinyVector<3,size_t>& dofShapeU = __dofShape[unknownNumber];
              real_t u_q =0;
              
              for (size_t i=0; i<(basis[0])->dimension(); ++i){
                const real_t q1 =alpha(nodes)* weight * T.volume() * baseValues[0][i];
                for (size_t j=0; j<(basis[1])->dimension(); ++j){
                  const real_t q2 = q1 * baseValues[1][j];
                  for (size_t k=0; k<(basis[2])->dimension(); ++k){
                    const real_t q3 = baseValues[2][k]* q2;
                    size_t dofNumberU  = i*dofShapeU[1]*dofShapeU[2] + j*dofShapeU[2] + k;
                    u_q += q3*u[__degreeOfFreedomSet(testFunctionNumber,dofNumberU)];
                    
                  }
                }
              }
              const TinyVector<3,size_t>& dofShapeV= __dofShape[testFunctionNumber];
              
              for (size_t m=0; m<(basis[0])->dimension(); ++m){
                const real_t v_qm = u_q*baseValues[0][m];
                for (size_t n=0; n<(basis[1])->dimension(); ++n){
                  const real_t v_mn = v_qm * baseValues[1][n];
                  for (size_t p=0; p<(basis[2])->dimension(); ++p){
                    size_t dofNumberV = m*dofShapeV[1]*dofShapeV[2] + n*dofShapeV[2] + p;
                    v[__degreeOfFreedomSet(unknownNumber,dofNumberV)]
                      += baseValues[2][p]* v_mn;
                  }
                }
              }
            }
          }       
          break;
        }
        case Mesh::surfaceMeshQuadrangles: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "non conforming boundary condition is not implemented in quadrangle meshes",
                             ErrorHandler::normal);
          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "not implemented",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "not implemented",
                         ErrorHandler::unexpected);
    }
    }
  }
}

---