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

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

template<class PolyT>
class MonomialFiniteElementBasis
{
public:
  typedef Quadrature::PointVector PointVector;
  
  MonomialFiniteElementBasis() 
    {
      ;
    }
  
  ~MonomialFiniteElementBasis() {}
  
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, points);
    }
    if (basisFunction.computeShapeDeriv())
    {
      typename BasisFunction<Dim>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
      basisFunction.resizeShapeFunctionDeriv(basisIndex.size(), points.rows());
      evaluateDerivative(shapeDerivMatrix, basisIndex, points);
    }
  }

  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 PointVector& points);
  static void evaluateFunction(typename BasisFunction<2>::ShapeMatrix& shapeMatrix, const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, 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 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 PointVector & points);
  static void evaluateDerivative(typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix, 
                                 const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, 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 PointVector & points);
};

template <class PolyT> 
void 
MonomialFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<1>::ShapeMatrix& shapeMatrix, const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>&  basisIndex, const PointVector& points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeMatrix(z,shapeQPindex) = 1;
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeMatrix(z,shapeQPindex) *= p1;
    }
  }
}

template <class PolyT> 
void 
MonomialFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<2>::ShapeMatrix& shapeMatrix, const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>&  basisIndex, const PointVector& points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    double p2 = points(shapeQPindex,1);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeMatrix(z,shapeQPindex) = 1;
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeMatrix(z,shapeQPindex) *= p1;
      for (unsigned int j = 0; j < basisIndex(z)(1); j++)
	shapeMatrix(z,shapeQPindex) *= p2;
    }
  }
}

template <class PolyT> 
void 
MonomialFiniteElementBasis<PolyT>:: 
evaluateFunction(typename BasisFunction<3>::ShapeMatrix& shapeMatrix, const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>&  basisIndex, const PointVector& points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    double p2 = points(shapeQPindex,1);
    double p3 = points(shapeQPindex,2);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeMatrix(z,shapeQPindex) = 1;
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeMatrix(z,shapeQPindex) *= p1;
      for (unsigned int j = 0; j < basisIndex(z)(1); j++)
	shapeMatrix(z,shapeQPindex) *= p2;
      for (unsigned int k = 0; k < basisIndex(z)(2); k++)
	shapeMatrix(z,shapeQPindex) *= p3;
    }
  }
}

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

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

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

template <class PolyT> 
void
MonomialFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix, 
                          const Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>& basisIndex, const PointVector & points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeDerivMatrix(0)(z,shapeQPindex) = basisIndex(z)(0);
      for (unsigned int i = 1; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p1;
    }
  }
}

template <class PolyT> 
void
MonomialFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix, 
                          const Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>& basisIndex, const PointVector & points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    double p2 = points(shapeQPindex,1);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeDerivMatrix(0)(z,shapeQPindex) = basisIndex(z)(0);
      for (unsigned int i = 1; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p1;
      for (unsigned int j = 0; j < basisIndex(z)(1); j++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p2;

      shapeDerivMatrix(1)(z,shapeQPindex) = basisIndex(z)(1);
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(1)(z,shapeQPindex) *= p1;
      for (unsigned int j = 1; j < basisIndex(z)(1); j++)
	shapeDerivMatrix(1)(z,shapeQPindex) *= p2;
    }
  }
}

template <class PolyT> 
void
MonomialFiniteElementBasis<PolyT>:: 
evaluateDerivative(typename BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix, 
                          const Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>& basisIndex, const PointVector & points)
{
  for (unsigned int shapeQPindex = 0; shapeQPindex < points.rows(); shapeQPindex++)
  {
    double p1 = points(shapeQPindex,0);
    double p2 = points(shapeQPindex,1);
    double p3 = points(shapeQPindex,2);
    for (unsigned int z = 0; z < basisIndex.size(); z++)
    {
      shapeDerivMatrix(0)(z,shapeQPindex) = basisIndex(z)(0);
      for (unsigned int i = 1; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p1;
      for (unsigned int j = 0; j < basisIndex(z)(1); j++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p2;
      for (unsigned int k = 0; k < basisIndex(z)(2); k++)
	shapeDerivMatrix(0)(z,shapeQPindex) *= p3;

      shapeDerivMatrix(1)(z,shapeQPindex) = basisIndex(z)(1);
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(1)(z,shapeQPindex) *= p1;
      for (unsigned int j = 1; j < basisIndex(z)(1); j++)
	shapeDerivMatrix(1)(z,shapeQPindex) *= p2;
      for (unsigned int k = 0; k < basisIndex(z)(2); k++)
	shapeDerivMatrix(1)(z,shapeQPindex) *= p3;
      
      shapeDerivMatrix(2)(z,shapeQPindex) = basisIndex(z)(2);
      for (unsigned int i = 0; i < basisIndex(z)(0); i++)
	shapeDerivMatrix(2)(z,shapeQPindex) *= p1;
      for (unsigned int j = 0; j < basisIndex(z)(1); j++)
	shapeDerivMatrix(2)(z,shapeQPindex) *= p2;
      for (unsigned int k = 1; k < basisIndex(z)(2); k++)
	shapeDerivMatrix(2)(z,shapeQPindex) *= p3;
    }
  } 
}
#endif 
