SpectralMethod.cpp

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: SpectralMethod.cpp,v 1.13 2008/02/13 22:05:38 delpinux Exp $

#include <SpectralMethod.hpp>

#include <SolverInformationCenter.hpp>

#include <SpectralLegendreDiscretizationConform.hpp>
#include <BoundaryConditionDiscretizationSpectralConform.hpp>

#include <SpectralLegendreDiscretizationNonConform.hpp>
#include <BoundaryConditionDiscretizationSpectralNonConform.hpp>

#include <PDESolution.hpp>

#include <PDEProblem.hpp>

#include <SpectralMesh.hpp>

#include <KrylovSolver.hpp>

#include <MatrixManagement.hpp>

#include <SparseMatrix.hpp>

#include <Timer.hpp>

#include <ErrorHandler.hpp>

void SpectralMethod::__discretizeOnConformingMesh()
{
  // overwritting memory matrix
  ParameterCenter::instance().set("memory::matrix","none");

  MemoryManager MM;

  MM.ReserveMatrix(__A,
                   problem().numberOfUnknown(),
                   __degreeOfFreedomSet.size());

  MM.ReserveVector(__b,
                   problem().numberOfUnknown(),
                   __degreeOfFreedomSet.size());

  ffout(2) << "Spectral method: disretization...\n";

  ReferenceCounting<SpectralLegendreDiscretizationConform> spectralMethod
    = new SpectralLegendreDiscretizationConform(problem(),
                                                dynamic_cast<const SpectralMesh&>(mesh()),
                                                *__A,*__b, __degreeOfFreedomSet,
                                                __discretizationType);

  spectralMethod->assembleMatrix();
  spectralMethod->assembleSecondMember();

  ffout(2) << "- discretizing boundary conditions\n";

  ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization
    = new BoundaryConditionDiscretizationSpectralConform(problem(),
                                                         dynamic_cast<const SpectralMesh&>(mesh()),
                                                         __degreeOfFreedomSet,
                                                         __discretizationType);

  ffout(2) << "- second member modification\n";
  bcDiscretization->setSecondMember(__A,__b);

  ffout(2) << "- matrix modification\n";
  bcDiscretization->setMatrix(__A,__b);

  ffout(2) << "Spectral method: disretization done\n";
}



void SpectralMethod::__discretizeOnOctreeMesh()
{
  // overwritting memory matrix
  ParameterCenter::instance().set("memory::matrix","none");

  MemoryManager MM;
  
  MM.ReserveMatrix(__A,
                   problem().numberOfUnknown(),
                   __degreeOfFreedomSet.size());

  MM.ReserveVector(__b,
                   problem().numberOfUnknown(),
                   __degreeOfFreedomSet.size());

  ffout(2) << "Spectral method: disretization...\n";
  
  ReferenceCounting<SpectralLegendreDiscretizationNonConform> spectralMethod
    = new SpectralLegendreDiscretizationNonConform(problem(),
                                                   dynamic_cast<const OctreeMesh&>(mesh()),
                                                   *__A,*__b, __degreeOfFreedomSet,
                                                   __discretizationType);

  spectralMethod->assembleMatrix();
  spectralMethod->assembleSecondMember();

  ffout(2) << "- discretizing boundary conditions\n";

  BoundaryConditionDiscretizationSpectralNonConform* bcd
    = new BoundaryConditionDiscretizationSpectralNonConform(problem(),
                                                            dynamic_cast<const OctreeMesh&>(mesh()),
                                                            __degreeOfFreedomSet,
                                                            __discretizationType);
  bcd->associatesMeshesToBoundaryConditions();
  ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization = bcd;

  ffout(2) << "- second member modification\n";
  bcDiscretization->setSecondMember(__A,__b);

  ffout(2) << "- matrix modification\n";
  bcDiscretization->setMatrix(__A,__b);

  ffout(2) << "Spectral method: disretization done\n";
}

void SpectralMethod::__discretize()
{
  switch (mesh().type()) {
  case Mesh::spectralMesh: {
    this->__discretizeOnConformingMesh();
    break;
  }
  case Mesh::octreeMesh: {
    this->__discretizeOnOctreeMesh();
    break;
  }
  default: {
    throw ErrorHandler(__FILE__, __LINE__,
                       "Cannot use '"+mesh().typeName()+"' for spectral method computations",
                       ErrorHandler::normal);
  }
  }
}

void SpectralMethod::Discretize (ConstReferenceCounting<Problem> Pb)
{
  __problem = Pb;

  switch(__discretizationType[0].type()) {
  case ScalarDiscretizationTypeBase::spectralLegendre: {
    this->__discretize();
    return;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Discretization type not implemented",
                       ErrorHandler::normal);
  }
  }
}

void SpectralMethod::Compute (Solution& U)
{
  PDESolution& u = static_cast<PDESolution&>(U);
  KrylovSolver K(*__A, *__b, __degreeOfFreedomSet);
  K.solve(problem(), u.values());
}

SpectralMethod::
SpectralMethod(const DiscretizationType& discretizationType,
               ConstReferenceCounting<Mesh> mesh,
               const DegreeOfFreedomSet& dOfFreedom)
  : Method(discretizationType),
    __mesh(mesh),
    __degreeOfFreedomSet(dOfFreedom)
{
  SolverInformationCenter::instance().pushMesh(mesh);
  SolverInformationCenter::instance().pushDiscretizationType(&discretizationType);
}

SpectralMethod::
~SpectralMethod()
{
  SolverInformationCenter::instance().pop();
}

---