SpectralLegendreDiscretizationNonConform.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2007 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: SpectralLegendreDiscretizationNonConform.hpp,v 1.4 2008/02/13 21:20:25 delpinux Exp $

#ifndef SPECTRAL_LEGENDRE_DISCRETIZATION_NON_CONFORM_HPP
#define SPECTRAL_LEGENDRE_DISCRETIZATION_NON_CONFORM_HPP

#include <Problem.hpp>

#include <DegreeOfFreedomSet.hpp>
#include <OctreeMesh.hpp>

#include <DoubleHashedMatrix.hpp>
#include <UnAssembledMatrix.hpp>

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

#include <Mesh.hpp>

#include <SpectralLegendreDiscretizer.hpp>

#include <ErrorHandler.hpp>
#include <Timer.hpp>

class SpectralLegendreDiscretizationNonConform
  : public Discretization
{
private:
  // Set of degrees of freedom
  const DegreeOfFreedomSet& __degreeOfFreedomSet;

  const OctreeMesh& __mesh;
  const DiscretizationType& __discretizationType;
  const TinyVector<3,size_t> __quadratureDegree;

  static TinyVector<3,size_t>
  _getQuadratureDegree(const DiscretizationType& discretizationType)
  {
#ifndef NDEBUG
    for (size_t i=0; i<discretizationType.number(); ++i) {
      ASSERT(discretizationType[i].type() == ScalarDiscretizationTypeBase::spectralLegendre);
    }
#endif // NDEBUG

    TinyVector<3,size_t> quadratureDegree
      = dynamic_cast <const ScalarDiscretizationTypeSpectral&> (discretizationType[0]).degrees();

    for (size_t i=1; i<discretizationType.number(); ++i) {
      TinyVector<3,size_t> degree
        = dynamic_cast <const ScalarDiscretizationTypeSpectral&>(discretizationType[i]).degrees();
      for (size_t j=0; j<3; ++j) {
        quadratureDegree[j] = std::max(quadratureDegree[j], degree[j]);
      }
    }

    return quadratureDegree;
  }

  ReferenceCounting<SpectralLegendreDiscretizer>
  _getLocalDiscretizer(const CartesianHexahedron& h) const
  {
      const TinyVector<3,real_t>& a = h(0);
      const TinyVector<3,real_t>& b = h(6);

      Structured3DMeshShape s3dM(__quadratureDegree,a,b);
      ReferenceCounting<VerticesCorrespondance> correspondance
        = new VerticesCorrespondance((__quadratureDegree[0]+1)*
                                     (__quadratureDegree[1]+1)*
                                     (__quadratureDegree[2]+1));

      SpectralMesh quadratureMesh(s3dM,correspondance);

      return new SpectralLegendreDiscretizer(this->problem(),
                                             quadratureMesh,
                                             __A,
                                             __b,
                                             __degreeOfFreedomSet,
                                             __discretizationType);
  }

public:  
  void assembleMatrix()
  {
    switch ((this->__A).type()) {
    case BaseMatrix::doubleHashedMatrix: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Spectral Method cannot be used with assembled matrices",
                         ErrorHandler::unexpected);
    }
    case BaseMatrix::unAssembled: {
      UnAssembledMatrix& A = dynamic_cast<UnAssembledMatrix&>(this->__A);
      A.setDiscretization(this);
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected matrix type",
                         ErrorHandler::unexpected);
    }
    }
  }
  
  void timesX(const BaseVector& u, BaseVector& v) const
  {
    dynamic_cast<Vector<real_t>&>(v)=0;
    for (OctreeMesh::const_iterator icell(__mesh);
         not icell.end(); ++icell) {
      const CartesianHexahedron& h = *icell;

      this->_getLocalDiscretizer(h)->timesX(u,v);
    }
  }
  
  void transposedTimesX(const BaseVector& u, BaseVector& v) const
  {
    dynamic_cast<Vector<real_t>&>(v) = 0;
    for (OctreeMesh::const_iterator icell(__mesh);
         not icell.end(); ++icell) {
      const CartesianHexahedron& h = *icell;

      this->_getLocalDiscretizer(h)->transposedTimesX(u,v);
    }
  }

  void getDiagonal(BaseVector& z) const
  {
    dynamic_cast<Vector<real_t>&>(z)=0;
    for (OctreeMesh::const_iterator icell(__mesh);
         not icell.end(); ++icell) {
      const CartesianHexahedron& h = *icell;

      this->_getLocalDiscretizer(h)->getDiagonal(z);
    }
  }

  void assembleSecondMember()
  {
    Timer t;
    t.start();
    
    ffout(2) << "- assembling second member\n";
    
    Vector<real_t>& b = (static_cast<Vector<real_t>&>(this->__b));
    b = 0;

    for (OctreeMesh::const_iterator icell(__mesh);
         not icell.end(); ++icell) {
      const CartesianHexahedron& h = *icell;

      this->_getLocalDiscretizer(h)->assembleSecondMember();
    }

    ffout(2) << "- assembling second member: done";
    t.stop();
    ffout(3) << " [cost: " << t << ']';
    ffout(2) << '\n';
  }
  
public:
  
  SpectralLegendreDiscretizationNonConform(const Problem& p,
                                           const OctreeMesh& m,
                                           BaseMatrix& a,
                                           BaseVector& bb,
                                           const DegreeOfFreedomSet& dof,
                                           const DiscretizationType& discretizationType)
    : Discretization(discretizationType,
                     p, a, bb),
      __degreeOfFreedomSet(dof),
      __mesh(m),
      __discretizationType(discretizationType),
      __quadratureDegree(SpectralLegendreDiscretizationNonConform::_getQuadratureDegree(discretizationType))
  {
    ;
  }

  ~SpectralLegendreDiscretizationNonConform()
  {
    ;
  }
};

    
#endif // SPECTRAL_LEGENDRE_DISCRETIZATION_NON_CONFORM_HPP

---