// 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 __HHOPhysicaSpaceFEManager_h__
#define __HHOPhysicaSpaceFEManager_h__

#include "BasisFunction.hpp"
#include "MappingBasisFunction.hpp"
#include "ItemMappingFunction.hpp"
#include "FEManager.hpp"
#include "FEManagerFaceHHO.hpp"
#include "Mesh.hpp"
#include "MeshFEDataAccessor.hpp"
#include "PhysicalSpaceFiniteElement.hpp"
#include "PhysicalSpaceScaledRotatedFiniteElement.hpp"
#include "MonomialFiniteElementBasis.hpp"
#include "JacobiGaussFiniteElementQuadrature.hpp"
#include "ReferenceSpaceMappingFunction.hpp"
#include "ExplicitLagrangianFiniteElementBasis.hpp"
#include "DimensionBasedPolynomial.hpp"

template<unsigned int dim>
class HHOPhysicalSpaceFEManager : public FEManager<BasisFunction<dim>, ItemMappingFunction<MappingBasisFunction<dim,dim> > >, 
                                  public FEManagerFaceHHO<BasisFunction<dim-1>, BasisFunction<dim>, ItemMappingFunction<MappingBasisFunction<dim-1,dim> > >
{
public:
  typedef BasisFunction<dim> BasisFunctionType;
  typedef BasisFunction<dim-1> FaceBasisFunctionType;
  typedef ItemMappingFunction<MappingBasisFunction<dim,dim> > MappingFunctionType;
  typedef ItemMappingFunction<MappingBasisFunction<dim-1,dim> > FaceMappingFunctionType;
  typedef typename BasisFunctionType::ShapeMatrix ShapeMatrixType;
  typedef typename BasisFunctionType::ShapeDerivMatrix ShapeDerivMatrixType;
  typedef typename FaceBasisFunctionType::ShapeMatrix FaceShapeMatrixType;
  typedef typename FaceBasisFunctionType::ShapeDerivMatrix FaceShapeDerivMatrixType;
  typedef PhysicalSpaceScaledRotatedFiniteElement<ReferenceSpaceMappingFunction<ExplicitLagrangianFiniteElementBasis>,
                                                  MonomialFiniteElementBasis<DimBasedPoly>,
                                                  JacobiGaussFiniteElementQuadrature> FE;
  
  HHOPhysicalSpaceFEManager() : feDataAccessor(dim)
  {
    ;
  } 
      
  ~HHOPhysicalSpaceFEManager() { }

  void initVolumeFE(const SpaMesh::Mesh<dim>& mesh, const SpaMesh::ElementIterator& elemIt, const unsigned int p_order)
  {
    feDataAccessor.initData(mesh, elemIt);
    const unsigned int ppo_order = p_order + 1;
    FE::initFE(this->basisFunction, ppo_order, 
               this->mappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
               this->quadrature, ppo_order*2+1,
               elemGeoInfo, feDataAccessor);
  }
  
  void initFaceElemFaceFE (const SpaMesh::Mesh<dim>& mesh, const SpaMesh::ElementIterator& elemIt, const int& side, const unsigned int p_order)
  {
    feDataAccessor.initData(mesh, elemIt.value(), side, false, true);
    const unsigned int ppo_order = p_order + 1;
    FE::initFE(this->elementFaceBasisFunction, ppo_order, 
               this->faceMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, ppo_order*2+1, 
               elemGeoInfo, feDataAccessor);
    FE::initFE(this->faceBasisFunction, p_order, 
               this->faceMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, ppo_order*2+1, 
               faceGeoInfo, feDataAccessor);
  }
  
  void initFaceFE (const SpaMesh::Mesh<dim>& mesh, const SpaMesh::ElementIterator& elemIt, const int& side, const unsigned int p_order)
  {
    feDataAccessor.initData(mesh, elemIt.value(), side, false, true);
    const unsigned int ppo_order = p_order + 1;
    FE::initFE(this->faceBasisFunction, p_order, 
               this->faceMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, ppo_order*2+1, 
               faceGeoInfo, feDataAccessor);
  }
  
  double computeFESize(const SpaMesh::Mesh<dim>& mesh, const int& elem)
  {
    feDataAccessor.initData(mesh, elem);
    return FE::sizeFE(feDataAccessor, elemGeoInfo.getMapping(), elemGeoInfo.getQuadrature());
  }
  double computeFESize(const SpaMesh::Mesh<dim>& mesh, SpaMesh::ElementIterator& elemIt)
  {
    feDataAccessor.initData(mesh, elemIt);
    return FE::sizeFE(feDataAccessor, elemGeoInfo.getMapping(), elemGeoInfo.getQuadrature());
  }
  double computeFESize(const SpaMesh::Mesh<dim>& mesh, SpaMesh::NeighborIterator& elemIt)
  {
    feDataAccessor.initData(mesh, elemIt);
    return FE::sizeFE(feDataAccessor, elemGeoInfo.getMapping(), elemGeoInfo.getQuadrature());
  }
  
  double computeFEFaceSize(const SpaMesh::Mesh<dim>& mesh, const int& side)
  {
    feDataAccessor.initDataSide(mesh, side);
    return FE::sizeFE(feDataAccessor, faceGeoInfo.getMapping(), faceGeoInfo.getQuadrature());
  }
  
  MeshFEDataAccessor& getFEDataAccessor()
  {
    return feDataAccessor;
  }

private:
  MeshFEDataAccessor feDataAccessor;
  ElemGeometryInfo<MappingFunctionType> elemGeoInfo;
  ElemGeometryInfo<FaceMappingFunctionType> faceGeoInfo;
};

#endif
