//  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: FictitiousDomainMethod.cpp,v 1.7 2004/03/15 21:01:20 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 <FEMDiscretization.hpp>

#include <ConformTransformation.hpp>

#include <BoundaryConditionDiscretization.hpp>

void FictitiousDomainMethod::Compute(Solution& U)
{
  ffout(2) << "Solving the linear system\n";

  KrylovSolver K(*__A, *__b, *__degreeOfFreedomSet,
		 new Structured3DMeshShape(dynamic_cast<Structured3DMesh&>(mesh()).shape()));
  
  PDESolution& u = static_cast<PDESolution&>(U);
  Vector<real_t>& values = dynamic_cast<Vector<real_t>&>(u.values());

  for (size_t i=0; i<__givenDegreeOfFreedomSet.numberOfVertices(); ++i) {
    if ((*__degreeOfFreedomSet).isDOFVertex(i)) {
      for (size_t j=0; j<__givenDegreeOfFreedomSet.numberOfVariables(); ++j) {
	(*__u)[(*__degreeOfFreedomSet)(j,i)] = values[__givenDegreeOfFreedomSet(j,i)];
      }
    }
  }

  K.solve(*__u, problem());

  for (size_t i=0; i<__givenDegreeOfFreedomSet.numberOfVertices(); ++i) {
    if ((*__degreeOfFreedomSet).isDOFVertex(i)) {
      for (size_t j=0; j<__givenDegreeOfFreedomSet.numberOfVariables(); ++j) {
	values[__givenDegreeOfFreedomSet(j,i)] = (*__u)[(*__degreeOfFreedomSet)(j,i)];
      }
    } else {
      for (size_t j=0; j<__givenDegreeOfFreedomSet.numberOfVariables(); ++j) {
	values[__givenDegreeOfFreedomSet(j,i)] = 0;
      }
    }
  }
}


template <typename ElementGeometry>
real_t FictitiousDomainMethod::
__integrateCharacteristicFunction(const ElementGeometry& k, const Domain& d)
{
  fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
  std::exit(1);
  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 ... " << std::flush;
  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);
  size_t numberOfUsedVertices = 0;
  ReferenceCounting<DegreeOfFreedomFDMSet::Correspondance> pCorrespondance
    = new DegreeOfFreedomFDMSet::Correspondance(mesh().numberOfVertices());

  DegreeOfFreedomFDMSet::Correspondance& correspondance = *pCorrespondance;
  correspondance = -1;

  for (size_t i=0; i<m.numberOfCells();++i) {
    typename MeshType::ElementGeometry& k = m.cell(i);
    k.setFictitious(false); // reset fictitious state!

    unsigned numberIn = 0;

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

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

    // reduction of degrees of freedom
    if (not(k.isFictitious())) {
      for (unsigned j=0; j<k.numberOfVertices(); ++j) {
	const size_t n = m.vertexNumber(k(j));
	if (correspondance[n] == -1) {
	  correspondance[n] = numberOfUsedVertices;
	  numberOfUsedVertices++;
	}
      }
    }
  }

  ConstReferenceCounting<UserFunction> u = new FEM0UserFunction(__mesh, kiValues);
  __problem = (*givenProblem) * u;
  ffout(2) << "done\n";

  __degreeOfFreedomSet = new DegreeOfFreedomFDMSet(__givenDegreeOfFreedomSet.numberOfVariables(),
						   numberOfUsedVertices,
						   pCorrespondance);
  ffout(0) << "Using "<< numberOfUsedVertices << " vertices over "
	   << mesh().numberOfVertices() << " available!\n";
}

void FictitiousDomainMethod::
computesDegreesOfFreedom(ConstReferenceCounting<Problem> givenProblem)
{
  switch (mesh().type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__computesDegreesOfFreedom<Structured3DMesh>(givenProblem);
    break;
  }
  default: {
    fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
    fferr(0) << "this mesh type is not supported by FDM penalty\n";
    std::exit(1);
  }
  }
}


template <typename MeshType>
void FictitiousDomainMethod::__discretize()
{
  MemoryManager MM;

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

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

  __u = new Vector<real_t>((*__degreeOfFreedomSet).size());

  ffout(2) << "Discretizing the PDE Problem ... ";

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

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

  (*FEM).assembleSecondMember();

  ffout(2) << "done\n";

  ffout(2) << "Discretizing Boundary Conditions\n";

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

  ffout(2) << "Second Member Modification\n";

  (*bcDiscretization).setSecondMember(__A,__b);

  ffout(2) << "Matrix Modification\n";

  (*bcDiscretization).setMatrix(__A,__b);

  ffout(2) << "done\n";

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

    SparseMatrix* aa
      = new SparseMatrix(static_cast<DoubleHashedMatrix&>(*__A));
    __A = aa; // now use sparse matrix

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

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

  switch (mesh().type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__discretize<Structured3DMesh>();
    break;
  }
  default: {
    fferr(0) << __FILE__ << ':' << __LINE__ << ": Not implemented\n";
    fferr(0) << "this mesh type is not supported by FDM penalty\n";
    std::exit(1);
  }
  }
}