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

#include "FiniteElementMapping.hpp"
#include "MappingBasisFunction.hpp"
#include "Quadrature.hpp"
#include "FiniteElementMappingBasis.hpp"
#include <math.h>
#include "CollapsedMapping.hpp"
#include "ReferenceFiniteElement.hpp"
#include "Tuple.hpp"
#include "EnumElemType.hpp"
#include "Eigen/Dense"
#include "Mesh.hpp"

template <class BasisFunctionT>
class ReferenceSpaceMappingFunction 
{
public:
  typedef FiniteElementMappingBasis<BasisFunctionT> FEMappingFunctionType;
  typedef Quadrature::PointVector PointVector;
  typedef SpaMesh::Connectivity::TopoPointVector TopoPointVector;
  typedef SpaMesh::Connectivity::TopoIndexVector TopoIndexVector;

  ReferenceSpaceMappingFunction() 
    {
      ;
    }
  
  ~ReferenceSpaceMappingFunction() {}
  
public:
  template<unsigned int Dim, unsigned int DIm,
           template <unsigned int, unsigned int> class MappingBasisFunctionType>
       static void initMappingFunction(MappingBasisFunctionType<Dim,DIm>& mappingFunction, 
	                               const Quadrature& quadrature, 
				       const ElemType& elemType, 
				       const TopoPointVector& coord, 
				       const int morder)
  {
    assert(quadrature.nbPoint()>0);
    FEMappingFunctionType::init(mappingFunction, quadrature, elemType, morder);

    const unsigned int n_qPoints = quadrature.nbPoint();
    mappingFunction.resizeMapping(n_qPoints);
    computeMap(mappingFunction, coord, n_qPoints);
    mappingFunction.getJacobianXweight().noalias() = mappingFunction.getJacobian().cwiseProduct(quadrature.getWeight());
  } 
  
  template<unsigned int Dim, unsigned int DIm, 
           template <unsigned int, unsigned int> class MappingBasisFunctionType>
     static void computePhysicalSpacePointsX(MappingBasisFunctionType<Dim,DIm>& mappingFunction, 
	                                     const Quadrature::PointVector& refSpacePoints, 
	                                     const ElemType& elemType, 
					     const TopoPointVector& coord, 
					     const int morder)
  {
    mappingFunction.doComputeShapeDeriv(false);
    const typename BasisFunction<DIm>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    FEMappingFunctionType::reInit(mappingFunction, refSpacePoints, elemType, morder);
    assert(coord.rows() == shapeMatrix.rows());
    mappingFunction.resizeQpointsX(refSpacePoints.rows());
    Quadrature::PointVector& pointsX = mappingFunction.getQpointX();
    for (unsigned int qp = 0; qp < refSpacePoints.rows() ; qp++)
      for (unsigned int a = 0; a < DIm; a++)
        pointsX(qp,a) = shapeMatrix.col(qp).dot(coord.col(a));    
    mappingFunction.doComputeShapeDeriv(true);
  } 
  
  template<unsigned int Dim, unsigned int DIm, 
           template <unsigned int, unsigned int> class MappingBasisFunctionType>
     static void computePhysicalSpacePointsWithMapX(MappingBasisFunctionType<Dim,DIm>& mappingFunction, 
	                                            const Quadrature::PointVector& refSpacePoints, 
	                                            const ElemType& elemType, 
					            const TopoPointVector& coord, 
					            const int morder)
  {
    FEMappingFunctionType::reInit(mappingFunction, refSpacePoints, elemType, morder);
    
    const unsigned int n_qPoints = refSpacePoints.rows();
    mappingFunction.resizeMapping(n_qPoints);
    computeMap(mappingFunction, coord, n_qPoints);
  } 

  template<unsigned int DIm,
           template <unsigned int, unsigned int> class MappingBasisFunctionType>
     static void reInitMappingFunctionElemSide(MappingBasisFunctionType<DIm,DIm>& mappingFunction, 
	                                       const typename Quadrature::PointVector& qVec, 
					       const ElemType& elemType, 
					       const TopoPointVector& coord, 
					       const int morder)
  {
    assert(qVec.size()>0);
    mappingFunction.doComputeShape(false);  
    FEMappingFunctionType::reInit(mappingFunction, qVec, elemType, morder);
 
    mappingFunction.resizeMapping(qVec.rows());
    computeMapLite(mappingFunction, coord, qVec.rows()); 
  } 
  
  static void referenceSpaceSideMap(const ElemType& elemType, const ElemType& sideType, const TopoIndexVector& side_indices,
	                            const MappingBasisFunction<2,3>& mappingFunction, PointVector& elemQuadrature)
  {
    assert(mappingFunction.nbShapeMode() == side_indices.size());  
    const typename BasisFunction<2>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector vertices;
    ReferenceFiniteElement::getVerticesE(elemType, vertices);

    for (unsigned int qp = 0; qp < elemQuadrature.rows(); qp++)
    {
      elemQuadrature.row(qp).setZero();
      for (unsigned int i = 0; i < side_indices.size(); i++)
        for (unsigned int a = 0; a < 3; a++)
          elemQuadrature(qp,a) += shapeMatrix(i,qp) * vertices(side_indices[i],a);    
    }
  }

  static void referenceSpaceSideMap(const ElemType& elemType, const ElemType& sideType, const TopoIndexVector& side_indices,
	                            const MappingBasisFunction<1,2>& mappingFunction, PointVector& elemQuadrature)
  {
    assert(mappingFunction.nbShapeMode() == side_indices.size()); 
    const typename BasisFunction<1>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector vertices;
    ReferenceFiniteElement::getVerticesE(elemType, vertices);

    for (unsigned int qp = 0; qp < elemQuadrature.rows(); qp++)
    {
      elemQuadrature.row(qp).setZero();
      for (unsigned int i = 0; i < side_indices.size(); i++)
        for (unsigned int a = 0; a < 2; a++)
          elemQuadrature(qp,a) += shapeMatrix(i,qp) * vertices(side_indices[i],a);    
    }
  }
  
  static void referenceSpaceSideMap(const ElemType& elemType, const ElemType& sideType, const TopoIndexVector& side_indices,
	                            MappingBasisFunction<0,1>& mappingFunction, PointVector& elemQuadrature)
  {
    assert(mappingFunction.nbShapeMode() == side_indices.size()); 
    const typename BasisFunction<0>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector vertices;
    ReferenceFiniteElement::getVerticesE(elemType, vertices);
    assert(elemQuadrature.rows() == 1);
    assert(side_indices.size() == 1);

    elemQuadrature(0,0) = vertices(side_indices[0],0);    
   
    // compute Normals for the 1D case
    typename MappingBasisFunction<0,1>::NormalVector& normalVec = mappingFunction.getNormalVec();
    if (side_indices[0] == 0)
    {
      normalVec(0)(0) = -1.;
    }
    else
    {
      assert (side_indices[0] == 1);
      normalVec(0)(0) = 1.;
    }
  }

private:
  static void computeMap(MappingBasisFunction<0,0>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compude QpointX
    PointVector& qPointX = mappingFunction.getQpointX();
    qPointX.col(0) = coord.col(0);    
    // compute jacobian
    typename MappingBasisFunction<0,0>::JacVector& jac = mappingFunction.getJacobian();
    jac(0) = 1.;
  } 
  static void computeMap(MappingBasisFunction<0,1>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compude QpointX
    PointVector& qPointX = mappingFunction.getQpointX();
    qPointX.col(0) = coord.col(0);    
    // compute jacobian
    typename MappingBasisFunction<0,1>::JacVector& jac = mappingFunction.getJacobian();
    jac(0) = 1.;
  } 
  static void computeMap(MappingBasisFunction<1,1>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compude QpointX
    const typename BasisFunction<1>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector& qPointX = mappingFunction.getQpointX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      qPointX(qp,0) = shapeMatrix.col(qp).dot(coord.col(0));    
    //compute dXdE 
    const typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<1,1>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      dXdE(qp)(0,0) = shapeDerivMatrix(0).col(qp).dot(coord.col(0));
    typename MappingBasisFunction<1,1>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<1,1>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      // compute jacobian
      jac(qp) = dXdE(qp)(0,0);
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< jac(qp) << " in 1D element "<<std::endl;
        std::logic_error("halt");
      } 
      // compute dEdX
      typename MappingBasisFunction<1,1>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
      dEdX(qp)(0,0) = dXdE(qp)(0,0)/jac(qp)/jac(qp);
    }
  } 
  static void computeMap(MappingBasisFunction<1,2>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compude QpointX
    const typename BasisFunction<1>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector& qPointX = mappingFunction.getQpointX();
    for (unsigned int a = 0; a < 2; a++)
      for (unsigned int qp = 0; qp < n_qPoints; qp++)
        qPointX(qp,a) = shapeMatrix.col(qp).dot(coord.col(a));    
    //compute dXdE 
    const typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<1,2>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int a = 0; a < 2; a++)
        dXdE(qp)(a,0) = shapeDerivMatrix(0).col(qp).dot(coord.col(a));
    typename MappingBasisFunction<1,2>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<1,2>::NormalVector& normalVec = mappingFunction.getNormalVec();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      // compute jacobian
      for (unsigned int a = 0; a < 2; a++)
      {
        jac(qp) += dXdE(qp)(a,0)*dXdE(qp)(a,0);
      }
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< jac(qp) << " in 1D side element "<<std::endl;
        std::logic_error("halt");
      } 
      jac(qp) = sqrt(jac(qp));
      // compute Normals
      normalVec(qp)(0) =  dXdE(qp)(1,0);   
      normalVec(qp)(1) = -dXdE(qp)(0,0);
      normalVec(qp) /= normalVec(qp).norm();
    }
  } 
  static void computeMap(MappingBasisFunction<2,2>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());//,("Error in ReferenceSpaceMappingFunctionElem: nodes and shape functions size must agree"));
    //compude QpointX
    const typename BasisFunction<2>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector& qPointX = mappingFunction.getQpointX();
    for (unsigned int a = 0; a < 2; a++)
      for (unsigned int qp = 0; qp < n_qPoints; qp++)
        qPointX(qp,a) = shapeMatrix.col(qp).dot(coord.col(a));    
    //compute dXdE 
    const typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<2,2>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int b = 0; b < 2; b++)
        for (unsigned int a = 0; a < 2; a++)
          dXdE(qp)(a,b) = shapeDerivMatrix(b).col(qp).dot(coord.col(a));
    typename MappingBasisFunction<2,2>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<2,2>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      // compute jacobian
      jac(qp) = dXdE(qp)(0,0) * dXdE(qp)(1,1) - dXdE(qp)(0,1) * dXdE(qp)(1,0);
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< jac(qp) << " in 2D element "<<std::endl;
        std::logic_error("halt");
      } 
      // compute dEdX
      const double inv_jac = 1./jac(qp);
      dEdX(qp)(0,0) = dXdE(qp)(1,1) * inv_jac; 
      dEdX(qp)(0,1) = - dXdE(qp)(0,1) * inv_jac;
      dEdX(qp)(1,0) = - dXdE(qp)(1,0) * inv_jac;
      dEdX(qp)(1,1) = dXdE(qp)(0,0) * inv_jac; 
    }
  }
  static void computeMap(MappingBasisFunction<2,3>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compute QpointX
    const typename BasisFunction<2>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector& qPointX = mappingFunction.getQpointX();
    for (unsigned int a = 0; a < 3; a++)
      for (unsigned int qp = 0; qp < n_qPoints; qp++)
        qPointX(qp,a) = shapeMatrix.col(qp).dot(coord.col(a));    
    //compute dXdE 
    const typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<2,3>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int b = 0; b < 2; b++)
        for (unsigned int a = 0; a < 3; a++)
          dXdE(qp)(a,b) = shapeDerivMatrix(b).col(qp).dot(coord.col(a));
    typename MappingBasisFunction<2,3>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<2,3>::NormalVector& normalVec = mappingFunction.getNormalVec();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      const double 
               dXdE00 = dXdE(qp)(0,0), dXdE10 = dXdE(qp)(1,0), dXdE20 = dXdE(qp)(2,0),
               dXdE01 = dXdE(qp)(0,1), dXdE11 = dXdE(qp)(1,1), dXdE21 = dXdE(qp)(2,1);
      
      // compute jacobian
      const double g11 = (dXdE00 * dXdE00 +
                        dXdE10 * dXdE10 +
                        dXdE20 * dXdE20);
          
      const double g12 = (dXdE00 * dXdE01 +
                        dXdE10 * dXdE11 +
                        dXdE20 * dXdE21);
          
      const double g21 = g12;
          
      const double g22 = (dXdE01 * dXdE01 +
                        dXdE11 * dXdE11 +
                        dXdE21 * dXdE21);
  
      const double det = (g11*g22 - g12*g21);
      if (det <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< det << " in 2D side element "<<std::endl;
        std::logic_error("halt");
      } 
      jac(qp) = sqrt(det);
      
      // compute Normals
      normalVec(qp)(0) =   dXdE10*dXdE21 - dXdE20*dXdE11;
      normalVec(qp)(1) = - dXdE00*dXdE21 + dXdE20*dXdE01;
      normalVec(qp)(2) =   dXdE00*dXdE11 - dXdE10*dXdE01; 
      normalVec(qp) /= normalVec(qp).norm();
    }
  }
  static void computeMap(MappingBasisFunction<3,3>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compude QpointX
    const typename BasisFunction<3>::ShapeMatrix& shapeMatrix = mappingFunction.getShapeFunction();
    PointVector& qPointX = mappingFunction.getQpointX();
    for (unsigned int a = 0; a < 3; a++)
      for (unsigned int qp = 0; qp < n_qPoints; qp++)
        qPointX(qp,a) = shapeMatrix.col(qp).dot(coord.col(a));    
    //compute dXdE 
    const typename BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<3,3>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int b = 0; b < 3; b++)
        for (unsigned int a = 0; a < 3; a++)
          dXdE(qp)(a,b) = shapeDerivMatrix(b).col(qp).dot(coord.col(a));

    typename MappingBasisFunction<3,3>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<3,3>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      const double 
               dXdE00 = dXdE(qp)(0,0), dXdE10 = dXdE(qp)(1,0), dXdE20 = dXdE(qp)(2,0),
               dXdE01 = dXdE(qp)(0,1), dXdE11 = dXdE(qp)(1,1), dXdE21 = dXdE(qp)(2,1),
               dXdE02 = dXdE(qp)(0,2), dXdE12 = dXdE(qp)(1,2), dXdE22 = dXdE(qp)(2,2);
      
      // compute jacobian
      jac(qp) = 0;
      jac(qp) = (dXdE00*(dXdE11*dXdE22 - dXdE21*dXdE12)  +
                 dXdE10*(dXdE21*dXdE02 - dXdE01*dXdE22)  +
                 dXdE20*(dXdE01*dXdE12 - dXdE11*dXdE02));
  
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian "<< jac(qp) << " in 3D element "<<std::endl;
        std::logic_error("halt");
	exit(0);
      } 
      const double inv_jac  = 1./jac(qp);       
          
      // compute dEdX
      dEdX(qp)(0,0) = (dXdE11*dXdE22 - dXdE21*dXdE12)*inv_jac;
      dEdX(qp)(0,1) = (dXdE21*dXdE02 - dXdE01*dXdE22)*inv_jac;
      dEdX(qp)(0,2) = (dXdE01*dXdE12 - dXdE11*dXdE02)*inv_jac;
                   
      dEdX(qp)(1,0) = (dXdE20*dXdE12 - dXdE10*dXdE22)*inv_jac;
      dEdX(qp)(1,1) = (dXdE00*dXdE22 - dXdE20*dXdE02)*inv_jac;
      dEdX(qp)(1,2) = (dXdE10*dXdE02 - dXdE00*dXdE12)*inv_jac;
                   
      dEdX(qp)(2,0) = (dXdE10*dXdE21 - dXdE20*dXdE11)*inv_jac;
      dEdX(qp)(2,1) = (dXdE20*dXdE01 - dXdE00*dXdE21)*inv_jac;
      dEdX(qp)(2,2) = (dXdE00*dXdE11 - dXdE10*dXdE01)*inv_jac;
    }
  }

  static void computeMapLite(MappingBasisFunction<1,1>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compute dXdE 
    const typename BasisFunction<1>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<1,1>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      dXdE(qp)(0,0) = shapeDerivMatrix(0).col(qp).dot(coord.col(0));
    typename MappingBasisFunction<1,1>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<1,1>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      // compute jacobian
      jac(qp) = dXdE(qp)(0,0);
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< jac(qp) << " in 1D reference element side"<<std::endl;
        std::logic_error("halt");
      } 
      // compute dEdX
      typename MappingBasisFunction<1,1>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
      dEdX(qp)(0,0) = dXdE(qp)(0,0)/jac(qp)/jac(qp);
    }
  } 

  static void computeMapLite(MappingBasisFunction<2,2>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compute dXdE 
    const typename BasisFunction<2>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<2,2>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int b = 0; b < 2; b++)
        for (unsigned int a = 0; a < 2; a++)
          dXdE(qp)(a,b) = shapeDerivMatrix(b).col(qp).dot(coord.col(a));
    typename MappingBasisFunction<2,2>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<2,2>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      jac(qp) = dXdE(qp)(0,0) * dXdE(qp)(1,1) - dXdE(qp)(0,1) * dXdE(qp)(1,0);
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian"<< jac(qp) << " in 2D reference element side"<<std::endl;
        std::logic_error("halt");
      } 
      // compute dEdX
      const double inv_jac = 1./jac(qp);
      dEdX(qp)(0,0) = dXdE(qp)(1,1) * inv_jac; 
      dEdX(qp)(0,1) = - dXdE(qp)(0,1) * inv_jac;
      dEdX(qp)(1,0) = - dXdE(qp)(1,0) * inv_jac;
      dEdX(qp)(1,1) = dXdE(qp)(0,0) * inv_jac; 
    }
  }

  static void computeMapLite(MappingBasisFunction<3,3>& mappingFunction, const TopoPointVector& coord, const unsigned int n_qPoints)
  {
    assert(mappingFunction.nbShapeMode() == coord.rows());
    //compute dXdE 
    const typename BasisFunction<3>::ShapeDerivMatrix& shapeDerivMatrix = mappingFunction.getShapeFunctionDeriv(); 
    typename MappingBasisFunction<3,3>::MappingMatrix& dXdE = mappingFunction.getdXdE();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
      for (unsigned int b = 0; b < 3; b++)
        for (unsigned int a = 0; a < 3; a++)
          dXdE(qp)(a,b) = shapeDerivMatrix(b).col(qp).dot(coord.col(a));
    typename MappingBasisFunction<3,3>::JacVector& jac = mappingFunction.getJacobian();
    typename MappingBasisFunction<3,3>::MappingMatrixInv& dEdX = mappingFunction.getdEdX();
    for (unsigned int qp = 0; qp < n_qPoints; qp++)
    {
      const double 
               dXdE00 = dXdE(qp)(0,0), dXdE10 = dXdE(qp)(1,0), dXdE20 = dXdE(qp)(2,0),
               dXdE01 = dXdE(qp)(0,1), dXdE11 = dXdE(qp)(1,1), dXdE21 = dXdE(qp)(2,1),
               dXdE02 = dXdE(qp)(0,2), dXdE12 = dXdE(qp)(1,2), dXdE22 = dXdE(qp)(2,2);
      
      // compute jacobian
      jac(qp) = 0;
      jac(qp) = (dXdE00*(dXdE11*dXdE22 - dXdE21*dXdE12)  +
                 dXdE10*(dXdE21*dXdE02 - dXdE01*dXdE22)  +
                 dXdE20*(dXdE01*dXdE12 - dXdE11*dXdE02));
  
      if (jac(qp) <= 0.)
      {
        std::cerr <<" Error in ReferenceSpaceMappingFunction: negative Jacobian "<< jac(qp) << " in 3D reference element side"<<std::endl;
        std::logic_error("halt");
      } 
      const double inv_jac  = 1./jac(qp);       
          
      // compute dEdX
      dEdX(qp)(0,0) = (dXdE11*dXdE22 - dXdE21*dXdE12)*inv_jac;
      dEdX(qp)(0,1) = (dXdE21*dXdE02 - dXdE01*dXdE22)*inv_jac;
      dEdX(qp)(0,2) = (dXdE01*dXdE12 - dXdE11*dXdE02)*inv_jac;
                   
      dEdX(qp)(1,0) = (dXdE20*dXdE12 - dXdE10*dXdE22)*inv_jac;
      dEdX(qp)(1,1) = (dXdE00*dXdE22 - dXdE20*dXdE02)*inv_jac;
      dEdX(qp)(1,2) = (dXdE10*dXdE02 - dXdE00*dXdE12)*inv_jac;
                   
      dEdX(qp)(2,0) = (dXdE10*dXdE21 - dXdE20*dXdE11)*inv_jac;
      dEdX(qp)(2,1) = (dXdE20*dXdE01 - dXdE00*dXdE21)*inv_jac;
      dEdX(qp)(2,2) = (dXdE00*dXdE11 - dXdE10*dXdE01)*inv_jac;
    }
  }
};

#endif 
