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

#include "BasisFunction.hpp"
#include "EnumElemType.hpp"
#include "Tuple.hpp"
#include <stdexcept>
#include "Quadrature.hpp"
#include "CollapsedMapping.hpp"

template<class PolyT>
class BernsteinFiniteElementBasis
{
public:
  typedef Quadrature::PointVector PointVector;
  
  BernsteinFiniteElementBasis() 
    {
      ;
    }
  
  ~BernsteinFiniteElementBasis() {}
  
public:
  template<unsigned int Dim> static void initBasisFunction(BasisFunction<Dim>& basisFunction, const PointVector& points)
  {
    const ElemType itemType = basisFunction.getItemType();
    const unsigned int b_order = basisFunction.getOrder();
    Eigen::Matrix<Eigen::Matrix<unsigned int,Dim,1>,Eigen::Dynamic,1> basisIndex;
    PolyT::computeBasisIndex(basisIndex, itemType, b_order);
    if (basisFunction.computeShape())
    {
      typename BasisFunction<Dim>::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
      basisFunction.resizeShapeFunction(basisIndex.size(), points.rows());
      evaluateFunction(shapeMatrix, basisIndex, itemType, points);
    }
    if (basisFunction.computeShapeDeriv())
    {
      typename BasisFunction<Dim>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
      basisFunction.resizeShapeFunctionDeriv(basisIndex.size(), points.rows());
      evaluateDerivative(shapeDerivMatrix, basisIndex, itemType, points);
    }
  }

  static unsigned int factorial(unsigned int n)
  {
    unsigned int factorial_n = 1;
    
    if (n==0)
     return factorial_n;
    
    for (unsigned int i=1; i<n; i++)
     factorial_n *= i+1;
    
     return factorial_n;
  }

  static unsigned int dofsNumber(const unsigned int Dim, const unsigned int order, const ElemType itemType);
  
  static bool is_hierarchical();

private:
  //! Evaluate \$\varphi_i(x)\$
  static void evaluateFunction(typename BasisFunction<1>::ShapeMatrix& shapeMatrix, 
                               const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>& basisIndex, 
	                       const ElemType itemType,
		    	       const PointVector& points);
  static void evaluateFunction(typename BasisFunction<2>::ShapeMatrix& shapeMatrix, 
                               const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, 
	                       const ElemType itemType,
			       const PointVector& points);
  static void evaluateFunction(typename BasisFunction<3>::ShapeMatrix& shapeMatrix, 
                               const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>& basisIndex, 
	                       const ElemType itemType,
			       const PointVector& points);

  //! Evaluate \$\partial_j \varphi_i(x)\$
  static void evaluateDerivative(typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix, 
                                 const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>& basisIndex, 
			         const ElemType itemType,
			         const PointVector& points);
  static void evaluateDerivative(typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix, 
                                 const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, 
			         const ElemType itemType,
			         const PointVector& points);
  static void evaluateDerivative(typename BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix, 
                                 const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>& basisIndex, 
			         const ElemType itemType,
			         const PointVector& points);
};

template <class PolyT> 
void 
BernsteinFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<1>::ShapeMatrix& shapeMatrix, 
                 const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>&  basisIndex, 
	         const ElemType itemType,
	         const PointVector& points)
{
  assert(itemType == EDGE2 || itemType == EDGE3);//, ("Error in BernsteinFiniteElementBasisService: wrong item type"));
  const int p_order = basisIndex.size()-1;
  for (unsigned int j = 0; j < points.rows(); j++)
  {
    unsigned int binomial_p_i = 1;
    double p1 = points(j,0);
    shapeMatrix(0,j) = binomial_p_i * std::pow((1-p1)/2,static_cast<int>(p_order));
    shapeMatrix(1,j) = binomial_p_i * std::pow((1+p1)/2,static_cast<int>(p_order));
    for (unsigned int i = 2; i < basisIndex.size(); i++)
    {
      const int m = p_order-basisIndex(i)(0)+1;
      const int n = (basisIndex(i)(0)-1);
      binomial_p_i = factorial(p_order)/(factorial(n)*factorial(p_order-n));
      shapeMatrix(i,j) = binomial_p_i * std::pow((1+p1)/2,n) * std::pow((1-p1)/2,m);
    }
  }
}

template <class PolyT> 
void 
BernsteinFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<2>::ShapeMatrix& shapeMatrix, 
                 const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>&  basisIndex, 
                 const ElemType itemType, 
	         const PointVector& points)
{
  switch (itemType)
  {
    default:
    {
      std::cerr<< "Bernstein2DBasis: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

template <class PolyT> 
void 
BernsteinFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<3>::ShapeMatrix& shapeMatrix, 
                 const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>&  basisIndex, 
	         const ElemType itemType, 
	         const PointVector& points)
{
  switch (itemType)
  {
    default:
    {
      std::cerr<< "Bernstein3DBasis: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

template <class PolyT> 
unsigned int 
BernsteinFiniteElementBasis<PolyT>:: 
dofsNumber(const unsigned int Dim, const unsigned int order, const ElemType itemType) 
{
  return PolyT::modesNumber(Dim, order, itemType);
}

template <class PolyT> bool
BernsteinFiniteElementBasis<PolyT>:: 
is_hierarchical() 
{
  return true;
}

/*---------------------------------------------------------------------------*/

template <class PolyT> 
void
BernsteinFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix, 
                   const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>& basisIndex, 
		   const ElemType itemType,
	           const PointVector& points)
{
  assert(itemType == EDGE2 || itemType == EDGE3);//, ("Error in BernsteinFiniteElementBasisService: wrong item type"));
  const int p_order = basisIndex.size()-1;
  for (unsigned int j = 0; j < points.rows(); j++)
  {
    unsigned int binomial_p_i = 1;
    double p1 = points(j,0);
    shapeDerivMatrix(0)(0,j) = binomial_p_i * (-0.5) * p_order * std::pow((1-p1)/2,static_cast<int>(p_order-1));
    shapeDerivMatrix(0)(1,j) = binomial_p_i * 0.5    * p_order * std::pow((1+p1)/2,static_cast<int>(p_order-1));
    for (unsigned int i = 2; i < basisIndex.size(); i++)
    {
      const int m = p_order-basisIndex(i)(0)+1;
      const int n = (basisIndex(i)(0)-1);
      binomial_p_i = factorial(p_order)/(factorial(n)*factorial(p_order-n));
      shapeDerivMatrix(0)(i,j) = binomial_p_i * (1./2. * n * std::pow((1+p1)/2,n-1) * std::pow((1-p1)/2,m)
				               - 1./2. * m * std::pow((1+p1)/2,n)   * std::pow((1-p1)/2,m-1));
    }
  }
}

template <class PolyT> 
void
BernsteinFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix, 
                   const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, 
		   const ElemType itemType,
	           const PointVector& points)
{
  switch (itemType)
  {
    default:
    {
      std::cerr<< "Bernstein Basis Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

template <class PolyT> 
void
BernsteinFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix, 
                   const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>& basisIndex, 
		   const ElemType itemType,
	           const PointVector& points)
{
  switch (itemType)
  {
    default:
    {
      std::cerr<< "Bernstein Basis Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

#endif 
