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

#include "BasisFunction.hpp"
#include "MappingBasisFunction.hpp"
#include "CollectiveMappingFunction.hpp"
#include "Quadrature.hpp"
#include "EnumElemType.hpp"
#include "Eigen/Dense"
#include <stdexcept>

template<unsigned int DIm>
class LiftingOperator
{
public:
  typedef Eigen::MatrixXd LiftMatrix;
  typedef Eigen::Matrix<LiftMatrix, DIm, 1> LiftingMatrix;

  LiftingOperator()
    {
      ;
    } 
  
  ~LiftingOperator() { }

  const LiftingMatrix& r11() const 
  {
    return _r11; 
  }

  LiftingMatrix& r11()
  {
    return _r11; 
  }
  
  const LiftingMatrix& r12() const 
  {
    return _r12; 
  }

  LiftingMatrix& r12()
  {
    return _r12; 
  }
  
  const LiftingMatrix& r22() const 
  {
    return _r22; 
  }

  LiftingMatrix& r22()
  {
    return _r22; 
  }
  
  const LiftingMatrix& r21() const 
  {
    return _r21; 
  }

  LiftingMatrix& r21()
  {
    return _r21; 
  }

  template <template <unsigned int> class BasisFunctionType, class MappingFunctionType> 
    void compute(const BasisFunctionType<DIm>& elemBasisFunction, 
                 const BasisFunctionType<DIm>& neighborBasisFunction, 
                 const MappingFunctionType& sideMappingFunction, 
		 const unsigned int firstIndex, const unsigned int lastIndex)
  {
    const typename BasisFunctionType<DIm>::ShapeMatrix& elemShapeMatrix = elemBasisFunction.getShapeFunction();
    const typename BasisFunctionType<DIm>::ShapeMatrix& neighborShapeMatrix = neighborBasisFunction.getShapeFunction();
    for(unsigned int d = 0; d < DIm; d++)
    {
      _r11(d) = LiftMatrix::Zero(elemShapeMatrix.rows(),elemShapeMatrix.rows());
      _r12(d) = LiftMatrix::Zero(elemShapeMatrix.rows(),neighborShapeMatrix.rows());
      _r22(d) = LiftMatrix::Zero(neighborShapeMatrix.rows(),neighborShapeMatrix.rows());
      _r21(d) = LiftMatrix::Zero(neighborShapeMatrix.rows(),elemShapeMatrix.rows());
      for(unsigned int qp = firstIndex; qp < lastIndex; qp++)
      {
	double JWN = 0.5 * sideMappingFunction.getJacobianXweight()[qp] * sideMappingFunction.getNormalVec()(qp)(d);
        _r11(d).noalias() += JWN * (elemShapeMatrix.col(qp) * elemShapeMatrix.col(qp).transpose());
        _r12(d).noalias() -= JWN * (elemShapeMatrix.col(qp) *  neighborShapeMatrix.col(qp).transpose());
        _r22(d).noalias() -= JWN * (neighborShapeMatrix.col(qp) * neighborShapeMatrix.col(qp).transpose());
        _r21(d).noalias() += JWN * (neighborShapeMatrix.col(qp) * elemShapeMatrix.col(qp).transpose());
      }
    }
  }
  
  template <template <unsigned int> class BasisFunctionType, class MappingFunctionType> 
    void compute(const BasisFunctionType<DIm>& elemBasisFunction, 
                 const MappingFunctionType& sideMappingFunction, 
		 const unsigned int firstIndex, const unsigned int lastIndex)
  {
    const typename BasisFunctionType<DIm>::ShapeMatrix& elemShapeMatrix = elemBasisFunction.getShapeFunction();
    for(unsigned int d = 0; d < DIm; d++)
    {
      _r11(d) = LiftMatrix::Zero(elemShapeMatrix.rows(),elemShapeMatrix.rows());
      for(unsigned int qp = firstIndex; qp < lastIndex; qp++)
        _r11(d).noalias() += sideMappingFunction.getJacobianXweight()[qp] * sideMappingFunction.getNormalVec()(qp)(d) * (elemShapeMatrix.col(qp) * elemShapeMatrix.col(qp).transpose());
    }
  }

private:
  
  LiftingMatrix _r11;
  LiftingMatrix _r12;
  LiftingMatrix _r22;
  LiftingMatrix _r21;
  
};

#endif
