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

#include "BasisFunction.hpp"
#include "MappingBasisFunction.hpp"
#include "LagrangeNodalMappingFunction.hpp"
#include "FEManager.hpp"
#include "FEManagerFace.hpp"
#include "Mesh.hpp"
#include "MeshFEDataAccessor.hpp"
#include "ReferenceSpaceFiniteElement.hpp"
#include "MonomialFiniteElementBasis.hpp"
#include "BernsteinFiniteElementBasis.hpp"
#include "JacobiFiniteElementBasis.hpp"
#include "JacobiGaussFiniteElementQuadrature.hpp"
#include "ReferenceSpaceMappingFunction.hpp"
#include "ExplicitLagrangianFiniteElementBasis.hpp"
#include "LagrangianFiniteElementBasis.hpp"
#include "DimensionBasedPolynomial.hpp"
#include "TensorProductPolynomial.hpp"

template<unsigned int dim>
class ReferenceSpaceFEManager : public FEManager<BasisFunction<dim>, MappingBasisFunction<dim,dim> >, 
                                public FEManagerFace<BasisFunction<dim>, MappingBasisFunction<dim-1,dim> > 
{
public:
  typedef BasisFunction<dim> BasisFunctionType;
  typedef MappingBasisFunction<dim,dim> MappingFunctionType;
  typedef typename BasisFunction<dim>::ShapeMatrix ShapeMatrixType;
  typedef typename BasisFunction<dim>::ShapeDerivMatrix ShapeDerivMatrixType;
  typedef ReferenceSpaceFiniteElement<ReferenceSpaceMappingFunction<ExplicitLagrangianFiniteElementBasis>,
                                      JacobiFiniteElementBasis<DimBasedPoly>,
                                      JacobiGaussFiniteElementQuadrature> FE;
  
  ReferenceSpaceFEManager() : feDataAccessor(dim)
  {
    ;
  } 
      
  ~ReferenceSpaceFEManager() { }

  void initVolumeFE(const SpaMesh::Mesh<dim>& mesh, const SpaMesh::ElementIterator& elemIt, const unsigned int p_order, const unsigned int q_order, const unsigned int addQorder = 0)
  {
    feDataAccessor.initData(mesh, elemIt);
    FE::initFE(this->basisFunction, p_order, 
	       this->mappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
	       this->quadrature, q_order, 
	       auxFeManager, feDataAccessor, addQorder);
  }

  void initBoundaryFaceFE (const SpaMesh::Mesh<dim>& mesh, const SpaMesh::ElementIterator& elemIt, const int side_number, 
                           const unsigned int p_order, const unsigned int q_order)
  {
    feDataAccessor.initData(mesh, elemIt, side_number, true);
    FE::initFE(this->boundaryBasisFunction, p_order, 
               this->boundaryMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, q_order, 
	       auxFeManager, feDataAccessor);
  }
  
  void initBoundaryFaceFE (const SpaMesh::Mesh<dim>& mesh, const int& elem, const SpaMesh::FaceIterator& sideIt, const unsigned int p_order, const unsigned int q_order)
  {
    feDataAccessor.initData(mesh, elem, sideIt.value(), true);
    FE::initFE(this->boundaryBasisFunction, p_order, 
               this->boundaryMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, q_order, 
	       auxFeManager, feDataAccessor);
  }
  
  void initInternalFaceFE (const SpaMesh::Mesh<dim>& mesh, 
                           const SpaMesh::ElementIterator& elemIt, const SpaMesh::NeighborIterator& neighborIt, const unsigned int p_order, const unsigned int q_order)
  {
    feDataAccessor.initData(mesh, elemIt, neighborIt, true);
    FE::initFE(this->internalElemBasisFunction, p_order, 
               this->internalNeighborBasisFunction, p_order, 
    	       this->internalMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, q_order, 
	       auxFeManager, feDataAccessor);
  }
  void initInternalFaceFE (const SpaMesh::Mesh<dim>& mesh, 
                           const int& elem, const int& neighbor, 
			   const SpaMesh::FaceIterator& sideIt, const unsigned int p_order, const unsigned int q_order)
  {
    feDataAccessor.initData(mesh, elem, neighbor, sideIt, true);
    FE::initFE(this->internalElemBasisFunction, p_order, 
               this->internalNeighborBasisFunction, p_order, 
    	       this->internalMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, q_order, 
	       auxFeManager, feDataAccessor);
  }

  double computeFESize(const SpaMesh::Mesh<dim>& mesh, const int& elem)
  {
    feDataAccessor.initData(mesh, elem);
    return FE::sizeFE(feDataAccessor, sizeMapping, sizeQuadrature);
  }
  double computeFESize(const SpaMesh::Mesh<dim>& mesh, SpaMesh::ElementIterator& elemIt)
  {
    feDataAccessor.initData(mesh, elemIt);
    return FE::sizeFE(feDataAccessor, sizeMapping, sizeQuadrature);
  }
  double computeFESize(const SpaMesh::Mesh<dim>& mesh, SpaMesh::NeighborIterator& elemIt)
  {
    feDataAccessor.initData(mesh, elemIt);
    return FE::sizeFE(feDataAccessor, sizeMapping, sizeQuadrature);
  }
  
  MeshFEDataAccessor& getFEDataAccessor()
  {
    return feDataAccessor;
  }
  
private:
  MeshFEDataAccessor feDataAccessor;
  MappingFunctionType sizeMapping;
  FEManager<BasisFunctionType, MappingFunctionType> auxFeManager;
  Quadrature sizeQuadrature;
};

#endif
