function L=GetLStokes(triaelem,gauss_l1l2l3)
%GETLSTOKES - compute L matrix for stokes diagnostic
%
%   Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
%   For grid i, Li can be expressed in the basic coordinate system by:
%   Li_basic=[ h    0    0   0]
%	          [ 0    h    0   0]
%            [ 0    0    h   0]
%            [ 0    0    h   0]
%            [ 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    h   0]
%            [ h    0    0   0]
%            [ 0    h    0   0]
%            [ 0    0    h   0]
%   where h is the interpolation function for grid i.
% 
%   Usage:
%
%   See also GETL, GETLPRIMESTOKES, CREATEKMATRIX

numdof=4;
num_grids=3;

%Get l1l2l3 in basic coordinate system: 
l1l2l3=GetNodalFunctions(triaelem,gauss_l1l2l3);

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

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