FiniteElementMethod.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: FiniteElementMethod.cpp,v 1.24 2008/02/13 21:24:24 delpinux Exp $

#include <FiniteElementMethod.hpp>

#include <PDESolution.hpp>

#include <PDEProblem.hpp>

#include <Structured3DMesh.hpp>
#include <MeshOfHexahedra.hpp>
#include <MeshOfTetrahedra.hpp>

#include <SpectralMesh.hpp>

#include <FEMDiscretization.hpp>

#include <BoundaryConditionDiscretizationFEM.hpp>

#include <map>
#include <list>

#include <KrylovSolver.hpp>

#include <MatrixManagement.hpp>

#include <SparseMatrix.hpp>
#include <PETScMatrix.hpp>

#include <Timer.hpp>

#include <ErrorHandler.hpp>

template <typename MeshType,
          ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FiniteElementMethod::__discretizeOnMesh()
{
  MemoryManager MM;

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

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

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

  ReferenceCounting<FEMDiscretization<MeshType,
                                      TypeOfDiscretization> > FEM
    = new FEMDiscretization<MeshType,
                            TypeOfDiscretization>(problem(),
                                                  dynamic_cast<const MeshType&>(mesh()),
                                                  __discretizationType,
                                                  *__A,*__b, __degreeOfFreedomSet);

  if (performAssembling) {
    FEM->assembleMatrix();
  } else {
    ffout(2) << "- keeping previous operator discretization\n";
  }

  FEM->assembleSecondMember();


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

  BoundaryConditionDiscretizationFEM<MeshType,
                                     TypeOfDiscretization>* bcd
    = new BoundaryConditionDiscretizationFEM<MeshType,
                                             TypeOfDiscretization>(problem(),
                                                                   dynamic_cast<const MeshType&>(mesh()),
                                                                   __degreeOfFreedomSet);
  bcd->associatesMeshesToBoundaryConditions();
  ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization = bcd;

  // Set Dirichlet information to the matrix
  FEM->setDirichletList(bcDiscretization->getDirichletList());

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

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

  ffout(2) << "Finite element method: disretization done\n";

  if (__A->type() == BaseMatrix::doubleHashedMatrix) {
    Timer t;
    t.start();

#warning temporary implementation
#ifdef    HAVE_PETSC
    PETScMatrix* aa
      = new PETScMatrix(static_cast<DoubleHashedMatrix&>(*__A));
    __A = aa; // now use sparse matrix
#else  // HAVE_PETSC
    SparseMatrix* aa
      = new SparseMatrix(static_cast<DoubleHashedMatrix&>(*__A));
    
    __A = aa; // now use sparse matrix
#endif // HAVE_PETSC

    t.stop();
    ffout(2) << "Matrix copy: " << t << '\n';
  }
}

template <ScalarDiscretizationTypeBase::Type TypeOfDiscretization>
void FiniteElementMethod::__discretize()
{
  switch (mesh().type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__discretizeOnMesh<Structured3DMesh, TypeOfDiscretization>();
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__discretizeOnMesh<MeshOfHexahedra, TypeOfDiscretization>();
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__discretizeOnMesh<MeshOfTetrahedra, TypeOfDiscretization>();
    break;
  }
  case Mesh::spectralMesh: {
    this->__discretizeOnMesh<SpectralMesh, TypeOfDiscretization>();
    break;
  }
  default: {
    throw ErrorHandler(__FILE__, __LINE__,
                       "Cannot use '"+mesh().typeName()+"' for finite element computations",
                       ErrorHandler::normal);
  }
  }
}

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

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

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

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

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

---