//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane 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: FictitiousDomainMethod.cpp,v 1.26 2007/05/20 23:15:36 delpinux Exp $

#include <FictitiousDomainMethod.hpp>

#include <KrylovSolver.hpp>
#include <PDESolution.hpp>

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

#include <PDESystem.hpp>

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

#include <Dirichlet.hpp>
#include <Fourrier.hpp>
#include <Neumann.hpp>

#include <MatrixManagement.hpp>

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

#include <FEMDiscretization.hpp>

#include <FEMFunctionBuilder.hpp>

#include <ConformTransformation.hpp>

#include <BoundaryConditionDiscretization.hpp>

void FictitiousDomainMethod::Compute(Solution& U)
{
  Structured3DMesh& m = dynamic_cast<Structured3DMesh&>(mesh());

  KrylovSolver K(*__A, *__b, __degreeOfFreedomSet,
		 new Structured3DMeshShape(m.shape()));
  
  PDESolution& u = static_cast<PDESolution&>(U);

  K.solve(problem(), u.values());
  for (size_t i=0; i<m.numberOfCells();++i) {
    Cell& k = m.cell(i);
    k.setFictitious(false); // reset fictitious state!
 }
}


template <typename CellType>
real_t FictitiousDomainMethod::
__integrateCharacteristicFunction(const CellType& k, const Domain& d)
{
  throw ErrorHandler(__FILE__,__LINE__,
		     "not implemented",
		       ErrorHandler::unexpected);
  return 0;
}

template <>
real_t FictitiousDomainMethod::
__integrateCharacteristicFunction(const CartesianHexahedron& k, const Domain& d)
{
  ConformTransformationQ1CartesianHexahedron T(k);
  return T.integrateCharacteristic(d);
}

template <typename MeshType>
void FictitiousDomainMethod::
__computesDegreesOfFreedom(ConstReferenceCounting<Problem> givenProblem)
{
  ffout(2) << "Adding Characteristic function of the domain to PDEs\n";
  Vector<bool> dofVertices;

  dofVertices.resize(mesh().numberOfVertices());
  dofVertices = false;

  const Domain& domain = *givenProblem->domain();

  for (size_t i=0; i<mesh().numberOfVertices(); ++i) {
    dofVertices[i] = domain.inside(mesh().vertex(i));
  }

  Vector<real_t> kiValues(mesh().numberOfCells()); // "mean" value of the characteristic function

  MeshType& m = static_cast<MeshType&>(*__mesh);

  for (size_t i=0; i<m.numberOfCells();++i) {
    typename MeshType::CellType& k = m.cell(i);
    unsigned numberIn = 0;

    // Computation of the caracteristic function
    for (unsigned j=0; j<MeshType::CellType::NumberOfVertices; ++j) {
      numberIn += (dofVertices(m.vertexNumber(k(j))))?1:0;
    }

    switch (numberIn) {
    case 0: {
      k.setFictitious(true);
      kiValues[i] = 0;
      break;
    }
    case MeshType::CellType::NumberOfVertices: {
      kiValues[i] = 1;
      break;
    }
    default:{
      kiValues[i] = __integrateCharacteristicFunction(k, domain);
    }
    }
  }

  DiscretizationType discretizationType(DiscretizationType::lagrangianFEM0);

  FEMFunctionBuilder femBuilder;
  femBuilder.build(discretizationType,
		   __mesh,
		   kiValues);
  ConstReferenceCounting<ScalarFunctionBase> u = femBuilder.getBuiltScalarFunction();

  __problem = (*givenProblem) * u;
  ffout(2) << "Adding Characteristic function of the domain to PDEs: done\n";

  ffout(2) << "Using "<< this->__degreeOfFreedomSet.numberOfUsedDOFPositions()
	   << " DOF over "
	   << this->__degreeOfFreedomSet.numberOfDOFPositions() << " available!\n";
}

void FictitiousDomainMethod::
computesDegreesOfFreedom(ConstReferenceCounting<Problem> givenProblem)
{
  switch (mesh().type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__computesDegreesOfFreedom<Structured3DMesh>(givenProblem);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "mesh type not supported by FDM",
		       ErrorHandler::normal);
  }
  }
}


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

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

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

  ffout(2) << "Fictitious domain method: discretization...\n";

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

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

  FEM->assembleSecondMember();

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

  ReferenceCounting<BoundaryConditionDiscretization> bcDiscretization
    = this->discretizeBoundaryConditions();

  // 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) << "Fictitious domain method: discretization 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 time: " << t << '\n';
  }
}

template <DiscretizationType::Type TypeOfDiscretization>
void FictitiousDomainMethod::__discretize()
{

  switch (this->mesh().type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__discretizeOnMesh<Structured3DMesh, TypeOfDiscretization>();
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "this mesh type is not supported by FDM",
		       ErrorHandler::normal);
  }
  }

}

void FictitiousDomainMethod::Discretize(ConstReferenceCounting<Problem> givenProblem)
{
  this->computesDegreesOfFreedom(givenProblem);

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