//  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: DegreeOfFreedomSetBuilder.cpp,v 1.12 2007/05/20 23:34:30 delpinux Exp $

#include <DegreeOfFreedomSetBuilder.hpp>
#include <Mesh.hpp>

#include <DegreeOfFreedomSetManager.hpp>

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

#include <FiniteElementTraits.hpp>

#include <ErrorHandler.hpp>

template <size_t ItemDOFNumber>
struct DegreeOfFreedomSetBuilder::ItemTraits
{
  typedef TinyVector<ItemDOFNumber,size_t> DOFList;
};

template <>
struct DegreeOfFreedomSetBuilder::ItemTraits<0>
{
  // This Traits is built to avoid the TinyVector<0,size_t> which is undefined
  // This typedef only acts for correct syntaxe
  typedef TinyVector<1,size_t> DOFList;
};

template <typename MeshType,
	  typename FiniteElementType>
size_t DegreeOfFreedomSetBuilder::
__buildFEMCorrespondance(const MeshType& mesh,
			 DegreeOfFreedomSet::Correspondance& correspondance)
{
  size_t numberOfUsedDOF = 0;
  size_t dofPositionNumber = 0;

  const VerticesCorrespondance& verticesCorrespondance = *mesh.verticesCorrespondance();

  // vertices dof management
  if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
    for (size_t i=0; i<mesh.numberOfVertices(); ++i) {
      correspondance[i] = verticesCorrespondance[i];
      numberOfUsedDOF = std::max(verticesCorrespondance[i],numberOfUsedDOF);
      dofPositionNumber++;
    }
  }

  // edge dof management
  if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
    ASSERT(mesh.hasEdges());

    typedef TinyVector<2,size_t> VerticesPair;

    typedef
      typename ItemTraits<FiniteElementType::numberOfEdgeDegreesOfFreedom>::DOFList
      EdgeDOFNumbers;

    typedef std::map<VerticesPair, EdgeDOFNumbers> SortedEdgeSet;
    SortedEdgeSet sortedEdgeSet;

    for (size_t i=0; i<mesh.numberOfEdges(); ++i) {
      const Edge& edge = mesh.edge(i);
      const size_t v0 = verticesCorrespondance[mesh.vertexNumber(edge(0))];
      const size_t v1 = verticesCorrespondance[mesh.vertexNumber(edge(1))];

      TinyVector<2,size_t> verticesPair;
      verticesPair[0] = std::min(v0,v1);
      verticesPair[1] = std::max(v0,v1);

      typename SortedEdgeSet::iterator iedge = sortedEdgeSet.find(verticesPair);
      if (iedge == sortedEdgeSet.end()) {
	EdgeDOFNumbers dofNumbers;
	for (size_t j=0; j<FiniteElementType::numberOfEdgeDegreesOfFreedom; ++j) {
	  numberOfUsedDOF++;
	  dofNumbers[j]=numberOfUsedDOF;
	}
	
	iedge = sortedEdgeSet.insert(iedge,std::make_pair(verticesPair, dofNumbers));
      }
      for (size_t j=0; j<FiniteElementType::numberOfEdgeDegreesOfFreedom; ++j) {
	correspondance[dofPositionNumber] = iedge->second[j];
	dofPositionNumber++;
      }
    }
  }

  // faces dof management
  if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
    ASSERT(mesh.hasFaces());

    typedef TinyVector<MeshType::FaceType::NumberOfVertices,size_t> VerticesList;

    typedef
      typename ItemTraits<FiniteElementType::numberOfFaceDegreesOfFreedom>::DOFList
      FaceDOFNumbers;

    typedef std::map<VerticesList, FaceDOFNumbers> SortedFaceSet;
    SortedFaceSet sortedFaceSet;

    for (size_t i=0; i<mesh.numberOfFaces(); ++i) {
      const typename MeshType::FaceType& face = mesh.face(i);
      std::set<size_t> verticesSet;

      for (size_t j=0; j<MeshType::FaceType::NumberOfVertices; ++j) {
	verticesSet.insert(verticesCorrespondance[mesh.vertexNumber(face(j))]);
      }
      VerticesList verticesList;
      {
	size_t j=0;
	for (std::set<size_t>::const_iterator k = verticesSet.begin();
	     k != verticesSet.end(); ++k, ++j) {
	  verticesList[j] = *k;
	}
      }

      typename SortedFaceSet::iterator iface = sortedFaceSet.find(verticesList);
      if (iface == sortedFaceSet.end()) {
	FaceDOFNumbers dofNumbers;
	for (size_t j=0; j<FiniteElementType::numberOfFaceDegreesOfFreedom; ++j) {
	  numberOfUsedDOF++;
	  dofNumbers[j]=numberOfUsedDOF;
	}
	
	iface = sortedFaceSet.insert(iface,std::make_pair(verticesList, dofNumbers));
      }
      for (size_t j=0; j<FiniteElementType::numberOfFaceDegreesOfFreedom; ++j) {
	correspondance[dofPositionNumber] = iface->second[j];
	dofPositionNumber++;
      }
    }
  }

  // cell dof management
  if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
    for (size_t i=0;i<mesh.numberOfCells(); ++i) {
      for (size_t j=0;j<FiniteElementType::numberOfVolumeDegreesOfFreedom; ++j) {
	correspondance[dofPositionNumber] = numberOfUsedDOF;
	dofPositionNumber++;
	numberOfUsedDOF++;
      }
    }
  }
  return numberOfUsedDOF+1;
}

template <typename MeshType>
size_t DegreeOfFreedomSetBuilder::
__buildCorrespondance(const MeshType& mesh,
		      const DiscretizationType& discretization,
		      DegreeOfFreedomSet::Correspondance& correspondance)
{
  if (not mesh.isPeriodic()) {
    for (size_t i=0; i<correspondance.size(); ++i) {
      correspondance[i]=i;
    }
    return correspondance.size();
  }

  typedef typename MeshType::CellType CellType;

  switch (discretization.type()) {
  case DiscretizationType::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM1>::Type
      FiniteElementType;

    return this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,correspondance);
  }
  case DiscretizationType::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM2>::Type
      FiniteElementType;

    return this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,correspondance);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented",
		       ErrorHandler::unexpected);
  }
  }
  return 0;
}

template <typename MeshType,
	  typename FiniteElementType>
void DegreeOfFreedomSetBuilder::
__buildFEMFictitious(const size_t& numberOfVariables,
		     const MeshType& mesh,
		     const Domain& domain)
{
  typedef DegreeOfFreedomSet::Correspondance Correspondance;

  ReferenceCounting<Correspondance> pCorrespondance
    = new Correspondance(__dofPositionSet.number());

  size_t numberOfUsedDegreeOfFreedom=0;

  Correspondance& correspondance = *pCorrespondance;
  numberOfUsedDegreeOfFreedom
    = this->__buildFEMCorrespondance<MeshType,FiniteElementType>(mesh,
								 correspondance);

  Vector<bool> nonFictitiousDOF(__dofPositionSet.number());
  nonFictitiousDOF = false;

  for (size_t i=0; i<__dofPositionSet.number(); ++i) {
    nonFictitiousDOF[i] = domain.inside(__dofPositionSet.vertex(i));
  }

  typedef typename MeshType::CellType CellType;

  for (size_t i=0; i<mesh.numberOfCells();++i) {
    size_t numberIn = 0;
    const CellType& k = mesh.cell(i);
    const size_t cellNumer = mesh.cellNumber(k);
    // Computation of the caracteristic function
    for (size_t j=0; j<FiniteElementType::numberOfDegreesOfFreedom; ++j) {
      numberIn += (nonFictitiousDOF(__dofPositionSet(cellNumer,j)))?1:0;
    }

    // reduction of degrees of freedom
    if (numberIn != 0) {
      for (size_t j=0; j<FiniteElementType::numberOfDegreesOfFreedom; ++j) {
	const size_t n = __dofPositionSet(cellNumer,j);
	nonFictitiousDOF[n] = true;
      }
    }
  }

  // numbering correspondant vertices remapping
  std::map<size_t, size_t> newVerticesCorrespondance;
  for (size_t i=0; i<__dofPositionSet.number(); ++i) {
    if (nonFictitiousDOF[i]) {
      newVerticesCorrespondance[correspondance[i]] = 0;
    }
  }

  size_t numberOfDLVertices = 0;
  for (std::map<size_t, size_t>::iterator i=newVerticesCorrespondance.begin();
       i != newVerticesCorrespondance.end(); ++i) {
    i->second = numberOfDLVertices++;
  }

  ReferenceCounting<Correspondance> pFDMCorrespondance
    = new Correspondance(__dofPositionSet.number());

  DegreeOfFreedomSet::Correspondance& fdmCorrespondance = *pFDMCorrespondance;
  fdmCorrespondance = -1;

  for (size_t i = 0; i<__dofPositionSet.number(); ++i) {
    if (nonFictitiousDOF[i]) {
      fdmCorrespondance[i] = newVerticesCorrespondance[correspondance[i]];
    }
  }
  
  __degreeOfFreedomSet
    = new DegreeOfFreedomSet(numberOfVariables,
			     numberOfDLVertices,
			     __dofPositionSet,
			     pFDMCorrespondance);
}

template <typename MeshType>
void DegreeOfFreedomSetBuilder::
__buildFictitious(const size_t& numberOfVariables,
		  const DiscretizationType& discretization,
		  const MeshType& mesh,
		  const Domain& domain)
{
  typedef typename MeshType::CellType CellType;

  switch (discretization.type()) {
  case DiscretizationType::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM1>::Type
      FiniteElementType;

    this->__buildFEMFictitious<MeshType,FiniteElementType>(numberOfVariables,
							   mesh,
							   domain);
    break;
  }
  case DiscretizationType::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM2>::Type
      FiniteElementType;

    this->__buildFEMFictitious<MeshType,FiniteElementType>(numberOfVariables,
							   mesh,
							   domain);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented",
		       ErrorHandler::unexpected);
  }
  }
}

DegreeOfFreedomSetBuilder::
DegreeOfFreedomSetBuilder(const size_t& numberOfVariables,
			  const DiscretizationType& discretization,
			  const Mesh& mesh)
  : __dofPositionSet(DegreeOfFreedomSetManager::instance().getDOFPositionsSet(mesh,
									      discretization))
{
  typedef DegreeOfFreedomSet::Correspondance Correspondance;

  ReferenceCounting<Correspondance> pCorrespondance
    = new Correspondance(__dofPositionSet.number());
  Correspondance& correspondance = *pCorrespondance;

  size_t numberOfUsedDegreeOfFreedom=0;

  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    numberOfUsedDegreeOfFreedom
      = this->__buildCorrespondance(static_cast<const Structured3DMesh&>(mesh),
				    discretization,
				    correspondance);
    break;
  }
  case Mesh::hexahedraMesh: {
    numberOfUsedDegreeOfFreedom
      = this->__buildCorrespondance(static_cast<const MeshOfHexahedra&>(mesh),
				    discretization,
				    correspondance);
    break;
  }
  case Mesh::tetrahedraMesh: {
    numberOfUsedDegreeOfFreedom
      = this->__buildCorrespondance(static_cast<const MeshOfTetrahedra&>(mesh),
				    discretization,
				    correspondance);
    break;
  }
  case Mesh::spectralMesh: {
    numberOfUsedDegreeOfFreedom
      = this->__buildCorrespondance(static_cast<const SpectralMesh&>(mesh),
				    discretization,
				    correspondance);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented yet",
		       ErrorHandler::unexpected);
  }
  }

  __degreeOfFreedomSet
    = new DegreeOfFreedomSet(numberOfVariables,
			     numberOfUsedDegreeOfFreedom,
			     __dofPositionSet,
			     pCorrespondance);    

}

DegreeOfFreedomSetBuilder::
DegreeOfFreedomSetBuilder(const size_t& numberOfVariables,
			  const DiscretizationType& discretizationType,
			  const Mesh& mesh,
			  const Domain& domain)
  : __dofPositionSet(DegreeOfFreedomSetManager::instance().getDOFPositionsSet(mesh,
									      discretizationType))
{
  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    this->__buildFictitious(numberOfVariables,
			    discretizationType,
			    static_cast<const Structured3DMesh&>(mesh),
			    domain);
    break;
  }
  case Mesh::tetrahedraMesh: {
    this->__buildFictitious(numberOfVariables,
			    discretizationType,
			    static_cast<const Structured3DMesh&>(mesh),
			    domain);
    break;
  }
  case Mesh::hexahedraMesh: {
    this->__buildFictitious(numberOfVariables,
			    discretizationType,
			    static_cast<const Structured3DMesh&>(mesh),
			    domain);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented yet",
		       ErrorHandler::unexpected);
  }
  }
}


DegreeOfFreedomSetBuilder::
~DegreeOfFreedomSetBuilder()
{
  DegreeOfFreedomSetManager::instance().unsubscribe(__dofPositionSet);
}