// 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/. */

#include "ExplicitLagrangianFiniteElementBasis.hpp"
#include "BasisFunction.hpp"
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

bool 
ExplicitLagrangianFiniteElementBasis::
is_hierarchical() 
{
  return false;
}

unsigned int 
ExplicitLagrangianFiniteElementBasis::
nbShapeFunction(const unsigned int order, const ElemType itemType) 
{
  switch (itemType)
  {
    case NODEELEM:
    {
      return 1;
    }
    case EDGE2:
    case EDGE3:
    {
      return order+1;
    }
    case TRI3:
    case TRI6:
    {
      return (order+1)*(order+2)/2;
    }
    case QUAD4:
    case QUAD9:
    {
      return (order+1)*(order+1);
    }
    case QUAD8:
    {
      switch (order)
      {
        case 1:  
          return 4;
        case 2:
          return 8;
        default:
        {
          std::cerr<< "ExplicitLagrangianBasis: Unsupported order "<<itemType<<std::endl;
          std::logic_error("halt");
          return 0;
        }
      }
    }
    case TET4:
    case TET10:
    {
      return (order+1)*(order+2)*(order+3)/6;
    }
    case HEX8:
    case HEX27:
    {
      return (order+1)*(order+1)*(order+1);
    }
    case HEX20:
    {
      switch (order)
      {
        case 1:  
          return 8;
        case 2:
          return 20;
        default:
        {
          std::cerr<< "ExplicitLagrangianBasis: Unsupported order "<<itemType<<std::endl;
          std::logic_error("halt");
          return 0;
        }
      }
    }
    case PRISM6:
    {
      return 6;
    }
    case PRISM15:
    {
      switch (order)
      {
        case 1:  
          return 6;
        case 2:
          return 15;
        default:
        {
          std::cerr<< "ExplicitLagrangianBasis: Unsupported order "<<itemType<<std::endl;
          std::logic_error("halt");
          return 0;
        }
      }
    }
    case PYRAMID5:
    {
      return 5;
    }
    case PYRAMID13:
    {
      switch (order)
      {
        case 1:  
          return 5;
        case 2:
          return 13;
        default:
        {
          std::cerr<< "ExplicitLagrangianBasis: Unsupported order "<<itemType<<std::endl;
          std::logic_error("halt");
          return 0;
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangianBasis: Unsupported Itemtype "<<itemType<<std::endl;
      std::logic_error("halt");
      return 0;
    }
  }
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateFunction(BasisFunction<0>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  assert(itemType == NODEELEM);//, ("Error in MonomialFiniteElementBasisService: wrong item type"));
  assert(points.rows() == 1);
  BasisFunction<0>::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
  basisFunction.resizeShapeFunction(1,1);
  for (unsigned int j = 0; j < 1; j++)
    shapeMatrix(0,j) = 1.;
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateFunction(BasisFunction<1>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<1>::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
  basisFunction.resizeShapeFunction(order+1, points.rows());
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {
        case EDGE2:
        case EDGE3:
        {
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            shapeMatrix(0,j) = .5*(1. - p1);
            shapeMatrix(1,j) = .5*(1. + p1);
	  }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {
        case EDGE3:
        {
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            shapeMatrix(0,j) = .5*p1*(p1 - 1.);
            shapeMatrix(1,j) = .5*p1*(p1 + 1.);
            shapeMatrix(2,j) = 1. - p1*p1;
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis: Unsupported order "<<order;
      std::logic_error("halt");
    }
  }
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateFunction(BasisFunction<2>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<2>::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {
        case TRI3:
        case TRI6:
        {
          basisFunction.resizeShapeFunction(3, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeMatrix(0,j) = .5*(- p1 - p2);
            shapeMatrix(1,j) = .5*(p1 + 1.);
            shapeMatrix(2,j) = .5*(p2 + 1.);
	  }
	  return;
        }
        case QUAD4:
        case QUAD8:
        case QUAD9:
        {
          basisFunction.resizeShapeFunction(4, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeMatrix(0,j) = .25*(1. - p1)*(1. - p2);
            shapeMatrix(1,j) = .25*(1. + p1)*(1. - p2);
            shapeMatrix(2,j) = .25*(1. + p1)*(1. + p2);
            shapeMatrix(3,j) = .25*(1. - p1)*(1. + p2);
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian 2D Basis: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {
        case TRI6:
        {
          basisFunction.resizeShapeFunction(6, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeMatrix(0,j) = .5*(p1 + p2)*(p1 + p2 + 1.);
            shapeMatrix(1,j) = .5 * p1 * (p1 + 1.);
            shapeMatrix(2,j) = .5 * p2 * (p2 + 1.);
            shapeMatrix(3,j) = (-p1 -p2) * (p1 + 1.);
            shapeMatrix(4,j) = (p1 + 1.) * (p2 + 1.);
            shapeMatrix(5,j) = (-p1 -p2) * (p2 + 1.);
          }
	  return;
	}
        case QUAD9:
        {
          basisFunction.resizeShapeFunction(9, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeMatrix(0,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 - 1.);
            shapeMatrix(1,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 - 1.);
            shapeMatrix(2,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 + 1.);
            shapeMatrix(3,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 + 1.);
            shapeMatrix(4,j) = .5 * (1. - p1*p1) * p2 * (p2 - 1.);
            shapeMatrix(5,j) = .5 * p1 * (p1 + 1) * (1. - p2*p2);
            shapeMatrix(6,j) = .5 * (1. - p1*p1) * p2 * (p2 + 1.);
            shapeMatrix(7,j) = .5 * p1 * (p1 - 1.) * (1. - p2*p2);
            shapeMatrix(8,j) = (1. - p1*p1) * (1. - p2*p2);
          }
	  return;
	}
	case QUAD8:
	{
          basisFunction.resizeShapeFunction(8, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeMatrix(0,j) = .25*(1. - p1)*(1. - p2)*(-1. - p1 - p2);
            shapeMatrix(1,j) = .25*(1. + p1)*(1. - p2)*(-1. + p1 - p2);
            shapeMatrix(2,j) = .25*(1. + p1)*(1. + p2)*(-1. + p1 + p2);
            shapeMatrix(3,j) = .25*(1. - p1)*(1. + p2)*(-1. - p1 + p2);
            shapeMatrix(4,j) = .5*(1. - p1*p1)*(1. - p2);
            shapeMatrix(5,j) = .5*(1. + p1)*(1. - p2*p2);
            shapeMatrix(6,j) = .5*(1. - p1*p1)*(1. + p2);
            shapeMatrix(7,j) = .5*(1. - p1)*(1. - p2*p2);
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis: Unsupported order "<<order;
      std::logic_error("halt");
    }
  }
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateFunction(BasisFunction<3>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<3>::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {
        case TET4:
        case TET10:
        {
          basisFunction.resizeShapeFunction(4, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeMatrix(0,j) = .5*(- p1 - p2 - p3 - 1);
            shapeMatrix(1,j) = .5*(1. + p1);
            shapeMatrix(2,j) = .5*(1. + p2);
            shapeMatrix(3,j) = .5*(1. + p3);
          }
	  return;
        }
        case HEX8:
        case HEX20:
        case HEX27:
        {
          basisFunction.resizeShapeFunction(8, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeMatrix(0,j) = .125*(1. - p1)*(1. - p2)*(1. - p3);
            shapeMatrix(1,j) = .125*(1. + p1)*(1. - p2)*(1. - p3);
            shapeMatrix(2,j) = .125*(1. + p1)*(1. + p2)*(1. - p3);
            shapeMatrix(3,j) = .125*(1. - p1)*(1. + p2)*(1. - p3);
            shapeMatrix(4,j) = .125*(1. - p1)*(1. - p2)*(1. + p3);
            shapeMatrix(5,j) = .125*(1. + p1)*(1. - p2)*(1. + p3);
            shapeMatrix(6,j) = .125*(1. + p1)*(1. + p2)*(1. + p3);
            shapeMatrix(7,j) = .125*(1. - p1)*(1. + p2)*(1. + p3);
	  }
	  return;
	}
        case PRISM6:
        case PRISM15:
        {
          basisFunction.resizeShapeFunction(6, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeMatrix(0,j) = .25*(1. - p3)*(-p1 - p2);
            shapeMatrix(1,j) = .25*(1. - p3)*(p1 + 1.);
            shapeMatrix(2,j) = .25*(1. - p3)*(p2 + 1.);
            shapeMatrix(3,j) = .25*(1. + p3)*(-p1 - p2);
            shapeMatrix(4,j) = .25*(1. + p3)*(p1 + 1.);
            shapeMatrix(5,j) = .25*(1. + p3)*(p2 + 1.);
          }
	  return;
        }
	case PYRAMID5:
	case PYRAMID13:
	{
          basisFunction.resizeShapeFunction(5, points.rows());
	  double eps = 1.e-30;
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
	    shapeMatrix(0,j) = (0.5 * (-p3 - p1) * (-p3 - p2)) / (1. - p3 + eps);
	    shapeMatrix(1,j) = (0.5 * (1. + p1) * (-p3 - p2))  / (1. - p3 + eps);
	    shapeMatrix(2,j) = (0.5 * (1. + p1) * (1. + p2))   / (1. - p3 + eps);
	    shapeMatrix(3,j) = (0.5 * (-p3 - p1) * (1. + p2))  / (1. - p3 + eps);
	    shapeMatrix(4,j) = 0.5 * (1. + p3);
	  }
	  return;
	}
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {
	case PYRAMID13:
	{
          basisFunction.resizeShapeFunction(13, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
	    double oneMp3_sq = (1.-p3)*(1.-p3)+1.e-30;
	    shapeMatrix(0,j) = (p1 + p3) * (-p2 -p3) * (p1 + p2 + ((1. - p3) * 0.5 * p3) + 2.)        / oneMp3_sq;
	    shapeMatrix(1,j) = (1. + p1) * (-p2 -p3) * (p1 - p2 - ((1. - p3) * (0.5 * p3 + 1.)))      / oneMp3_sq;
	    shapeMatrix(2,j) = (1. + p1) * (1. + p2) * (p1 + p2 - ((1. - p3) * (0.5 * p3 + 2.)) + 2.) / oneMp3_sq;
	    shapeMatrix(3,j) = (p1 + p3) * (1. + p2) * (p1 - p2 + ((1. - p3) * (0.5 * p3 + 1.)))      / oneMp3_sq;
	    shapeMatrix(4,j) = 0.5 * p3 * (1. + p3);
	    shapeMatrix(5,j) = -2. * (p1 + 1) * (p3 + p1) * (-p2 -p3) / oneMp3_sq; 
	    shapeMatrix(6,j) = -2. * (p2 + 1) * (p3 + p2) * (1. + p1) / oneMp3_sq;
	    shapeMatrix(7,j) = -2. * (p1 + 1) * (p3 + p1) * (1. + p2) / oneMp3_sq;
	    shapeMatrix(8,j) = -2. * (p2 + 1) * (p3 + p2) * (-p1 -p3) / oneMp3_sq;
	    shapeMatrix(9,j)  = (-p1 -p3) * (-p2 -p3) * (1. - p3*p3)  / oneMp3_sq;
	    shapeMatrix(10,j) = (1. + p1) * (-p2 -p3) * (1. - p3*p3)  / oneMp3_sq;
	    shapeMatrix(11,j) = (1. + p1) * (1. + p2) * (1. - p3*p3)  / oneMp3_sq;
	    shapeMatrix(12,j) = (-p1 -p3) * (1. + p2) * (1. - p3*p3)  / oneMp3_sq;
	  }
	  return;
	}
	case PRISM15:
	{
          basisFunction.resizeShapeFunction(15, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeMatrix(0,j) = 0.25 * (-p1 - p2) * (1. - p3) * (-2. - p1 - p2 - p3);
            shapeMatrix(1,j) = 0.25 * (1 + p1) * (1. - p3) * (-1. + p1 - p3);
            shapeMatrix(2,j) = 0.25 * (1 + p2) * (1. - p3) * (-1. + p2 - p3);
            shapeMatrix(3,j) = 0.25 * (-p1 - p2) * (1. + p3) * (-2. - p1 - p2 + p3);
            shapeMatrix(4,j) = 0.25 * (1. + p1) * (1. + p3) * (-1. + p1 + p3);
            shapeMatrix(5,j) = 0.25 * (1. + p2) * (1. + p3) * (-1. + p2 + p3);
	    shapeMatrix(6,j) = 0.5 * (1. + p1) * (-p1 - p2) * (1 - p3);
	    shapeMatrix(7,j) = 0.5 * (1. + p1) * (1. + p2) * (1 - p3);
	    shapeMatrix(8,j) = 0.5 * (1. + p2) * (-p1 - p2) * (1 - p3);
	    shapeMatrix(9,j) = 0.5 * (1. + p3) * (-p1 - p2) * (1 - p3);
	    shapeMatrix(10,j) = 0.5 * (1. + p3) * (1. + p1) * (1 - p3);
	    shapeMatrix(11,j) = 0.5 * (1. + p3) * (1. + p2) * (1 - p3);
	    shapeMatrix(12,j) = 0.5 * (1. + p1) * (-p1 - p2) * (1 + p3);
	    shapeMatrix(13,j) = 0.5 * (1. + p1) * (1. + p2) * (1 + p3);
	    shapeMatrix(14,j) = 0.5 * (1. + p2) * (-p1 - p2) * (1 + p3);
	  }
	  return;
	}

	case TET10:
	{
          basisFunction.resizeShapeFunction(10, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
	    shapeMatrix(0,j) = 0.5 * (1. + p1 + p2 + p3) * (2. + p1 + p2 + p3);
            shapeMatrix(1,j) = 0.5 * (p1 + 1.) * p1;
            shapeMatrix(2,j) = 0.5 * (p2 + 1.) * p2;
            shapeMatrix(3,j) = 0.5 * (p3 + 1.) * p3;
	    shapeMatrix(4,j) = (-1. - p1 - p2 - p3) * (p1 + 1.);
            shapeMatrix(5,j) = (p1 + 1.) * (p2 + 1.);
	    shapeMatrix(6,j) = (-1. - p1 - p2 - p3) * (p2 + 1.);
	    shapeMatrix(7,j) = (-1. - p1 - p2 - p3) * (p3 + 1.);
            shapeMatrix(8,j) = (p1 + 1.) * (p3 + 1.);
            shapeMatrix(9,j) = (p2 + 1.) * (p3 + 1.);
          }
	  return;
	}
        case HEX20:
        {
          basisFunction.resizeShapeFunction(20, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
	    shapeMatrix(0,j) = 0.125 * (1. - p1) * (1. - p2) * (1. - p3) * (-2. - p1 - p2 - p3);
	    shapeMatrix(1,j) = 0.125 * (1. + p1) * (1. - p2) * (1. - p3) * (-2. + p1 - p2 - p3);
	    shapeMatrix(2,j) = 0.125 * (1. + p1) * (1. + p2) * (1. - p3) * (-2. + p1 + p2 - p3);
	    shapeMatrix(3,j) = 0.125 * (1. - p1) * (1. + p2) * (1. - p3) * (-2. - p1 + p2 - p3);
	    shapeMatrix(4,j) = 0.125 * (1. - p1) * (1. - p2) * (1. + p3) * (-2. - p1 - p2 + p3);
	    shapeMatrix(5,j) = 0.125 * (1. + p1) * (1. - p2) * (1. + p3) * (-2. + p1 - p2 + p3);
	    shapeMatrix(6,j) = 0.125 * (1. + p1) * (1. + p2) * (1. + p3) * (-2. + p1 + p2 + p3);
	    shapeMatrix(7,j) = 0.125 * (1. - p1) * (1. + p2) * (1. + p3) * (-2. - p1 + p2 + p3);
	    shapeMatrix(8,j) =  0.25 * (1. - p1*p1) * (1 - p2) * (1 - p3);
	    shapeMatrix(9,j) =  0.25 * (1. - p2*p2) * (1 + p1) * (1 - p3);
	    shapeMatrix(10,j) = 0.25 * (1. - p1*p1) * (1 + p2) * (1 - p3);
	    shapeMatrix(11,j) = 0.25 * (1. - p2*p2) * (1 - p1) * (1 - p3);
	    shapeMatrix(12,j) = 0.25 * (1. - p3*p3) * (1 - p1) * (1 - p2);
	    shapeMatrix(13,j) = 0.25 * (1. - p3*p3) * (1 + p1) * (1 - p2);
	    shapeMatrix(14,j) = 0.25 * (1. - p3*p3) * (1 + p1) * (1 + p2);
	    shapeMatrix(15,j) = 0.25 * (1. - p3*p3) * (1 - p1) * (1 + p2);
	    shapeMatrix(16,j) = 0.25 * (1. - p1*p1) * (1 - p2) * (1 + p3);
	    shapeMatrix(17,j) = 0.25 * (1. - p2*p2) * (1 + p1) * (1 + p3);
	    shapeMatrix(18,j) = 0.25 * (1. - p1*p1) * (1 + p2) * (1 + p3);
	    shapeMatrix(19,j) = 0.25 * (1. - p2*p2) * (1 - p1) * (1 + p3);
          }
	  return;
        }
        case HEX27:
        {
          basisFunction.resizeShapeFunction(27, points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeMatrix(0,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeMatrix(1,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeMatrix(2,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeMatrix(3,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeMatrix(4,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeMatrix(5,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeMatrix(6,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeMatrix(7,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeMatrix(8,j) = .25 * (1. - p1*p1) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeMatrix(9,j) = .25 * p1 * (p1 + 1.) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeMatrix(10,j) = .25 * (1. - p1*p1) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeMatrix(11,j) = .25 * p1 * (p1 - 1.) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeMatrix(12,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeMatrix(13,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeMatrix(14,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeMatrix(15,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeMatrix(16,j) = .25 * (1. - p1*p1) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeMatrix(17,j) = .25 * p1 * (p1 + 1.) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeMatrix(18,j) = .25 * (1. - p1*p1) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeMatrix(19,j) = .25 * p1 * (p1 - 1.) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeMatrix(20,j) = .5 * (1. - p1*p1) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeMatrix(21,j) = .5 * p1 * (p1 + 1.) * (1. - p2*p2) * (1. - p3*p3);
            shapeMatrix(22,j) = .5 * (1. - p1*p1) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeMatrix(23,j) = .5 * p1 * (p1 - 1.) * (1. - p2*p2) * (1. - p3*p3);
            shapeMatrix(24,j) = .5 * (1. - p1*p1) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeMatrix(25,j) = .5 * (1. - p1*p1) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeMatrix(26,j) = (1. - p1*p1) * (1. - p2*p2) * (1. - p3*p3);
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis Service: Unsupported order "<<order;
      std::logic_error("halt");
    }
  }
}

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


void
ExplicitLagrangianFiniteElementBasis::
evaluateDerivative(BasisFunction<0>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  assert(itemType == NODEELEM);//, ("Error in MonomialFiniteElementBasisService: wrong item type"));
  assert(points.rows() == 1);
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateDerivative(BasisFunction<1>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {      
        case EDGE2:
        case EDGE3:
        {
          basisFunction.resizeShapeFunctionDeriv(2,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
             shapeDerivMatrix(0)(0,j) = -0.5;
             shapeDerivMatrix(0)(1,j) = .5;
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {      
        // Lagrange linear shape function derivatives
        case EDGE3:
        {
          basisFunction.resizeShapeFunctionDeriv(3,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
	    shapeDerivMatrix(0)(0,j) = p1-.5;
	    shapeDerivMatrix(0)(1,j) = p1+.5;
	    shapeDerivMatrix(0)(2,j) = -2.*p1;
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis Service: Unsupported order "<<order;
      std::logic_error("halt");
    }
  }
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateDerivative(BasisFunction<2>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {      
        case TRI3:
        case TRI6:
        {
          basisFunction.resizeShapeFunctionDeriv(3,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            shapeDerivMatrix(0)(0,j) = -.5;
            shapeDerivMatrix(0)(1,j) = .5;
            shapeDerivMatrix(0)(2,j) = 0. ;
            shapeDerivMatrix(1)(0,j) = -.5;
            shapeDerivMatrix(1)(1,j) = 0.;
            shapeDerivMatrix(1)(2,j) = 0.5;
          }
	  return;
        }
        case QUAD4:
        case QUAD8:
        case QUAD9:
        {
          basisFunction.resizeShapeFunctionDeriv(4,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeDerivMatrix(0)(0,j) = -0.25 * (1-p2);
            shapeDerivMatrix(0)(1,j) = 0.25 * (1-p2);
            shapeDerivMatrix(0)(2,j) = 0.25 * (1+p2);
            shapeDerivMatrix(0)(3,j) = -0.25 * (1+p2);
            shapeDerivMatrix(1)(0,j) = -0.25 * (1-p1);
            shapeDerivMatrix(1)(1,j) = -0.25 * (1+p1);
            shapeDerivMatrix(1)(2,j) = 0.25 * (1+p1);
            shapeDerivMatrix(1)(3,j) = 0.25 * (1-p1);
	  }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {      
        case TRI6:
        {
          basisFunction.resizeShapeFunctionDeriv(6,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeDerivMatrix(0)(0,j) = p1 + p2 + 0.5;
            shapeDerivMatrix(0)(1,j) = p1 + .5;
            shapeDerivMatrix(0)(2,j) = 0.;
            shapeDerivMatrix(0)(3,j) = - 2. * p1 - p2 - 1.;
            shapeDerivMatrix(0)(4,j) = p2 + 1.;
            shapeDerivMatrix(0)(5,j) = -p2 - 1.;
            shapeDerivMatrix(1)(0,j) = p1 + p2 + 0.5;
            shapeDerivMatrix(1)(1,j) = 0.;
            shapeDerivMatrix(1)(2,j) = p2 + .5;
            shapeDerivMatrix(1)(3,j) = - p1 - 1.;
            shapeDerivMatrix(1)(4,j) = p1 + 1.;
            shapeDerivMatrix(1)(5,j) = -2.*p2 - p1 - 1.;
	  }
	  return;
        }
        case QUAD9:
        {
          basisFunction.resizeShapeFunctionDeriv(9,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            shapeDerivMatrix(0)(0,j) = .25 * (2.*p1 - 1.) * p2 * (p2 - 1.);
            shapeDerivMatrix(0)(1,j) = .25 * (2.*p1 + 1.) * p2 * (p2 - 1.);
            shapeDerivMatrix(0)(2,j) = .25 * (2.*p1 + 1.) * p2 * (p2 + 1.);
            shapeDerivMatrix(0)(3,j) = .25 * (2.*p1 - 1.) * p2 * (p2 + 1.);
            shapeDerivMatrix(0)(4,j) = .5 * (- 2.*p1) * p2 * (p2 - 1.);
            shapeDerivMatrix(0)(5,j) = .5 * (2.*p1 + 1.) * (1. - p2*p2);
            shapeDerivMatrix(0)(6,j) = .5 * (- 2.*p1) * p2 * (p2 + 1.);
            shapeDerivMatrix(0)(7,j) = .5 * (2.*p1 - 1.) * (1. - p2*p2);
            shapeDerivMatrix(0)(8,j) = (- 2.*p1) * (1. - p2*p2);
            shapeDerivMatrix(1)(0,j) = .25 * p1 * (p1 - 1.) * (2.*p2 - 1.);
            shapeDerivMatrix(1)(1,j) = .25 * p1 * (p1 + 1.) * (2.*p2 - 1.);
            shapeDerivMatrix(1)(2,j) = .25 * p1 * (p1 + 1.) * (2.*p2 + 1.);
            shapeDerivMatrix(1)(3,j) = .25 * p1 * (p1 - 1.) * (2.*p2 + 1.);
            shapeDerivMatrix(1)(4,j) = .5 * (1. - p1*p1) * (2.*p2 - 1.);
            shapeDerivMatrix(1)(5,j) = .5 * p1 * (p1 + 1.) * (- 2.*p2);
            shapeDerivMatrix(1)(6,j) = .5 * (1. - p1*p1) * (2.*p2 + 1.);
            shapeDerivMatrix(1)(7,j) = .5 * p1 * (p1 - 1.) * (- 2.*p2);
            shapeDerivMatrix(1)(8,j) = (1. - p1*p1) * (- 2.*p2);
          }
	  return;
        }
        case QUAD8:
        {
          basisFunction.resizeShapeFunctionDeriv(8,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
	    shapeDerivMatrix(0)(0,j) = .25*(1. - p2)*((1. - p1)*(-1.) + (-1.)*(-1. - p1 - p2));
	    shapeDerivMatrix(0)(1,j) = .25*(1. - p2)*((1. + p1)*(1.) + (1.)*(-1. + p1 - p2));
	    shapeDerivMatrix(0)(2,j) = .25*(1. + p2)*((1. + p1)*(1.) + (1.)*(-1. + p1 + p2));
	    shapeDerivMatrix(0)(3,j) = .25*(1. + p2)*((1. - p1)*(-1.) + (-1.)*(-1. - p1 + p2));
	    shapeDerivMatrix(0)(4,j) = .5*(-2.*p1)*(1. - p2);
	    shapeDerivMatrix(0)(5,j) = .5*(1.)*(1. - p2*p2);
	    shapeDerivMatrix(0)(6,j) = .5*(-2.*p1)*(1. + p2);
	    shapeDerivMatrix(0)(7,j) = .5*(-1.)*(1. - p2*p2);
            shapeDerivMatrix(1)(0,j) = .25*(1. - p1)*((1. - p2)*(-1.) + (-1.)*(-1. - p1 - p2));
            shapeDerivMatrix(1)(1,j) = .25*(1. + p1)*((1. - p2)*(-1.) + (-1.)*(-1. + p1 - p2));
            shapeDerivMatrix(1)(2,j) = .25*(1. + p1)*((1. + p2)*(1.) + (1.)*(-1. + p1 + p2));
            shapeDerivMatrix(1)(3,j) = .25*(1. - p1)*((1. + p2)*(1.) + (1.)*(-1. - p1 + p2));
            shapeDerivMatrix(1)(4,j) = .5*(1. - p1*p1)*(-1.);
            shapeDerivMatrix(1)(5,j) = .5*(1. + p1)*(-2.*p2);
            shapeDerivMatrix(1)(6,j) = .5*(1. - p1*p1)*(1.);
            shapeDerivMatrix(1)(7,j) = .5*(1. - p1)*(-2.*p2);
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis Service: Unsupported order "<<order;
      std::logic_error("halt");
    }
  }
}

void
ExplicitLagrangianFiniteElementBasis::
evaluateDerivative(BasisFunction<3>& basisFunction, const ElemType itemType, const unsigned int order, const PointVector& points)
{
  BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix = basisFunction.getShapeFunctionDeriv();
  switch (order)
  {
    case 1:
    {
      switch (itemType)
      {      
        case TET4:
        case TET10:
        {
          basisFunction.resizeShapeFunctionDeriv(4,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            shapeDerivMatrix(0)(0,j) = -0.5;
            shapeDerivMatrix(0)(1,j) = .5;
            shapeDerivMatrix(0)(2,j) = 0.;
            shapeDerivMatrix(0)(3,j) = 0.;
            shapeDerivMatrix(1)(0,j) = -0.5;
            shapeDerivMatrix(1)(1,j) = 0.;
            shapeDerivMatrix(1)(2,j) = 0.5;
            shapeDerivMatrix(1)(3,j) = 0.;
            shapeDerivMatrix(2)(0,j) = -0.5;
            shapeDerivMatrix(2)(1,j) = 0.;
            shapeDerivMatrix(2)(2,j) = 0.;
            shapeDerivMatrix(2)(3,j) = 0.5;
          }
	  return;
        }
        case HEX8:
        case HEX20:
        case HEX27:
        {
          basisFunction.resizeShapeFunctionDeriv(8,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeDerivMatrix(0)(0,j) = -0.125*(1. - p2)*(1. - p3);
            shapeDerivMatrix(0)(1,j) = .125*(1. - p2)*(1. - p3);
            shapeDerivMatrix(0)(2,j) = .125*(1. + p2)*(1. - p3);
            shapeDerivMatrix(0)(3,j) = -0.125*(1. + p2)*(1. - p3);
            shapeDerivMatrix(0)(4,j) = -0.125*(1. - p2)*(1. + p3);
            shapeDerivMatrix(0)(5,j) = .125*(1. - p2)*(1. + p3);
            shapeDerivMatrix(0)(6,j) = .125*(1. + p2)*(1. + p3);
            shapeDerivMatrix(0)(7,j) = -0.125*(1. + p2)*(1. + p3);
            shapeDerivMatrix(1)(0,j) = -0.125*(1. - p1)*(1. - p3);
            shapeDerivMatrix(1)(1,j) = -0.125*(1. + p1)*(1. - p3);
            shapeDerivMatrix(1)(2,j) = .125*(1. + p1)*(1. - p3);
            shapeDerivMatrix(1)(3,j) = .125*(1. - p1)*(1. - p3);
            shapeDerivMatrix(1)(4,j) = -0.125*(1. - p1)*(1. + p3);
            shapeDerivMatrix(1)(5,j) = -0.125*(1. + p1)*(1. + p3);
            shapeDerivMatrix(1)(6,j) = .125*(1. + p1)*(1. + p3);
            shapeDerivMatrix(1)(7,j) = .125*(1. - p1)*(1. + p3);
            shapeDerivMatrix(2)(0,j) = -0.125*(1. - p1)*(1. - p2);
            shapeDerivMatrix(2)(1,j) = -0.125*(1. + p1)*(1. - p2);
            shapeDerivMatrix(2)(2,j) = -0.125*(1. + p1)*(1. + p2);
            shapeDerivMatrix(2)(3,j) = -0.125*(1. - p1)*(1. + p2);
            shapeDerivMatrix(2)(4,j) = .125*(1. - p1)*(1. - p2);
            shapeDerivMatrix(2)(5,j) = .125*(1. + p1)*(1. - p2);
            shapeDerivMatrix(2)(6,j) = .125*(1. + p1)*(1. + p2);
            shapeDerivMatrix(2)(7,j) = .125*(1. - p1)*(1. + p2);
	  }
	  return;
        }
        case PRISM6:
        case PRISM15:
        {
          basisFunction.resizeShapeFunctionDeriv(6,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
            shapeDerivMatrix(0)(0,j) = -0.25*(1.- p3);
            shapeDerivMatrix(0)(1,j) = .25*(1. - p3);
            shapeDerivMatrix(0)(2,j) = .0;
            shapeDerivMatrix(0)(3,j) = -0.25*(1.+ p3);
            shapeDerivMatrix(0)(4,j) = .25*(1. + p3);
            shapeDerivMatrix(0)(5,j) = .0;
            shapeDerivMatrix(1)(0,j) = -0.25*(1.- p3);
            shapeDerivMatrix(1)(1,j) = .0;
            shapeDerivMatrix(1)(2,j) = .25*(1. - p3);
            shapeDerivMatrix(1)(3,j) = -0.25*(1.+ p3);
            shapeDerivMatrix(1)(4,j) = .0;
            shapeDerivMatrix(1)(5,j) = .25*(1. + p3);
            shapeDerivMatrix(2)(0,j) = .25*(+p1 + p2);
            shapeDerivMatrix(2)(1,j) = .25*(-p1 - 1.);
            shapeDerivMatrix(2)(2,j) = .25*(-p2 - 1.);
            shapeDerivMatrix(2)(3,j) = .25*(-p1 - p2);
            shapeDerivMatrix(2)(4,j) = .25*(p1 + 1.);
            shapeDerivMatrix(2)(5,j) = .25*(p2 + 1.);
	  }
	  return;
        }
        case PYRAMID5:
        case PYRAMID13:
        {
          basisFunction.resizeShapeFunctionDeriv(5,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
            double p1 = points(j,0);
            double p2 = points(j,1);
            double p3 = points(j,2);
	    double oneMp3 = 1.-p3+1.e-30;
	    double oneMp3_sq = (1.-p3)*(1.-p3)+1.e-30;
            shapeDerivMatrix(0)(0,j) =  0.5 * (p2 + p3) / oneMp3;
            shapeDerivMatrix(0)(1,j) = -0.5 * (p2 + p3) / oneMp3;
            shapeDerivMatrix(0)(2,j) =  0.5 * (1. + p2) / oneMp3;
            shapeDerivMatrix(0)(3,j) = -0.5 * (1. + p2) / oneMp3;  
            shapeDerivMatrix(0)(4,j) =  0.; 
            shapeDerivMatrix(1)(0,j) =  0.5 * (p1 + p3) / oneMp3;
            shapeDerivMatrix(1)(1,j) = -0.5 * (1. + p1) / oneMp3;
            shapeDerivMatrix(1)(2,j) =  0.5 * (1. + p1) / oneMp3;
            shapeDerivMatrix(1)(3,j) = -0.5 * (p1 + p3) / oneMp3;
            shapeDerivMatrix(1)(4,j) =  0.; 
            shapeDerivMatrix(2)(0,j) = -0.5 * (-p1 -p3) * (-p2 -p3) / oneMp3_sq;
            shapeDerivMatrix(2)(1,j) = -0.5 * (1. + p1) * (-p2 -p3) / oneMp3_sq;
            shapeDerivMatrix(2)(2,j) = -0.5 * (1. + p1) * (1. + p2) / oneMp3_sq;
            shapeDerivMatrix(2)(3,j) = -0.5 * (-p1 -p3) * (1. + p2) / oneMp3_sq;
            shapeDerivMatrix(2)(4,j) =  0.5; 
	  }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    case 2:
    {
      switch (itemType)
      {      
        case TET10:
        {
          basisFunction.resizeShapeFunctionDeriv(10,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
	    const double zeta1 = 0.5*(points(j,0)+1.);
	    const double zeta2 = 0.5*(points(j,1)+1.);
	    const double zeta3 = 0.5*(points(j,2)+1.);
	    const double zeta0 = 1. - zeta1 - zeta2 - zeta3;
	    shapeDerivMatrix(0)(0,j) = -2.*zeta0 + 0.5;
	    shapeDerivMatrix(0)(1,j) = 2.*zeta1 - 0.5;
	    shapeDerivMatrix(0)(2,j) = 0.;
	    shapeDerivMatrix(0)(3,j) = 0.;
	    shapeDerivMatrix(0)(4,j) = 2.*(zeta0 - zeta1);
	    shapeDerivMatrix(0)(5,j) = 2.*zeta2;
	    shapeDerivMatrix(0)(6,j) = -2.*zeta2;
	    shapeDerivMatrix(0)(7,j) = -2.*zeta3;
	    shapeDerivMatrix(0)(8,j) = 2.*zeta3;
	    shapeDerivMatrix(0)(9,j) = 0.;
	    shapeDerivMatrix(1)(0,j) = -2.*zeta0 + 0.5;
	    shapeDerivMatrix(1)(1,j) = 0.;
	    shapeDerivMatrix(1)(2,j) = 2.*zeta2 - 0.5;
	    shapeDerivMatrix(1)(3,j) = 0.;
	    shapeDerivMatrix(1)(4,j) = -2.*zeta1;
	    shapeDerivMatrix(1)(5,j) = 2.*zeta1;
	    shapeDerivMatrix(1)(6,j) = 2.*(zeta0 - zeta2);
	    shapeDerivMatrix(1)(7,j) = -2.* zeta3;
	    shapeDerivMatrix(1)(8,j) = 0.;
	    shapeDerivMatrix(1)(9,j) = 2.*zeta3;
	    shapeDerivMatrix(2)(0,j) = -2.*zeta0 + 0.5;
	    shapeDerivMatrix(2)(1,j) = 0.;
	    shapeDerivMatrix(2)(2,j) = 0.;
	    shapeDerivMatrix(2)(3,j) = 2.*zeta3 - 0.5;
	    shapeDerivMatrix(2)(4,j) = -2.*zeta1;
	    shapeDerivMatrix(2)(5,j) = 0.;
	    shapeDerivMatrix(2)(6,j) = -2.*zeta2;
	    shapeDerivMatrix(2)(7,j) = 2.*(zeta0 - zeta3);
	    shapeDerivMatrix(2)(8,j) = 2.*zeta1;
	    shapeDerivMatrix(2)(9,j) = 2.*zeta2;
          }
	  return;
        }
        case PYRAMID13:
        {
          basisFunction.resizeShapeFunctionDeriv(13,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
	    const double p1 = points(j,0);
	    const double p2 = points(j,1);
	    const double p3 = points(j,2);
	    double oneMp3 = 1.-p3+1.e-30;
	    double oneMp3_sq = (1.-p3)*(1.-p3)+1.e-30;
	    double oneMp3_th = (1.-p3)*(1.-p3)*(1.-p3)+1.e-30;
            shapeDerivMatrix(0)(0,j) = -0.5*( (2.-p3)*(-p2-p3)/oneMp3 + (-4.*(1.+p1)-2.*(1.+p2))*(-p2-p3)/oneMp3_sq);
            shapeDerivMatrix(0)(1,j) =  0.5*((-2.-p3)*(-p2-p3)/oneMp3 + ( 4.*(1.+p1)-2.*(1.+p2))*(-p2-p3)/oneMp3_sq);
            shapeDerivMatrix(0)(2,j) =  0.5*((-4.-p3)*(1.+p2) /oneMp3 + ( 4.*(1.+p1)+2.*(1.+p2))*(1.+p2) /oneMp3_sq);
            shapeDerivMatrix(0)(3,j) = -0.5*(    -p3 *(1.+p2) /oneMp3 + (-4.*(1.+p1)+2.*(1.+p2))*(1.+p2) /oneMp3_sq);
	    shapeDerivMatrix(0)(4,j) = 0.;
            shapeDerivMatrix(0)(5,j) = 2.*((-p2-p3)/oneMp3 - 2*(1+p1)*(-p2-p3)/oneMp3_sq);
            shapeDerivMatrix(0)(6,j) = -2.*(1+p2)*(p2+p3)/oneMp3_sq;
            shapeDerivMatrix(0)(7,j) = 2.*((1.+p2) /oneMp3 - 2*(1+p1)*(1.+p2) /oneMp3_sq);
            shapeDerivMatrix(0)(8,j) =  2.*(1+p2)*(p2+p3)/oneMp3_sq;
            shapeDerivMatrix(0)(9,j)  = -(-p2-p3)*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(0)(10,j) =  (-p2-p3)*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(0)(11,j) =  (1.+p2 )*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(0)(12,j) = -(1.+p2 )*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(1)(0,j) = -0.5*( (2.-p3)*(-p1-p3)/oneMp3 + (-4.*(1.+p2)-2.*(1.+p1))*(-p1-p3)/oneMp3_sq);
            shapeDerivMatrix(1)(1,j) = -0.5*(    -p3 *(1.+p1) /oneMp3 + (-4.*(1.+p2)+2.*(1.+p1))*(1.+p1) /oneMp3_sq);
            shapeDerivMatrix(1)(2,j) =  0.5*((-4.-p3)*(1.+p1) /oneMp3 + ( 4.*(1.+p2)+2.*(1.+p1))*(1.+p1) /oneMp3_sq);
            shapeDerivMatrix(1)(3,j) =  0.5*((-2.-p3)*(-p1-p3)/oneMp3 + ( 4.*(1.+p2)-2.*(1.+p1))*(-p1-p3)/oneMp3_sq);
	    shapeDerivMatrix(1)(4,j) = 0.;
            shapeDerivMatrix(1)(5,j) =  2.*(1+p1)*(p1+p3)/oneMp3_sq;
            shapeDerivMatrix(1)(6,j) = 2.*((1.+p1) /oneMp3 - 2*(1+p2)*(1.+p1) /oneMp3_sq);
            shapeDerivMatrix(1)(7,j) = -2.*(1+p1)*(p1+p3)/oneMp3_sq;
            shapeDerivMatrix(1)(8,j) = 2.*((-p1-p3)/oneMp3 - 2*(1+p2)*(-p1-p3)/oneMp3_sq);
            shapeDerivMatrix(1)(9,j)  = -(-p1-p3)*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(1)(10,j) = -(1.+p1) *(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(1)(11,j) =  (1.+p1) *(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(1)(12,j) =  (-p1-p3)*(1.-p3*p3)/oneMp3_sq;
            shapeDerivMatrix(2)(0,j) = -0.5*((-2.*p3+1)*(-p1-p3)*(-p2-p3)/oneMp3_sq + (-2.*(1.+p1)-2.*(1.+p2))*(-p1-p3)*(-p2-p3)/oneMp3_th);
            shapeDerivMatrix(2)(1,j) = -0.5*((-2.*p3-1)*(1.+p1) *(-p2-p3)/oneMp3_sq + ( 2.*(1.+p1)-2.*(1.+p2))*(1.+p1) *(-p2-p3)/oneMp3_th);
            shapeDerivMatrix(2)(2,j) = -0.5*((-2.*p3-3)*(1.+p1) *(1.+p2) /oneMp3_sq + ( 2.*(1.+p1)+2.*(1.+p2))*(1.+p1) *(1.+p2) /oneMp3_th);
            shapeDerivMatrix(2)(3,j) = -0.5*((-2.*p3-1)*(-p1-p3)*(1.+p2) /oneMp3_sq + (-2.*(1.+p1)+2.*(1.+p2))*(-p1-p3)*(1.+p2) /oneMp3_th);
	    shapeDerivMatrix(2)(4,j) = 0.5 + p3;
            shapeDerivMatrix(2)(5,j) = 2.*((-p2-p3)*(1+p1)*(p1+p3))/oneMp3_th;
            shapeDerivMatrix(2)(6,j) = 2.*((1.+p1) *(1+p2)*(p2+p3))/oneMp3_th;
            shapeDerivMatrix(2)(7,j) = 2.*((1.+p2) *(1+p1)*(p1+p3))/oneMp3_th;
            shapeDerivMatrix(2)(8,j) = 2.*((-p1-p3)*(1+p2)*(p2+p3))/oneMp3_th;
            shapeDerivMatrix(2)(9,j)  = -2.*(-p1 -p3) * (-p2 -p3) * p3 / oneMp3_sq;
            shapeDerivMatrix(2)(10,j) = -2.*(1. + p1) * (-p2 -p3) * p3 / oneMp3_sq;
            shapeDerivMatrix(2)(11,j) = -2.*(1. + p1) * (1. + p2) * p3 / oneMp3_sq;
            shapeDerivMatrix(2)(12,j) = -2.*(-p1 -p3) * (1. + p2) * p3 / oneMp3_sq;
          }
	  return;
        }
        case PRISM15:
        {
          basisFunction.resizeShapeFunctionDeriv(15,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
	    const double p1 = points(j,0);
	    const double p2 = points(j,1);
	    const double p3 = points(j,2);
            shapeDerivMatrix(0)(0,j) = 0.25 * (1. - p3) * (2. + 2.*p1 +2.*p2 + p3);
            shapeDerivMatrix(0)(1,j) = 0.25 * (1. - p3) * (2.*p1 - p3);
	    shapeDerivMatrix(0)(2,j) = 0.;
            shapeDerivMatrix(0)(3,j) = 0.25 * (1. + p3) * (2. + 2.*p1 +2.*p2 - p3);
            shapeDerivMatrix(0)(4,j) = 0.25 * (1. + p3) * (2.*p1 + p3);
            shapeDerivMatrix(0)(5,j) = 0.;
	    shapeDerivMatrix(0)(6,j) = 0.5 * (1. - p3) * (-1. - 2.*p1 - p2);
	    shapeDerivMatrix(0)(7,j) = 0.5 * (1. + p2) * (1 - p3);
	    shapeDerivMatrix(0)(8,j) = -0.5 * (1. + p2) * (1 - p3);
	    shapeDerivMatrix(0)(9,j) = -0.5 * (1. + p3) * (1 - p3);
	    shapeDerivMatrix(0)(10,j) = 0.5 * (1. + p3) * (1 - p3);
	    shapeDerivMatrix(0)(11,j) = 0.;
	    shapeDerivMatrix(0)(12,j) = 0.5 * (1. + p3) * (-1 -2.*p1 - p2);
	    shapeDerivMatrix(0)(13,j) = 0.5 * (1. + p2) * (1 + p3);
	    shapeDerivMatrix(0)(14,j) = -0.5 * (1. + p2) * (1 + p3);
            shapeDerivMatrix(1)(0,j) = 0.25 * (1. - p3) * (2. + 2.*p1 +2.*p2 + p3);
            shapeDerivMatrix(1)(1,j) = 0.;
            shapeDerivMatrix(1)(2,j) = 0.25 * (1. - p3) * (2.*p2 - p3);
            shapeDerivMatrix(1)(3,j) = 0.25 * (1. + p3) * (2. + 2.*p1 +2.*p2 - p3);
            shapeDerivMatrix(1)(4,j) = 0.;
            shapeDerivMatrix(1)(5,j) = 0.25 * (1. + p3) * (2.*p2 + p3);
            shapeDerivMatrix(1)(6,j) = -0.5 * (1. + p1) * (1. - p3);
	    shapeDerivMatrix(1)(7,j) = 0.5 * (1. + p1) * (1 - p3);
            shapeDerivMatrix(1)(8,j) = 0.5 * (1. - p3) * (-1. - p1 - 2.*p2);
	    shapeDerivMatrix(1)(9,j) = -0.5 * (1. + p3) * (1 - p3);
            shapeDerivMatrix(1)(10,j) = 0.;
	    shapeDerivMatrix(1)(11,j) = 0.5 * (1. + p3) * (1 - p3);
	    shapeDerivMatrix(1)(12,j) = -0.5 * (1. + p1) * (1 + p3);
	    shapeDerivMatrix(1)(13,j) = 0.5 * (1. + p1) * (1 + p3);
            shapeDerivMatrix(1)(14,j) = 0.5 * (1. + p3) * (-1. - p1 - 2.*p2);
            shapeDerivMatrix(2)(0,j) = 0.25 * (-p1 - p2) * (1. + p1 + p2 + 2.*p3);
            shapeDerivMatrix(2)(1,j) = 0.25 * (1. + p1) * (2.*p3 - p1);
            shapeDerivMatrix(2)(2,j) = 0.25 * (1. + p2) * (2.*p3 - p2);
            shapeDerivMatrix(2)(3,j) = 0.25 * (-p1 - p2) * (-1. - p1 - p2 + 2.*p3);
            shapeDerivMatrix(2)(4,j) = 0.25 * (1. + p1) * (2.*p3 + p1);
            shapeDerivMatrix(2)(5,j) = 0.25 * (1. + p2) * (2.*p3 + p2);
            shapeDerivMatrix(2)(6,j) = -0.5 * (1. + p1) * (-p1 - p2);
	    shapeDerivMatrix(2)(7,j) = -0.5 * (1. + p1) * (1 + p2);
	    shapeDerivMatrix(2)(8,j) = -0.5 * (1. + p2) * (-p1 - p2);
	    shapeDerivMatrix(2)(9,j) = 0.5 * (-p1 - p2) * (- 2.*p3);
	    shapeDerivMatrix(2)(10,j) = 0.5 * (1. + p1) * (- 2.*p3);
	    shapeDerivMatrix(2)(11,j) = 0.5 * (1. + p2) * (- 2.*p3);
            shapeDerivMatrix(2)(12,j) = 0.5 * (1. + p1) * (-p1 - p2);
	    shapeDerivMatrix(2)(13,j) = 0.5 * (1. + p1) * (1 + p2);
	    shapeDerivMatrix(2)(14,j) = 0.5 * (1. + p2) * (-p1 - p2);
          }
	  return;
        }
        case HEX20:
        {
          basisFunction.resizeShapeFunctionDeriv(20,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
	    const double x = .5*(points(j,0) + 1.);
	    const double y = .5*(points(j,1) + 1.);
	    const double z = .5*(points(j,2) + 1.);
	    shapeDerivMatrix(0)(0,j) = .5*(1. - y)*(1. - z)*((1. - x)*(-2.) + (-1.)*(1. - 2.*x - 2.*y - 2.*z));
	    shapeDerivMatrix(0)(1,j) = .5*(1. - y)*(1. - z)*(x*(2.) + (1.)*(2.*x - 2.*y - 2.*z - 1.));
	    shapeDerivMatrix(0)(2,j) = .5*y*(1. - z)*(x*(2.) + (1.)*(2.*x + 2.*y - 2.*z - 3.));
	    shapeDerivMatrix(0)(3,j) = .5*y*(1. - z)*((1. - x)*(-2.) + (-1.)*(2.*y - 2.*x - 2.*z - 1.));
	    shapeDerivMatrix(0)(4,j) = .5*(1. - y)*z*((1. - x)*(-2.) + (-1.)*(2.*z - 2.*x - 2.*y - 1.));
	    shapeDerivMatrix(0)(5,j) = .5*(1. - y)*z*(x*(2.) + (1.)*(2.*x - 2.*y + 2.*z - 3.));
	    shapeDerivMatrix(0)(6,j) = .5*y*z*(x*(2.) + (1.)*(2.*x + 2.*y + 2.*z - 5.));
	    shapeDerivMatrix(0)(7,j) = .5*y*z*((1. - x)*(-2.) + (-1.)*(2.*y - 2.*x + 2.*z - 3.));
	    shapeDerivMatrix(0)(8,j) = 2.*(1. - y)*(1. - z)*(1. - 2.*x);
	    shapeDerivMatrix(0)(9,j) = 2.*y*(1. - y)*(1. - z);
	    shapeDerivMatrix(0)(10,j) =	2.*y*(1. - z)*(1. - 2.*x);
	    shapeDerivMatrix(0)(11,j) =	2.*y*(1. - y)*(1. - z)*(-1.);
	    shapeDerivMatrix(0)(12,j) =	2.*(1. - y)*z*(1. - z)*(-1.);
	    shapeDerivMatrix(0)(13,j) =	2.*(1. - y)*z*(1. - z);
	    shapeDerivMatrix(0)(14,j) =	2.*y*z*(1. - z);
	    shapeDerivMatrix(0)(15,j) = 2.*y*z*(1. - z)*(-1.);
	    shapeDerivMatrix(0)(16,j) = 2.*(1. - y)*z*(1. - 2.*x);
	    shapeDerivMatrix(0)(17,j) = 2.*y*(1. - y)*z;
	    shapeDerivMatrix(0)(18,j) = 2.*y*z*(1. - 2.*x);
	    shapeDerivMatrix(0)(19,j) = 2.*y*(1. - y)*z*(-1.);
	    shapeDerivMatrix(1)(0,j) = .5*(1. - x)*(1. - z)*((1. - y)*(-2.) + (-1.)*(1. - 2.*x - 2.*y - 2.*z));
	    shapeDerivMatrix(1)(1,j) = .5*x*(1. - z)*((1. - y)*(-2.) + (-1.)*(2.*x - 2.*y - 2.*z - 1.));
	    shapeDerivMatrix(1)(2,j) = .5*x*(1. - z)*(y*(2.) + (1.)*(2.*x + 2.*y - 2.*z - 3.));
	    shapeDerivMatrix(1)(3,j) = .5*(1. - x)*(1. - z)*(y*(2.) + (1.)*(2.*y - 2.*x - 2.*z - 1.));
	    shapeDerivMatrix(1)(4,j) = .5*(1. - x)*z*((1. - y)*(-2.) + (-1.)*(2.*z - 2.*x - 2.*y - 1.));
	    shapeDerivMatrix(1)(5,j) = .5*x*z*((1. - y)*(-2.) + (-1.)*(2.*x - 2.*y + 2.*z - 3.));
	    shapeDerivMatrix(1)(6,j) = .5*x*z*(y*(2.) + (1.)*(2.*x + 2.*y + 2.*z - 5.));
	    shapeDerivMatrix(1)(7,j) = .5*(1. - x)*z*(y*(2.) + (1.)*(2.*y - 2.*x + 2.*z - 3.));
	    shapeDerivMatrix(1)(8,j) = 2.*x*(1. - x)*(1. - z)*(-1.);
	    shapeDerivMatrix(1)(9,j) = 2.*x*(1. - z)*(1. - 2.*y);
	    shapeDerivMatrix(1)(10,j) = 2.*x*(1. - x)*(1. - z);
	    shapeDerivMatrix(1)(11,j) = 2.*(1. - x)*(1. - z)*(1. - 2.*y);
	    shapeDerivMatrix(1)(12,j) = 2.*(1. - x)*z*(1. - z)*(-1.);
	    shapeDerivMatrix(1)(13,j) = 2.*x*z*(1. - z)*(-1.);
	    shapeDerivMatrix(1)(14,j) = 2.*x*z*(1. - z);
	    shapeDerivMatrix(1)(15,j) = 2.*(1. - x)*z*(1. - z);
	    shapeDerivMatrix(1)(16,j) = 2.*x*(1. - x)*z*(-1.);
	    shapeDerivMatrix(1)(17,j) = 2.*x*z*(1. - 2.*y);
	    shapeDerivMatrix(1)(18,j) = 2.*x*(1. - x)*z;
	    shapeDerivMatrix(1)(19,j) = 2.*(1. - x)*z*(1. - 2.*y);
	    shapeDerivMatrix(2)(0,j) = .5*(1. - x)*(1. - y)*((1. - z)*(-2.) + (-1.)*(1. - 2.*x - 2.*y - 2.*z));
	    shapeDerivMatrix(2)(1,j) = .5*x*(1. - y)*((1. - z)*(-2.) + (-1.)*(2.*x - 2.*y - 2.*z - 1.));
	    shapeDerivMatrix(2)(2,j) = .5*x*y*((1. - z)*(-2.) + (-1.)*(2.*x + 2.*y - 2.*z - 3.));
	    shapeDerivMatrix(2)(3,j) = .5*(1. - x)*y*((1. - z)*(-2.) + (-1.)*(2.*y - 2.*x - 2.*z - 1.));
	    shapeDerivMatrix(2)(4,j) = .5*(1. - x)*(1. - y)*(z*(2.) + (1.)*(2.*z - 2.*x - 2.*y - 1.));
	    shapeDerivMatrix(2)(5,j) = .5*x*(1. - y)*(z*(2.) + (1.)*(2.*x - 2.*y + 2.*z - 3.));
	    shapeDerivMatrix(2)(6,j) = .5*x*y*(z*(2.) + (1.)*(2.*x + 2.*y + 2.*z - 5.));
	    shapeDerivMatrix(2)(7,j) = .5*(1. - x)*y*(z*(2.) + (1.)*(2.*y - 2.*x + 2.*z - 3.));
	    shapeDerivMatrix(2)(8,j) = 2.*x*(1. - x)*(1. - y)*(-1.);
	    shapeDerivMatrix(2)(9,j) = 2.*x*y*(1. - y)*(-1.);
	    shapeDerivMatrix(2)(10,j) = 2.*x*(1. - x)*y*(-1.);
	    shapeDerivMatrix(2)(11,j) = 2.*(1. - x)*y*(1. - y)*(-1.);
	    shapeDerivMatrix(2)(12,j) = 2.*(1. - x)*(1. - y)*(1. - 2.*z);
	    shapeDerivMatrix(2)(13,j) = 2.*x*(1. - y)*(1. - 2.*z);
	    shapeDerivMatrix(2)(14,j) = 2.*x*y*(1. - 2.*z);
	    shapeDerivMatrix(2)(15,j) = 2.*(1. - x)*y*(1. - 2.*z);
	    shapeDerivMatrix(2)(16,j) = 2.*x*(1. - x)*(1. - y);
	    shapeDerivMatrix(2)(17,j) = 2.*x*y*(1. - y);
	    shapeDerivMatrix(2)(18,j) = 2.*x*(1. - x)*y;
	    shapeDerivMatrix(2)(19,j) = 2.*(1. - x)*y*(1. - y);
          }
	  return;
        }
        case HEX27:
        {
          basisFunction.resizeShapeFunctionDeriv(27,points.rows());
          for (unsigned int j = 0; j < points.rows(); j++)
          {
	    const double p1 = points(j,0);
	    const double p2 = points(j,1);
	    const double p3 = points(j,2);
            shapeDerivMatrix(0)(0,j) = .125 * (2.*p1 - 1.) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(1,j) = .125 * (2.*p1 + 1.) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(2,j) = .125 * (2.*p1 + 1.) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(3,j) = .125 * (2.*p1 - 1.) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(4,j) = .125 * (2.*p1 - 1.) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(5,j) = .125 * (2.*p1 + 1.) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(6,j) = .125 * (2.*p1 + 1.) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(7,j) = .125 * (2.*p1 - 1.) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(8,j) = .25 * (- 2.*p1) * p2 * (p2 - 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(9,j) = .25 * (2.*p1 + 1.) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(10,j) = .25 * (- 2.*p1) * p2 * (p2 + 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(11,j) = .25 * (2.*p1 - 1.) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(12,j) = .25 * (2.*p1 - 1.) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(13,j) = .25 * (2.*p1 + 1.) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(14,j) = .25 * (2.*p1 + 1.) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(15,j) = .25 * (2.*p1 - 1.) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(16,j) = .25 * (- 2.*p1) * p2 * (p2 - 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(17,j) = .25 * (2.*p1 + 1.) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(18,j) = .25 * (- 2.*p1) * p2 * (p2 + 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(19,j) = .25 * (2.*p1 - 1.) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(20,j) = .5 * (- 2.*p1) * p2 * (p2 - 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(21,j) = .5 * (2.*p1 + 1.) * (1. - p2*p2) * (1. - p3*p3);
            shapeDerivMatrix(0)(22,j) = .5 * (- 2.*p1) * p2 * (p2 + 1.) * (1. - p3*p3);
            shapeDerivMatrix(0)(23,j) = .5 * (2.*p1 - 1.) * (1. - p2*p2) * (1. - p3*p3);
            shapeDerivMatrix(0)(24,j) = .5 * (- 2.*p1) * (1. - p2*p2) * p3 * (p3 - 1.);
            shapeDerivMatrix(0)(25,j) = .5 * (- 2.*p1) * (1. - p2*p2) * p3 * (p3 + 1.);
            shapeDerivMatrix(0)(26,j) = (- 2.*p1) * (1. - p2*p2) * (1. - p3*p3);
            shapeDerivMatrix(1)(0,j) = .125 * p1 * (p1 - 1.) * (2.*p2 - 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(1)(1,j) = .125 * p1 * (p1 + 1.) * (2.*p2 - 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(1)(2,j) = .125 * p1 * (p1 + 1.) * (2.*p2 + 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(1)(3,j) = .125 * p1 * (p1 - 1.) * (2.*p2 + 1.) * p3 * (p3 - 1.);
            shapeDerivMatrix(1)(4,j) = .125 * p1 * (p1 - 1.) * (2.*p2 - 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(1)(5,j) = .125 * p1 * (p1 + 1.) * (2.*p2 - 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(1)(6,j) = .125 * p1 * (p1 + 1.) * (2.*p2 + 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(1)(7,j) = .125 * p1 * (p1 - 1.) * (2.*p2 + 1.) * p3 * (p3 + 1.);
            shapeDerivMatrix(1)(8,j) = .25 * (1. - p1*p1) * (2.*p2 - 1.) * p3 * (p3 - 1.);   
            shapeDerivMatrix(1)(9,j) = .25 * p1 * (p1 + 1.) * (- 2.*p2) * p3 * (p3 - 1.);    
            shapeDerivMatrix(1)(10,j) = .25 * (1. - p1*p1) * (2.*p2 + 1.) * p3 * (p3 - 1.);   
            shapeDerivMatrix(1)(11,j) = .25 * p1 * (p1 - 1.) * (- 2.*p2) * p3 * (p3 - 1.);    
            shapeDerivMatrix(1)(12,j) = .25 * p1 * (p1 - 1.) * (2.*p2 - 1.) * (1. - p3*p3);   
            shapeDerivMatrix(1)(13,j) = .25 * p1 * (p1 + 1.) * (2.*p2 - 1.) * (1. - p3*p3);   
            shapeDerivMatrix(1)(14,j) = .25 * p1 * (p1 + 1.) * (2.*p2 + 1.) * (1. - p3*p3);   
            shapeDerivMatrix(1)(15,j) = .25 * p1 * (p1 - 1.) * (2.*p2 + 1.) * (1. - p3*p3);   
            shapeDerivMatrix(1)(16,j) = .25 * (1. - p1*p1) * (2.*p2 - 1.) * p3 * (p3 + 1.);   
            shapeDerivMatrix(1)(17,j) = .25 * p1 * (p1 + 1.) * (- 2.*p2) * p3 * (p3 + 1.);    
            shapeDerivMatrix(1)(18,j) = .25 * (1. - p1*p1) * (2.*p2 + 1.) * p3 * (p3 + 1.);   
            shapeDerivMatrix(1)(19,j) = .25 * p1 * (p1 - 1.) * (- 2.*p2) * p3 * (p3 + 1.);    
            shapeDerivMatrix(1)(20,j) = .5 * (1. - p1*p1) * (2.*p2 - 1.) * (1. - p3*p3);      
            shapeDerivMatrix(1)(21,j) = .5 * p1 * (p1 + 1.) * (- 2.*p2) * (1. - p3*p3);       
            shapeDerivMatrix(1)(22,j) = .5 * (1. - p1*p1) * (2.*p2 + 1.) * (1. - p3*p3);      
            shapeDerivMatrix(1)(23,j) = .5 * p1 * (p1 - 1.) * (-2.*p2) * (1. - p3*p3);        
            shapeDerivMatrix(1)(24,j) = .5 * (1. - p1*p1) * (- 2.*p2) * p3 * (p3 - 1);        
            shapeDerivMatrix(1)(25,j) = .5 * (1. - p1*p1) * (-2.*p2) * p3 * (p3 + 1.);         
            shapeDerivMatrix(1)(26,j) = (1. - p1*p1) * (- 2.*p2) * (1. - p3*p3);
            shapeDerivMatrix(2)(0,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(1,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(2,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(3,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(4,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(5,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(6,j) = .125 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(7,j) = .125 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(8,j) = .25 * (1. - p1*p1) * p2 * (p2 - 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(9,j) = .25 * p1 * (p1 + 1.) * (1. - p2*p2) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(10,j) = .25 * (1. - p1*p1) * p2 * (p2 + 1.) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(11,j) = .25 * p1 * (p1 - 1.) * (1. - p2*p2) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(12,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 - 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(13,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 - 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(14,j) = .25 * p1 * (p1 + 1.) * p2 * (p2 + 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(15,j) = .25 * p1 * (p1 - 1.) * p2 * (p2 + 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(16,j) = .25 * (1. - p1*p1) * p2 * (p2 - 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(17,j) = .25 * p1 * (p1 + 1.) * (1. - p2*p2) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(18,j) = .25 * (1. - p1*p1) * p2 * (p2 + 1.) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(19,j) = .25 * p1 * (p1 - 1.) * (1. - p2*p2) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(20,j) = .5 * (1. - p1*p1) * p2 * (p2 - 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(21,j) = .5 * p1 * (p1 + 1.) * (1. - p2*p2) * (- 2.*p3);
            shapeDerivMatrix(2)(22,j) = .5 * (1. - p1*p1) * p2 * (p2 + 1.) * (- 2.*p3);
            shapeDerivMatrix(2)(23,j) = .5 * p1 * (p1 - 1.) * (1. - p2*p2) * (- 2.*p3);
            shapeDerivMatrix(2)(24,j) = .5 * (1. - p1*p1) * (1. - p2*p2) * (2.*p3 - 1.);
            shapeDerivMatrix(2)(25,j) = .5 * (1. - p1*p1) * (1. - p2*p2) * (2.*p3 + 1.);
            shapeDerivMatrix(2)(26,j) = (1. - p1*p1) * (1. - p2*p2) * (- 2.*p3);
          }
	  return;
        }
        default:
        {
          std::cerr<< "ExplicitLagrangian Basis Service: Unsupported Itemtype "<<itemType;
          std::logic_error("halt");
        }
      }
    }
    default:
    {
      std::cerr<< "ExplicitLagrangian Basis Service: Unsupported order "<<order;
      std::logic_error("halt");
    }
  } 
}


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