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

#include "BasisFunction.hpp"
#include "MappingBasisFunction.hpp"
#include "ItemMappingFunction.hpp"
#include "FEManager.hpp"
#include "FEManagerFace.hpp"
#include "MeshFEDataAccessor.hpp"
#include "PhysicalSpaceScaledRotatedFiniteElement.hpp"
#include "PhysicalSpaceFiniteElement.hpp"
#include "OrthoNormalPhysicalSpaceFiniteElement.hpp"
#include "MonomialFiniteElementBasis.hpp"
#include "JacobiGaussFiniteElementQuadrature.hpp"
#include "ReferenceSpaceMappingFunction.hpp"
#include "ExplicitLagrangianFiniteElementBasis.hpp"
#include "DimensionBasedPolynomial.hpp"
#include "LiftingOperator.hpp"
#include "ElementIterator.hpp"

template<unsigned int dim>
class OrthoNormalPhysicalSpaceFEManager : public FEManager<BasisFunction<dim>, ItemMappingFunction<MappingBasisFunction<dim,dim> > >, 
                                          public FEManagerFace<BasisFunction<dim>, ItemMappingFunction<MappingBasisFunction<dim-1,dim> > >
{
public:
  typedef BasisFunction<dim> BasisFunctionType;
  typedef ItemMappingFunction<MappingBasisFunction<dim,dim> > MappingFunctionType;
  typedef typename BasisFunctionType::ShapeMatrix ShapeMatrixType;
  typedef typename BasisFunctionType::ShapeDerivMatrix ShapeDerivMatrixType;
  typedef OrthoNormalPhysicalSpaceFiniteElement<PhysicalSpaceScaledRotatedFiniteElement<ReferenceSpaceMappingFunction<ExplicitLagrangianFiniteElementBasis>,
                                                                                        MonomialFiniteElementBasis<DimBasedPoly>,
                                                                                        JacobiGaussFiniteElementQuadrature> > FE;
  typedef FE::OrthoMatrixType OrthoMatrixType;

  OrthoNormalPhysicalSpaceFEManager() : feDataAccessor(dim)
  {
    ;
  } 
      
  ~OrthoNormalPhysicalSpaceFEManager() { }

  const typename LiftingOperator<dim>::LiftingMatrix& getLifting11() const 
  {
    return lifting.r11(); 
  }
  const typename LiftingOperator<dim>::LiftingMatrix& getLifting12() const 
  {
    return lifting.r12(); 
  }
  const typename LiftingOperator<dim>::LiftingMatrix& getLifting22() const 
  {
    return lifting.r22(); 
  }
  const typename LiftingOperator<dim>::LiftingMatrix& getLifting21() const 
  {
    return lifting.r21(); 
  }

  void computeInternalFaceLifting()
  {
    assert(this->sideQuadrature.nbPoint() == this->internalMappingFunction.getQpointX().rows());
    assert(this->sideQuadrature.nbPoint() == this->internalElemBasisFunction.nbQuadraturePoint());
    assert(this->sideQuadrature.nbPoint() == this->internalNeighborBasisFunction.nbQuadraturePoint());
    lifting.compute(this->internalElemBasisFunction, this->internalNeighborBasisFunction, this->internalMappingFunction, 0, this->sideQuadrature.nbPoint());
  }
  
  void computeBoundaryFaceLifting()
  {
    assert(this->sideQuadrature.nbPoint() == this->boundaryMappingFunction.getQpointX().rows());
    assert(this->sideQuadrature.nbPoint() == this->boundaryBasisFunction.nbQuadraturePoint());
    lifting.compute(this->boundaryBasisFunction, this->boundaryMappingFunction, 0, this->sideQuadrature.nbPoint());
  }
  
  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);
    OrthoMatrixType ortCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, ortCoeff);
    FE::initFE(this->basisFunction, p_order, 
	      this->mappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
	      this->quadrature, q_order, 
	      elemGeoInfo, feDataAccessor, ortCoeff, 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, false);
    OrthoMatrixType ortCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, ortCoeff);
    FE::initFE(this->boundaryBasisFunction, p_order, 
              this->boundaryMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	      this->sideQuadrature, q_order, 
              elemGeoInfo, feDataAccessor, ortCoeff);
  }

  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(), false);
    OrthoMatrixType ortCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, ortCoeff);
    FE::initFE(this->boundaryBasisFunction, p_order, 
               this->boundaryMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	       this->sideQuadrature, q_order, 
               elemGeoInfo, feDataAccessor, ortCoeff);
  }

  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);
    OrthoMatrixType elemOrtCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, elemOrtCoeff);
    feDataAccessor.initData(mesh, neighborIt);
    OrthoMatrixType neighborOrtCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, neighborOrtCoeff);
    feDataAccessor.initData(mesh, elemIt, neighborIt, false);
    FE::initFE(this->internalElemBasisFunction, p_order, 
              this->internalNeighborBasisFunction, p_order, 
    	      this->internalMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	      this->sideQuadrature, q_order, 
              elemGeoInfo, feDataAccessor, elemOrtCoeff, neighborOrtCoeff);
  }

  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);
    OrthoMatrixType elemOrtCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, elemOrtCoeff);
    feDataAccessor.initData(mesh, neighbor);
    OrthoMatrixType neighborOrtCoeff;
    FE::initAuxFESpace(p_order, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), q_order, 
                      auxFeManager, elemGeoInfo, feDataAccessor, neighborOrtCoeff);
    
    feDataAccessor.initData(mesh, elem, neighbor, sideIt, false);
    FE::initFE(this->internalElemBasisFunction, p_order, 
              this->internalNeighborBasisFunction, p_order, 
    	      this->internalMappingFunction, ReferenceFiniteElement::getGeoOrder(feDataAccessor.getElemType()), 
    	      this->sideQuadrature, q_order, 
              elemGeoInfo, feDataAccessor, elemOrtCoeff, neighborOrtCoeff);
  }
  
  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());
  }

  MeshFEDataAccessor& getFEDataAccessor()
  {
    return feDataAccessor;
  }

private:
  MeshFEDataAccessor feDataAccessor;
  LiftingOperator<dim> lifting;
  FEManager<BasisFunctionType, MappingFunctionType> auxFeManager;
  ElemGeometryInfo<MappingFunctionType> elemGeoInfo;
};

#endif
