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

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

template<class PolyT>
class JacobiFiniteElementBasis
{
public:
  typedef Quadrature::PointVector PointVector;
  
  JacobiFiniteElementBasis() 
    {
      ;
    }
  
  ~JacobiFiniteElementBasis() {}
  
public:
  template<unsigned int Dim> static void initBasisFunction(BasisFunction<Dim>& basisFunction, const PointVector& q_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;
    PointVector points(q_points.rows(), 3);
    CollapsedMapping::mapEN(q_points, points, itemType);
    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 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 
JacobiFiniteElementBasis<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 JacobiFiniteElementBasisService: wrong item type"));
  for (unsigned int j = 0; j < points.rows(); j++)
  {
    double p1 = points(j,0);
    for (unsigned int i = 0; i < basisIndex.size(); i++)
    {
      unsigned int n1 = basisIndex(i)(0);
      shapeMatrix(i,j) = JacobiPoly::evalJacobi(0,0,n1,p1) * sqrt(n1+0.5);
    }
  }
}

template <class PolyT> 
void 
JacobiFiniteElementBasis<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)
  {
    case TRI3:
    case TRI6:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1) * sqrt((n1+0.5)*(n1+n2+1));
        }
      }
      return;
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(0.,0.,n2,p2) * sqrt((n1+0.5)*(n2+0.5));
        }
      }
      return;
    }
    default:
    {
      std::cerr<< "Jacobi2DBasis: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

template <class PolyT> 
void 
JacobiFiniteElementBasis<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)
  {
    case TET4:
    case TET10:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1) * 
                             JacobiPoly::evalJacobi(2.0*(n1+n2+1.0),0.,n3,p3) *pow(0.5*(1.0-p3),n1+n2) * sqrt((n1+0.5)*(n1+n2+1.)*(n1+n2+n3+1.5));
        }
      }
      return;
    }
    case PRISM6:
    case PRISM15:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2)  * pow(0.5*(1.0-p2),n1) * 
			     JacobiPoly::evalJacobi(0.,0.,n3,p3) * sqrt((n1+0.5)*(n1+n2+1)*(n3+0.5));
        }
      }
      return;
    }
    case HEX8:
    case HEX20:
    case HEX27:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(0.,0.,n2,p2) * 
                             JacobiPoly::evalJacobi(0.,0.,n3,p3) * sqrt((n1+0.5)*(n2+0.5)*(n3+0.5));
        }
      }
      return;
    }
    case PYRAMID5:
    case PYRAMID13:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          shapeMatrix(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                             JacobiPoly::evalJacobi(0.,0.,n2,p2) * 
                             JacobiPoly::evalJacobi(2.0*(n1+n2+1.0),0.,n3,p3) *pow(0.5*(1.0-p3),n1+n2) * sqrt((n1+0.5)*(n2+0.5)*(n1+n2+n3+1.5));
        }
      }
      return;
    }
    default:
    {
      std::cerr<< "Jacobi3DBasis: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

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

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

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

template <class PolyT> 
void
JacobiFiniteElementBasis<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 JacobiFiniteElementBasisService: wrong item type"));
  for (unsigned int j = 0; j < points.rows(); j++)
  {
    double p1 = points(j,0);
    for (unsigned int i = 0; i < basisIndex.size(); i++)
    {
      unsigned int n1 = basisIndex(i)(0);
      shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0,0,n1,p1) * sqrt(n1+0.5);
    }
  }
}

template <class PolyT> 
void
JacobiFiniteElementBasis<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)
  {
    case TRI3:
    case TRI6:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          double jacN1 = JacobiPoly::evalJacobi(0.,0.,n1,p1);
          double jacN2 = JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 
                                     JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1-1.0) * sqrt((n1+0.5)*(n1+n2+1));
          shapeDerivMatrix(1)(i,j) = (JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * (1.0+p1) * 0.5 * 
                                      jacN2 * pow(0.5*(1.0-p2),n1-1.0) -
                                      jacN1 * 
                                      jacN2 * 0.5 * n1 * pow(0.5*(1.0-p2),n1-1.0) +
                                      jacN1 * 
                                      JacobiPoly::evalJacobiDeriv(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1)) * sqrt((n1+0.5)*(n1+n2+1));
        }
      }
      return;
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 
                                     JacobiPoly::evalJacobi(0.,0.,n2,p2) * sqrt((n1+0.5)*(n2+0.5));
          shapeDerivMatrix(1)(i,j) = JacobiPoly::evalJacobi(0.,0.,n1,p1) * 
                                     JacobiPoly::evalJacobiDeriv(0.,0.,n2,p2) * sqrt((n1+0.5)*(n2+0.5));
        }
      }
      return;
    }
    default:
    {
      std::cerr<< "Jacobi Basis Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

template <class PolyT> 
void
JacobiFiniteElementBasis<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)
  {
    case TET4:
    case TET10:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          double jacN1 = JacobiPoly::evalJacobi(0.,0.,n1,p1);
          double jacN2 = JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2);
          double jacN3 = JacobiPoly::evalJacobi(2.0*(n1+n2+1.),0.,n3,p3);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) *
                                     jacN2 * pow(0.5*(1.0-p2),n1-1.0) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) * sqrt((n1+0.5)*(n1+n2+1.)*(n1+n2+n3+1.5));
          shapeDerivMatrix(1)(i,j) = (JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * (1.0+p1) * 0.5 *
                                     jacN2 * pow(0.5*(1.0-p2),n1-1.0) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) -
                                     jacN1 * 
                                     jacN2 * 0.5 * n1 * pow(0.5*(1.0-p2),n1-1.0) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) + 
                                     jacN1 * 
                                     JacobiPoly::evalJacobiDeriv(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0)) * sqrt((n1+0.5)*(n1+n2+1)*(n1+n2+n3+1.5));
          shapeDerivMatrix(2)(i,j) = (JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * (1.0+p1) * 0.5 *
                                     jacN2 * pow(0.5*(1.0-p2),n1-1.0) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) -
                                     jacN1 * (1.0+p2) * 0.25 *
                                     jacN2 * n1 * pow(0.5*(1.0-p2),n1-1.0) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) +
                                     jacN1 * (1.0+p2) * 0.5 *
                                     JacobiPoly::evalJacobiDeriv(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1) * 
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1.0) -
                                     jacN1 * 
                                     jacN2 * pow(0.5*(1.0-p2),n1) * 
                                     jacN3 * 0.5 * (n1+n2) *pow(0.5*(1.0-p3),n1+n2-1.0) +
                                     jacN1 * 
                                     jacN2 * pow(0.5*(1.0-p2),n1) * 
                                     JacobiPoly::evalJacobiDeriv(2.0*(n1+n2+1.0),0.,n3,p3) * pow(0.5*(1.0-p3),n1+n2)) * sqrt((n1+0.5)*(n1+n2+1)*(n1+n2+n3+1.5));
        }
      }
      return;
    }
    case HEX8:
    case HEX20:
    case HEX27:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          double jacN1 = JacobiPoly::evalJacobi(0.,0.,n1,p1);
          double jacN2 = JacobiPoly::evalJacobi(0.,0.,n2,p2);
          double jacN3 = JacobiPoly::evalJacobi(0.,0.,n3,p3);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 
                                     jacN2 *
                                     jacN3 * sqrt((n1+0.5)*(n2+0.5)*(n3+0.5));
          shapeDerivMatrix(1)(i,j) = jacN1 * 
                                     JacobiPoly::evalJacobiDeriv(0.,0.,n2,p2) *
                                     jacN3 * sqrt((n1+0.5)*(n2+0.5)*(n3+0.5));
          shapeDerivMatrix(2)(i,j) = jacN1 * 
                                     jacN2 *
                                     JacobiPoly::evalJacobiDeriv(0.,0.,n3,p3) * sqrt((n1+0.5)*(n2+0.5)*(n3+0.5));
        }
      } 
      return;
    }
    case PRISM6:
    case PRISM15:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          double jacN1 = JacobiPoly::evalJacobi(0.,0.,n1,p1);
          double jacN2 = JacobiPoly::evalJacobi(2.0*n1+1.,0.,n2,p2);
          double jacN3 = JacobiPoly::evalJacobi(0.,0.,n3,p3);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 
                                     jacN2 * pow(0.5*(1.0-p2),n1-1.0) * 
			             jacN3 * sqrt((n1+0.5)*(n1+n2+1)*(n3+0.5));
          shapeDerivMatrix(1)(i,j) = (JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * (1.0+p1) * 0.5 * 
                                     jacN2 * pow(0.5*(1.0-p2),n1-1.0) *
			             jacN3 -
                                     jacN1 * 
                                     jacN2 * 0.5 * n1 * pow(0.5*(1.0-p2),n1-1.0) *
			             jacN3 +
                                     jacN1 * 
                                     JacobiPoly::evalJacobiDeriv(2.0*n1+1.,0.,n2,p2) * pow(0.5*(1.0-p2),n1) *
			             jacN3) * sqrt((n1+0.5)*(n1+n2+1)*(n3+0.5));
          shapeDerivMatrix(2)(i,j) = jacN1 * 
                                     jacN2 * pow(0.5*(1.0-p2),n1) * 
			             JacobiPoly::evalJacobiDeriv(0.,0.,n3,p3) * sqrt((n1+0.5)*(n1+n2+1)*(n3+0.5));
        }
      }
      return;
    }
    case PYRAMID5:
    case PYRAMID13:
    {
      for (unsigned int j = 0; j < points.rows(); j++)
      {
        double p1 = points(j,0);
        double p2 = points(j,1);
        double p3 = points(j,2);
        for (unsigned int i = 0; i < basisIndex.size(); i++)
        {
          unsigned int n1 = basisIndex(i)(0);
          unsigned int n2 = basisIndex(i)(1);
          unsigned int n3 = basisIndex(i)(2);
          double jacN1 = JacobiPoly::evalJacobi(0.,0.,n1,p1);
          double jacN2 = JacobiPoly::evalJacobi(0.,0.,n2,p2);
          double jacN3 = JacobiPoly::evalJacobi(2.0*(n1+n2+1.),0.,n3,p3);
          shapeDerivMatrix(0)(i,j) = JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 
                                     jacN2 *
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1) * sqrt((n1+0.5)*(n2+0.5)*(n1+n2+n3+1.5));
          shapeDerivMatrix(1)(i,j) = jacN1 * 
                                     JacobiPoly::evalJacobiDeriv(0.,0.,n2,p2) *
                                     jacN3 * pow(0.5*(1.0-p3),n1+n2-1) * sqrt((n1+0.5)*(n2+0.5)*(n1+n2+n3+1.5));
          shapeDerivMatrix(2)(i,j) = (JacobiPoly::evalJacobiDeriv(0.,0.,n1,p1) * 0.5 * (1. + p1) *
                                     jacN2 *
				     jacN3 * pow(0.5*(1.0-p3),n1+n2-1) +
				     jacN1 *
				     JacobiPoly::evalJacobiDeriv(0.,0.,n2,p2) * 0.5 * (1. + p2) *
				     jacN3 * pow(0.5*(1.0-p3),n1+n2-1) -
				     jacN1 *
				     jacN2 *
				     jacN3 * 0.5 * (n1 + n2) * pow(0.5*(1.0-p3),n1+n2-1) +
				     jacN1 * 
				     jacN2 * 
				     JacobiPoly::evalJacobiDeriv(2.0*(n1+n2+1.),0.,n3,p3) * pow(0.5*(1.0-p3),n1+n2)) * sqrt((n1+0.5)*(n2+0.5)*(n1+n2+n3+1.5));
        }
      } 
      return;
    }
    default:
    {
      std::cerr<< "Jacobi Basis Service: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

#endif 
