function J=GetJacobian(pentaelem,xyz_list,gauss_coord)
%GETJACOBIAN - computes the jacobian for a pentaelem
%
%   Usage:
%      J=GetJacobian(pentaelem,xyz_list,gauss_coord)
%
%   See also GETJACOBIANDETERMINANT, GETJACOBIANINVERT

%The Jacobian is constant over the element, discard the gaussian points. 

J=zeros(3,3);
sqrt3=sqrt(3.0);

%grids coordinates	
coord1=xyz_list(1,:)';
coord2=xyz_list(2,:)';
coord3=xyz_list(3,:)';
coord4=xyz_list(4,:)';
coord5=xyz_list(5,:)';
coord6=xyz_list(6,:)';

%Gauss point coordinate in reference element
r=gauss_coord(1,2) - gauss_coord(1,1);
s= - 3.0/sqrt3 *(gauss_coord(1,1)+gauss_coord(1,2) - 2.0/3.0);
zeta=gauss_coord(1,4);

%Computation of the jacobian on this gauss point
J(:,1)=1.0/4.0* (-coord1+coord2-coord4+coord5) + 1.0/4.0 *(coord1-coord2-coord4+coord5)*zeta;
J(:,2)=sqrt3/12.0 *(-coord1-coord2+2*coord3-coord4-coord5+2*coord6) +  sqrt3/12.0* (coord1+coord2-2*coord3-coord4-coord5+2*coord6)*zeta;
J(:,3)=1.0/4.0* (coord1-coord2-coord4+coord5)*r +  sqrt3/12.0 *(coord1+coord2-2*coord3-coord4-coord5+2*coord6)*s + 1/6.0 *(-coord1-coord2-coord3+coord4+coord5+coord6) ;

%Faster than :
%J(1,1)=1.0/4.0* (-x1+x2-x4+x5) + 1.0/4.0 *(x1-x2-x4+x5)*zeta;
%J(2,1)=1.0/4.0 *(-y1+y2-y4+y5) + 1.0/4.0 *(y1-y2-y4+y5)*zeta;
%J(3,1)=1.0/4.0 *(-z1+z2-z4+z5) + 1.0/4.0 *(z1-z2-z4+z5)*zeta;

%J(1,2)=sqrt3/12.0 *(-x1-x2+2*x3-x4-x5+2*x6) +  sqrt3/12.0* (x1+x2-2*x3-x4-x5+2*x6)*zeta;
%J(2,2)=sqrt3/12.0 *(-y1-y2+2*y3-y4-y5+2*y6) +  sqrt3/12.0* (y1+y2-2*y3-y4-y5+2*y6)*zeta;
%J(3,2)=sqrt3/12.0 *(-z1-z2+2*z3-z4-z5+2*z6) +  sqrt3/12.0* (z1+z2-2*z3-z4-z5+2*z6)*zeta;

%J(1,3)=1.0/4.0* (x1-x2-x4+x5)*r +  sqrt3/12.0 *(x1+x2-2*x3-x4-x5+2*x6)*s + 1/6.0 *(-x1-x2-x3+x4+x5+x6) ;
%J(2,3)=1.0/4.0* (y1-y2-y4+y5)*r +  sqrt3/12.0 *(y1+y2-2*y3-y4-y5+2*y6)*s + 1/6.0 *(-y1-y2-y3+y4+y5+y6) ;
%J(3,3)=1.0/4.0* (z1-z2-z4+z5)*r +  sqrt3/12.0 *(z1+z2-2*z3-z4-z5+2*z6)*s + 1/6.0 *(-z1-z2-z3+z4+z5+z6) ;


