//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001-2005 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: DegreeOfFreedomPositionsSet.cpp,v 1.9 2007/06/17 20:46:10 delpinux Exp $

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

#include <FiniteElementTraits.hpp>

#include <DegreeOfFreedomPositionsSet.hpp>
#include <ConnectivityBuilder.hpp>

/**
 * @file   DegreeOfFreedomPositionsSet.cpp
 * @author Stephane Del Pino
 * @date   Sat Nov 11 17:24:59 2006
 * 
 * @brief  Builds the DegreeOfFreedomPositionsSet
 */
class DegreeOfFreedomPositionsSet::Builder
{
private:
  DegreeOfFreedomPositionsSet&
  __dofPositionsSet;		/**< Reference to the
				   DegreeOfFreedomPositionsSet */

  /** 
   * Computes the number of degrees of freedom on a mesh for a given
   * finite element type
   * 
   * @param mesh given mesh
   * 
   * @return the number of degrees of freedom
   */
  template <typename FiniteElementType,
	    typename MeshType>
  size_t __computeFEMNbDOF(const MeshType& mesh)
  {
    size_t nbDOF = 0;
    size_t neededConnectivities = 0;

    // adds vertices dof
    if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
      nbDOF
	+= FiniteElementType::numberOfVertexDegreesOfFreedom
	*  mesh.numberOfVertices();
    }
    
    // adds edges dof
    if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
      if (not(mesh.hasEdges())) {
	const_cast<MeshType&>(mesh).buildEdges();
      }
      if (not(mesh.connectivity().hasCellToEdges())) {
	neededConnectivities += Connectivity<MeshType>::CellToEdges;
      }

      nbDOF
	+= FiniteElementType::numberOfEdgeDegreesOfFreedom
	*  mesh.numberOfEdges();
    }

    // adds faces dof
    if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
      if (not(mesh.hasFaces())) {
	const_cast<MeshType&>(mesh).buildFaces();
      }
      if (not(mesh.connectivity().hasCellToFaces())) {
	neededConnectivities += Connectivity<MeshType>::CellToFaces;
      }
      nbDOF
	+= FiniteElementType::numberOfFaceDegreesOfFreedom
	*  mesh.numberOfFaces();
    }

    // adds cells dof
    if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
      nbDOF
	+= FiniteElementType::numberOfVolumeDegreesOfFreedom
	*  mesh.numberOfCells();
    }

    if (neededConnectivities > 0) {
      ConnectivityBuilder<MeshType> builder(mesh);
      builder.generates(neededConnectivities);
    }

    ffout(3) << "- number of degrees of freedom positions: " << nbDOF << '\n';
    return nbDOF;
  }

  /** 
   * Builds the degree of freedom set for a finite element type and a
   * given mesh
   * 
   * @param mesh given mesh
   */
  template<typename MeshType,
	   typename FiniteElementType>
  void __buildFEM(const MeshType& mesh)
  {
    typedef typename MeshType::CellType CellType;

    __dofPositionsSet.__numberOfDOFPerCell = FiniteElementType::numberOfDegreesOfFreedom;
    __dofPositionsSet.__dofNumber.resize(mesh.numberOfCells()*FiniteElementType::numberOfDegreesOfFreedom);
    // fills the correspondance table with wrong values
    __dofPositionsSet.__dofNumber = std::numeric_limits<size_t>::max();

    // Computes the number of degrees of freedom
    const size_t nbDOF = __computeFEMNbDOF<FiniteElementType>(mesh);

    __dofPositionsSet.__positions.resize(nbDOF);

    size_t begin = 0;
    size_t end = 0;

    // vertices
    if (FiniteElementType::numberOfVertexDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfVertices();
      for (size_t i=begin; i < end; ++i) {
	__dofPositionsSet.__positions[i] = mesh.vertex(i);
      }
      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
	const CellType& cell = mesh.cell(i);
	for (size_t j=0; j<CellType::NumberOfVertices; ++j) {
	  const size_t vertexNumber = mesh.vertexNumber(cell(j));
	  __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+j]
	    = vertexNumber;
	}
      }
    }

    // edges
    if (FiniteElementType::numberOfEdgeDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfEdges();
      for (size_t i=begin; i < end; ++i) {
	const Edge& edge = mesh.edge(i-begin);
	__dofPositionsSet.__positions[i] = 0.5*(edge(0)+edge(1));
      }

      // Take into account vertices degrees of freedom in numbering
      const size_t inputShift
	= FiniteElementType::numberOfVertexDegreesOfFreedom
	* CellType::NumberOfVertices;
      const size_t outputShift
	= mesh.numberOfVertices()
	* FiniteElementType::numberOfVertexDegreesOfFreedom;

      const Connectivity<MeshType>&
	connectivity = mesh.connectivity();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
	const CellType& cell = mesh.cell(i);
	const typename Connectivity<MeshType>::CellToEdgesType&
	  cellEdges = connectivity.edges(cell);
	for (size_t j=0; j<CellType::NumberOfEdges; ++j) {
	  const Edge& edge = *(cellEdges[j]);
	  const size_t edgeNumber = mesh.edgeNumber(edge);
	  __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift+j]
	    = edgeNumber + outputShift;
	}
      }
    }

    // faces
    if (FiniteElementType::numberOfFaceDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfFaces();
      for (size_t i=begin; i < end; ++i) {
	typedef typename MeshType::CellType::FaceType FaceType;

	const FaceType& face = mesh.face(i-begin);
	__dofPositionsSet.__positions[i] = 0;
	for (size_t j=0; j<FaceType::NumberOfVertices; ++j) {
	  __dofPositionsSet.__positions[i] += face(j);
	}
	__dofPositionsSet.__positions[i] *= 1./FaceType::NumberOfVertices;
      }

      // Take into account vertices and edges degrees of freedom in
      // numbering
      const size_t inputShift
	= FiniteElementType::numberOfVertexDegreesOfFreedom
	* CellType::NumberOfVertices
	+ FiniteElementType::numberOfEdgeDegreesOfFreedom
	* CellType::NumberOfEdges;
      const size_t outputShift
	= mesh.numberOfVertices()
	* FiniteElementType::numberOfVertexDegreesOfFreedom
	+ (mesh.hasEdges() ? mesh.numberOfEdges() : 0)
	* FiniteElementType::numberOfEdgeDegreesOfFreedom;

      const Connectivity<MeshType>&
	connectivity = mesh.connectivity();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
	const CellType& cell = mesh.cell(i);
	const typename Connectivity<MeshType>::CellToFacesType&
	  cellFaces = connectivity.faces(cell);

	for (size_t j=0; j<CellType::NumberOfFaces; ++j) {
	  const typename CellType::FaceType& face = *(cellFaces[j]);
	  const size_t faceNumber = mesh.faceNumber(face);
	  __dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift+j]
	    = faceNumber + outputShift;
	}
      }
    }

    // cells
    if (FiniteElementType::numberOfVolumeDegreesOfFreedom > 0) {
      begin = end;
      end += mesh.numberOfCells();

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
	const CellType& cell = mesh.cell(i);
	TinyVector<3, real_t> position = 0;
	for (size_t l = 0; l<MeshType::CellType::NumberOfVertices; ++l) {
	  position += cell(l);
	}
	position /= MeshType::CellType::NumberOfVertices;
	__dofPositionsSet.__positions[i+begin] = position;
      }

      // Take into account vertices, edges and face degrees of freedom
      // in numbering
      const size_t inputShift
	= FiniteElementType::numberOfVertexDegreesOfFreedom
	* CellType::NumberOfVertices
	+ FiniteElementType::numberOfEdgeDegreesOfFreedom
	* CellType::NumberOfEdges
	+ FiniteElementType::numberOfFaceDegreesOfFreedom
	* CellType::NumberOfFaces;
      const size_t outputShift
	= mesh.numberOfVertices()
	* FiniteElementType::numberOfVertexDegreesOfFreedom
	+ (mesh.hasEdges() ? mesh.numberOfEdges() : 0)
	* FiniteElementType::numberOfEdgeDegreesOfFreedom
	+ (mesh.hasFaces() ? mesh.numberOfFaces() : 0)
	* FiniteElementType::numberOfFaceDegreesOfFreedom;

      for (size_t i=0; i<mesh.numberOfCells(); ++i) {
	__dofPositionsSet.__dofNumber[i*FiniteElementType::numberOfDegreesOfFreedom+inputShift]
	  = i + outputShift;
      }
    }
  }

  template<typename MeshType>
  void __buildLegendre(const MeshType& mesh)
  {
    throw ErrorHandler(__FILE__,__LINE__,
		       "Cannot build spectral degrees of freefom on '"
		       +mesh.typeName()+"'",
		       ErrorHandler::unexpected);
  }

  /** 
   * Builds the degrees of freedom position set for a discretization
   * type and a given mesh
   * 
   * @param discretizationType discretization type
   * @param mesh given mesh
   */
  template <typename MeshType>
  void __build(const DiscretizationType& discretizationType,
	       const MeshType& mesh);
public:
  /** 
   * Constructor
   * 
   * @param dofPositionsSet degrees of freedom position set to built
   */
  Builder(DegreeOfFreedomPositionsSet& dofPositionsSet)
    : __dofPositionsSet(dofPositionsSet)
  {
    ;
  }

  /** 
   * Destructor
   * 
   */
  ~Builder()
  {
    ;
  }

  /** 
   * Build a DegreeOfFreedomPositionsSet
   * 
   * @param discretizationType given discretizationType
   * @param mesh given mesh
   */
  void build(const DiscretizationType& discretizationType,
	     const Mesh& mesh);
};


template<>
void DegreeOfFreedomPositionsSet::Builder::
__buildLegendre<SpectralMesh>(const SpectralMesh& mesh)
{
  // All the degrees of freedom are contained in 1 spectral cell...
  __dofPositionsSet.__numberOfDOFPerCell = (mesh.degree(0)+1)*(mesh.degree(1)+1)*(mesh.degree(2)+1);
  __dofPositionsSet.__dofNumber.resize(__dofPositionsSet.__numberOfDOFPerCell);

  for (size_t i=0; i<__dofPositionsSet.__dofNumber.size(); ++i) {
    __dofPositionsSet.__dofNumber[i]=i;
  }

  // legendre degrees of freedom do not really have a position...
  __dofPositionsSet.__positions.resize(__dofPositionsSet.__numberOfDOFPerCell);
  __dofPositionsSet.__positions = TinyVector<3,real_t>(0,0,0);
}


template <typename MeshType>
void DegreeOfFreedomPositionsSet::Builder::
__build(const DiscretizationType& discretizationType,
	const MeshType& mesh)
{
  typedef typename MeshType::CellType CellType;
  switch(discretizationType.type()) {
  case DiscretizationType::lagrangianFEM0: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM0>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case DiscretizationType::lagrangianFEM1: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM1>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case DiscretizationType::lagrangianFEM2: {
    typedef
      typename FiniteElementTraits<CellType,
      DiscretizationType::lagrangianFEM2>::Type
      FiniteElementType;
    __buildFEM<MeshType, FiniteElementType>(mesh);
    break;
  }
  case DiscretizationType::spectralLegendre: {
    __buildLegendre<MeshType>(mesh);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented",
		       ErrorHandler::unexpected);
  }
  }
}


void DegreeOfFreedomPositionsSet::Builder::
build(const DiscretizationType& discretizationType,
      const Mesh& mesh)
{
  ffout(2) << "Computing degrees of freedom positions...\n";
  switch(mesh.type()) {
  case Mesh::cartesianHexahedraMesh: {
    __build(discretizationType, static_cast<const Structured3DMesh&>(mesh));
    break;
  }
  case Mesh::hexahedraMesh: {
    __build(discretizationType, static_cast<const MeshOfHexahedra&>(mesh));
    break;
  }
  case Mesh::spectralMesh: {
    __build(discretizationType, static_cast<const SpectralMesh&>(mesh));
    break;
  }
  case Mesh::tetrahedraMesh: {
    __build(discretizationType, static_cast<const MeshOfTetrahedra&>(mesh));
    break;
  }
  case Mesh::trianglesMesh: {
    __build(discretizationType, static_cast<const MeshOfTriangles&>(mesh));
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "not implemented",
		       ErrorHandler::unexpected);
  }
  }
  ffout(2) << "Computing degrees of freedom positions: done\n";
}

void DegreeOfFreedomPositionsSet::
__build(const DiscretizationType& discretizationType,
	const Mesh& mesh)
{
  Builder builder(*this);

  builder.build(discretizationType, mesh);
}