BoundaryConditionDiscretizationSpectralConform.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: BoundaryConditionDiscretizationSpectralConform.hpp,v 1.2 2008/03/10 21:58:11 delpinux Exp $

#ifndef BOUNDARY_CONDITION_DISCRETIZATION_SPECTRAL_CONFORM_HPP
#define BOUNDARY_CONDITION_DISCRETIZATION_SPECTRAL_CONFORM_HPP

#include <Problem.hpp>
#include <DegreeOfFreedomSet.hpp>

#include <BoundaryConditionDiscretization.hpp>

#include <DiscretizationType.hpp>

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

#include <SpectralMesh.hpp>
#include <Interval.hpp>

#include <GaussLobattoManager.hpp>

#include <LegendreBasis.hpp>
#include <SpectralConformTransformation.hpp>

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

class BoundaryConditionDiscretizationSpectralConform
  : public BoundaryConditionDiscretization
{
private:
  const DegreeOfFreedomSet& __degreeOfFreedomSet;
  const SpectralMesh& __mesh;
  
  Vector<TinyVector<3,size_t> >  __dofShape;

  TinyVector<3, ConstReferenceCounting<Interval> > __interval;
  Vector<TinyVector<3, ConstReferenceCounting<Interval> > > __intervalU;
  TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > __transform;
  Vector<TinyVector<3, ConstReferenceCounting<SpectralConformTransformation> > > __transformU;
  TinyVector<3, ConstReferenceCounting<GaussLobatto> > __gaussLobatto;
  Vector<TinyVector<3, ConstReferenceCounting<LegendreBasis>  > > __basisU;
  
  
public:
  void getDiagonal(BaseVector& X) const;
  
  void setMatrix(ReferenceCounting<BaseMatrix> A,
                 ReferenceCounting<BaseVector> b) const;
  
  void setSecondMember(ReferenceCounting<BaseMatrix> A,
                       ReferenceCounting<BaseVector> b) const;
  
  void timesX(const BaseVector& X, BaseVector& Z) const;
  
  void transposedTimesX(const BaseVector& X, BaseVector& Z) const;
  
  
  BoundaryConditionDiscretizationSpectralConform(const Problem& problem,
                                                 const SpectralMesh& mesh,
                                                 const DegreeOfFreedomSet& dof,
                                                 const DiscretizationType& discretizationType)
    : BoundaryConditionDiscretization(problem,dof),
      __degreeOfFreedomSet(dof),
      __mesh(mesh),
      __dofShape(this->problem().numberOfUnknown()),
      __interval(new Interval(__mesh.shape().a()[0],__mesh.shape().b()[0]),
                 new Interval(__mesh.shape().a()[1],__mesh.shape().b()[1]),
                 new Interval(__mesh.shape().a()[2],__mesh.shape().b()[2])),
      
      __intervalU(this->problem().numberOfUnknown()),
      __transform(new SpectralConformTransformation(*__interval[0]),
                  new SpectralConformTransformation(*__interval[1]),
                  new SpectralConformTransformation(*__interval[2])),
      
      __transformU(this->problem().numberOfUnknown()),
      
      __gaussLobatto(GaussLobattoManager::instance().getReference(__mesh.degree(0)+1),
                     GaussLobattoManager::instance().getReference(__mesh.degree(1)+1),
                     GaussLobattoManager::instance().getReference(__mesh.degree(2)+1)),
      __basisU(this->problem().numberOfUnknown())
  {
    for (size_t i=0; i < this->problem().numberOfUnknown(); i++) {
      if (discretizationType[i].type() != ScalarDiscretizationTypeBase::spectralLegendre) {
        throw ErrorHandler(__FILE__,__LINE__,
                           "discretization '"
                           +ScalarDiscretizationTypeBase::name(discretizationType[i])
                           +"' is incompatible with spectral method",
                           ErrorHandler::unexpected);
      }
      const ScalarDiscretizationTypeSpectral& discretization
        = dynamic_cast<const ScalarDiscretizationTypeSpectral&>(discretizationType[i]);
      
      __dofShape[i] = discretization.degrees()+TinyVector<3,size_t>(1,1,1);
      
      __intervalU[i][0]= new Interval(discretization.a()[0],discretization.b()[0]);
      __intervalU[i][1]= new Interval(discretization.a()[1],discretization.b()[1]);
      __intervalU[i][2]= new Interval(discretization.a()[2],discretization.b()[2]);
      __transformU[i][0] = new SpectralConformTransformation(*__intervalU[i][0]);
      __transformU[i][1] = new SpectralConformTransformation(*__intervalU[i][1]);
      __transformU[i][2] = new SpectralConformTransformation(*__intervalU[i][2]); 
      __basisU[i][0]= new LegendreBasis(discretization.degrees()[0]);
      __basisU[i][1]= new LegendreBasis(discretization.degrees()[1]);
      __basisU[i][2]= new LegendreBasis(discretization.degrees()[2]);
    }
    
  }
  
  ~BoundaryConditionDiscretizationSpectralConform()
  {
    ;
  }
};

#endif // BOUNDARY_CONDITION_DISCRETIZATION_SPECTRAL_CONFORM_HPP

---