// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */ 
/* License as published by the Free Software Foundation either */ 
/* version 3.0 of the License, or (at your option) any later version. */

/* This software 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 */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/. */

#ifndef __reference_space_finite_element_h__
#define __reference_space_finite_element_h__

#include "BasisFunction.hpp"
#include "FiniteElementBasis.hpp"
#include "FiniteElementMapping.hpp"
#include "MappingBasisFunction.hpp"
#include "LagrangeNodalMappingFunction.hpp"
#include "CollectiveMappingFunction.hpp"
#include "Quadrature.hpp"
#include "FiniteElementQuadrature.hpp"
#include "ReferenceFiniteElement.hpp"
#include "Mesh.hpp"
#include "MeshFEDataAccessor.hpp"
#include "FEManager.hpp"

template<class MappingT, class BasisT, class QuadratureT>
class ReferenceSpaceFiniteElement 
{
public:
  typedef FiniteElementMapping<MappingT> FiniteElementMappingType;
  typedef FiniteElementBasis<BasisT> FiniteElementBasisType;
  typedef FiniteElementQuadrature<QuadratureT> FiniteElementQuadratureType;
  typedef SpaMesh::Connectivity::TopoPointVector TopoPointVector;
  typedef SpaMesh::Connectivity::TopoIndexVector TopoIndexVector;

  ReferenceSpaceFiniteElement() 
    {
      ;
    }
  
  ~ReferenceSpaceFiniteElement() {}
  
public:
  template <unsigned int Dim, unsigned int DIm, 
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void initFE(BasisFunction<DIm>& elemBasisFunction, const unsigned int elemBorder,
                       BasisFunction<DIm>& neighborBasisFunction, const unsigned int neighborBorder,
                       MappingBasisFunctionType<Dim,DIm>& sideMappingFunction, const unsigned int m_order,
                       Quadrature& quadrature, const unsigned int q_order,
		       FEManager<BasisFunction<DIm>, MappingBasisFunctionType<DIm,DIm> >& auxFEManager,
		       const MeshFEDataAccessor& feTOmoab)
  {
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    const TopoPointVector& neighbor_coord = feTOmoab.getNeighborCoord();
    const TopoPointVector& side_coord = feTOmoab.getSideCoord();

    const TopoIndexVector elem_side_indices = feTOmoab.getElemSideIndices();
    const TopoIndexVector neighbor_side_indices = feTOmoab.getNeighborSideIndices();

    const ElemType elemType = feTOmoab.getElemType();
    const ElemType neighborType = feTOmoab.getNeighborType();
    const ElemType sideType = feTOmoab.getSideType();

    FiniteElementQuadratureType::init(quadrature, sideType, q_order);
    FiniteElementMappingType::init(sideMappingFunction, quadrature, sideType, side_coord, m_order);
    typename Quadrature::PointVector& _auxQuadrature = auxFEManager.quadrature.getQpoint();
    _auxQuadrature.resize(quadrature.nbPoint(), Eigen::NoChange);
    
    // map to find side quadrature
    FiniteElementMappingType::referenceSpaceSideMap(neighborType, sideType, neighbor_side_indices,
					       sideMappingFunction, _auxQuadrature);
    FiniteElementBasisType::reInit(neighborBasisFunction, _auxQuadrature, neighborType, neighborBorder);
    if (neighborBasisFunction.computeShapeDeriv())
    {
      computeSideBasisFunction(neighborBasisFunction, neighborType, neighbor_coord, m_order, auxFEManager);
    }
    
    FiniteElementMappingType::referenceSpaceSideMap(elemType, sideType, elem_side_indices, 
					       sideMappingFunction, _auxQuadrature);
    FiniteElementBasisType::reInit(elemBasisFunction, _auxQuadrature, elemType, elemBorder);
    if (elemBasisFunction.computeShapeDeriv())
    {
      computeSideBasisFunction(elemBasisFunction, elemType, elem_coord, m_order, auxFEManager);
    }

  }
  
  template <unsigned int Dim, unsigned int DIm,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void initFE(BasisFunction<DIm>& elemBasisFunction, const unsigned int elemBorder,
                       MappingBasisFunctionType<Dim,DIm>& sideMappingFunction, const unsigned int m_order,
                       Quadrature& quadrature, const unsigned int q_order,
		       FEManager<BasisFunction<DIm>, MappingBasisFunctionType<DIm,DIm> >& auxFEManager,
		       const MeshFEDataAccessor& feTOmoab)
  {
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    const TopoPointVector& side_coord = feTOmoab.getSideCoord();

    const TopoIndexVector& elem_side_indices = feTOmoab.getElemSideIndices();

    const ElemType elemType = feTOmoab.getElemType();
    const ElemType sideType = feTOmoab.getSideType();

    FiniteElementQuadratureType::init(quadrature, sideType, q_order);
    FiniteElementMappingType::init(sideMappingFunction, quadrature, sideType, side_coord, m_order);
    typename Quadrature::PointVector& _auxQuadrature = auxFEManager.quadrature.getQpoint();
    _auxQuadrature.resize(quadrature.nbPoint(), Eigen::NoChange);

    // map to find side quadrature
    FiniteElementMappingType::referenceSpaceSideMap(elemType, sideType, elem_side_indices, 
					       sideMappingFunction, _auxQuadrature);
    FiniteElementBasisType::reInit(elemBasisFunction, _auxQuadrature, elemType, elemBorder);
    if (elemBasisFunction.computeShapeDeriv())
    {
      computeSideBasisFunction(elemBasisFunction, elemType, elem_coord, m_order, auxFEManager);
    }
  }
  
  template <unsigned int Dim,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void initFE(BasisFunction<Dim>& basisFunction, const unsigned int b_order,
                       MappingBasisFunctionType<Dim,Dim>& mappingFunction, const unsigned int m_order,
                       Quadrature& quadrature, const unsigned int q_order,
	               FEManager<BasisFunction<Dim>, MappingBasisFunctionType<Dim,Dim> >& auxFEManager,
		       const MeshFEDataAccessor& feTOmoab, const unsigned int additionalQorder = 0)
  {
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    const ElemType elemType = feTOmoab.getElemType();

    FiniteElementQuadratureType::init(quadrature, elemType, q_order+additionalQorder);
    FiniteElementMappingType::init(mappingFunction, quadrature, elemType, elem_coord, m_order);
    computeBasisFunction(basisFunction, b_order,  mappingFunction, quadrature, elemType, auxFEManager);
  }
  
  template <unsigned int Dim,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void initFE(BasisFunction<Dim>& basisFunction, const unsigned int b_order,
                       MappingBasisFunctionType<Dim,Dim>& mappingFunction, const unsigned int m_order,
                       Quadrature& quadrature, const unsigned int q_order,
		       const MeshFEDataAccessor& feTOmoab, const unsigned int additionalQorder = 0)
  {
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    const ElemType elemType = feTOmoab.getElemType();

    FiniteElementQuadratureType::init(quadrature, elemType, q_order+additionalQorder);
    FiniteElementMappingType::init(mappingFunction, quadrature, elemType, elem_coord, m_order);
    
    FEManager<BasisFunction<Dim>, MappingBasisFunctionType<Dim,Dim> > auxFEManager;
    computeBasisFunction(basisFunction, b_order,  mappingFunction, quadrature, elemType, auxFEManager);
  }
  
  template <unsigned int Dim, 
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static double sizeFE(MeshFEDataAccessor& feTOmoab, 
	                 MappingBasisFunctionType<Dim,Dim>& sizeMapping, 
			 Quadrature& sizeQuadrature)
  {
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    const ElemType elemType = feTOmoab.getElemType();

    const unsigned int q_order = ReferenceFiniteElement::getGeoOrder(feTOmoab.getElemType());
    FiniteElementQuadratureType::init(sizeQuadrature, elemType, q_order);
    FiniteElementMappingType::init(sizeMapping, sizeQuadrature, elemType, elem_coord, q_order);
    return sizeMapping.getJacobianXweight().sum();
  }
  
  template <unsigned int Dim,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void initFENodalPoints(BasisFunction<Dim>& basisFunction, const unsigned int b_order,
	                          MappingBasisFunctionType<Dim,Dim>& mappingFunction, const unsigned int m_order,
                                  const MeshFEDataAccessor& feTOmoab,
			          const Quadrature::PointVector& referenceSpaceNodes)
  {
    const ElemType elemType = feTOmoab.getElemType();
    const TopoPointVector& elem_coord = feTOmoab.getElemCoord();
    if (!basisFunction.computeShapeDeriv())
    {
      FiniteElementMappingType::computePhysicalSpacePoints(mappingFunction,referenceSpaceNodes,elemType,elem_coord,m_order);
      FiniteElementBasisType::init(basisFunction, referenceSpaceNodes, elemType, b_order);
    }
    else
    {
      FiniteElementMappingType::computePhysicalSpacePointsWithMap(mappingFunction,referenceSpaceNodes,elemType,elem_coord,m_order);
      BasisFunction<Dim> auxBasisFunction;
      FiniteElementBasisType::init(auxBasisFunction, referenceSpaceNodes, elemType, b_order);
      basisFunction.getShapeFunction() = auxBasisFunction.getShapeFunction();
      typename BasisFunction<Dim>::ShapeDerivMatrix& shapeDerivMatrix= basisFunction.getShapeFunctionDeriv();
      typename BasisFunction<Dim>::ShapeDerivMatrix& auxShapeDerivMatrix= auxBasisFunction.getShapeFunctionDeriv();
      const typename MappingBasisFunction<Dim,Dim>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
      for (unsigned int k = 0; k < Dim; k++)
      {
        shapeDerivMatrix(k).resize(auxShapeDerivMatrix(k).rows(),auxShapeDerivMatrix(k).cols());
        shapeDerivMatrix(k).setZero();
        for (unsigned int j = 0; j < Dim; j++)
          for (unsigned int qp = 0; qp < referenceSpaceNodes.rows(); qp++)
            shapeDerivMatrix(k).col(qp).noalias() += auxShapeDerivMatrix(j).col(qp) * dEdX(qp)(j,k);
      }
    }
  }
  
  static bool is_hierarchical ()
  {
    return FiniteElementBasisType::is_hierarchical();
  }

private:
  template <unsigned int Dim,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void computeBasisFunction(BasisFunction<Dim>& basisFunction, const unsigned int b_order, 
                                     const MappingBasisFunctionType<Dim,Dim>& mappingFunction, 
			             const Quadrature& quadrature, const ElemType& elemType,
			             FEManager<BasisFunction<Dim>, MappingBasisFunctionType<Dim, Dim> >& auxFEManager)
  {
    FiniteElementBasisType::init(auxFEManager.basisFunction, quadrature, elemType, b_order, basisFunction.getShapeFunction());
    typename BasisFunction<Dim>::ShapeDerivMatrix& shapeDerivMatrix= basisFunction.getShapeFunctionDeriv();
    typename BasisFunction<Dim>::ShapeDerivMatrix& auxShapeDerivMatrix= auxFEManager.basisFunction.getShapeFunctionDeriv();
    const typename MappingBasisFunction<Dim,Dim>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int k = 0; k < Dim; k++)
    {
      shapeDerivMatrix(k).resize(auxShapeDerivMatrix(k).rows(),auxShapeDerivMatrix(k).cols());
      shapeDerivMatrix(k).setZero();
      for (unsigned int j = 0; j < Dim; j++)
        for (unsigned int qp = 0; qp < quadrature.nbPoint(); qp++)
          shapeDerivMatrix(k).col(qp).noalias() += auxShapeDerivMatrix(j).col(qp) * dEdX(qp)(j,k);
    }
  }
  
  template <unsigned int Dim,
	    template <unsigned int, unsigned int> class MappingBasisFunctionType>
    static void computeSideBasisFunction(BasisFunction<Dim>& basisFunction, 
                                         const ElemType& elemType, 
				         const TopoPointVector& elem_coord, const unsigned int m_order,
			                 FEManager<BasisFunction<Dim>, MappingBasisFunctionType<Dim,Dim> >& auxFEManager)
  {
    FiniteElementMappingType::reInitSide(auxFEManager.mappingFunction, auxFEManager.quadrature.getQpoint(), elemType, elem_coord, m_order);
    typename BasisFunction<Dim>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
    const typename MappingBasisFunction<Dim,Dim>::MappingMatrixInv& dEdX = auxFEManager.mappingFunction.getdEdX();
    Eigen::Matrix<Eigen::Matrix<double, Eigen::Dynamic, 1>, Dim, 1> ishapeDerivBkp;
    for (unsigned int qp = 0; qp < shapeDerivMatrix(0).cols(); qp++)
    {
      for (unsigned int k = 0; k < Dim; k++)
      {
        ishapeDerivBkp(k) = shapeDerivMatrix(k).col(qp);
        shapeDerivMatrix(k).col(qp).setZero();
      }
      for (unsigned int k = 0; k < Dim; k++)
        for (unsigned int j = 0; j < Dim; j++)
          shapeDerivMatrix(k).col(qp).noalias() += ishapeDerivBkp(j) * dEdX(qp)(j,k);
    }
  }
};

#endif
