function Lprime=GetLprimeStokes(triaelem,xyz_list,gauss_coord_tria,gauss_coord_penta)
%GETLPRIMESTOKES - compute Lprime matrix for stokes diagnostic
%
%   Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
%   For grid i, Lpi can be expressed in the basic coordinate system by:
%   Lpi_basic=[ h    0    0   0]
%             [ 0    h    0   0]
%             [ h    0    0   0]
%             [ 0    h    0   0]
%             [ 0    0    h   0]
%             [ 0    0    h   0]
%             [ 0    0  dh/dz 0]
%             [ 0    0  dh/dz 0]
%             [ 0    0  dh/dz 0]
%             [dh/dz 0  dh/dx 0]
%             [ 0 dh/dz dh/dy 0]
%             [ 0    0    0   h]
%             [ 0    0    0   h]
%             [ 0    0    0   h]
%   where h is the interpolation function for grid i.
%
%   Usage:
%      Lprime=GetLprimeStokes(triaelem,xyz_list,gauss_coord_tria,gauss_coord_penta)
%
%   See also GETLSTOKES, GETL, CREATEKMATRIX

numdof=4;
num_grids=3;

%Get l1l2l3 in basic coordinate system: 
l1l2l3=GetNodalFunctions(triaelem,gauss_coord_tria);
dh1dh6=GetNodalFunctionsDerivativesBasic(pentaelem,xyz_list,gauss_coord_penta);

%Build Lprime:
Lprime=zeros(12,numdof*num_grids);

for i=1:num_grids,
	Lprime(1,numdof*(i-1)+1)=l1l2l3(i); 
	Lprime(2,numdof*(i-1)+2)=l1l2l3(i);
	Lprime(3,numdof*(i-1)+1)=l1l2l3(i); 
	Lprime(4,numdof*(i-1)+2)=l1l2l3(i); 
	Lprime(5,numdof*(i-1)+3)=l1l2l3(i); 
	Lprime(6,numdof*(i-1)+3)=l1l2l3(i); 
	Lprime(7,numdof*(i-1)+3)=dh1dh6(3,i); 
	Lprime(8,numdof*(i-1)+3)=dh1dh6(3,i); 
	Lprime(9,numdof*(i-1)+3)=dh1dh6(3,i); 
	Lprime(10,numdof*(i-1)+3)=dh1dh6(1,i); 
	Lprime(10,numdof*(i-1)+1)=dh1dh6(3,i); 
	Lprime(11,numdof*(i-1)+3)=dh1dh6(2,i); 
	Lprime(11,numdof*(i-1)+2)=dh1dh6(3,i); 
	Lprime(12,numdof*(i-1)+4)=l1l2l3(i); 
	Lprime(13,numdof*(i-1)+4)=l1l2l3(i); 
	Lprime(14,numdof*(i-1)+4)=l1l2l3(i); 
	
%other terms are 0
end
