function L=GetL(triaelem,gauss_l1l2l3,numdof)
%GETL - compute L matrix for a 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:
%
%   numdof=1: 
%   Li_basic=h;
%
%   numdof=2:
%   Li_basic=[ h    0  ]
%            [   0   h ]
%   numdof=4:
%   Li_basic=[ h    0     0    0 ]
%            [ 0    h     0    0 ]
%            [ 0    0     0    0 ]
%            [ 0    0     0    0 ]
%   where h is the interpolation function for grid i.
%
%   Usage:
%      L=GetL(triaelem,gauss_l1l2l3,numdof)
%
%   See also CREATEKMATRIX, GETLSTOKES, GETLPRIMESTOKES

num_grids=3;

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

%Build L:
if numdof==1,
	L=zeros(1,num_grids);
	for i=1:num_grids,
		L(i)=l1l2l3(i);
	end
elseif numdof==2,
	L=zeros(numdof,numdof*num_grids);

	for i=1:num_grids,
		L(1,numdof*(i-1)+1)=l1l2l3(i); 
		L(1,numdof*(i-1)+2)=0;

		L(2,numdof*(i-1)+1)=0;
		L(2,numdof*(i-1)+2)=l1l2l3(i);
	end
elseif numdof==4,
	L=zeros(4,numdof*num_grids);

	for i=1:num_grids,
		L(1,numdof*(i-1)+1)=l1l2l3(i); 
		L(1,numdof*(i-1)+2)=0;
		L(1,numdof*(i-1)+3)=0;
		L(1,numdof*(i-1)+4)=0;

		L(2,numdof*(i-1)+1)=0;
		L(2,numdof*(i-1)+2)=l1l2l3(i);
		L(2,numdof*(i-1)+3)=0;
		L(2,numdof*(i-1)+4)=0;

		L(3,numdof*(i-1)+1)=0;
		L(3,numdof*(i-1)+2)=0;
		L(3,numdof*(i-1)+3)=0;
		L(3,numdof*(i-1)+4)=0;

		L(4,numdof*(i-1)+1)=0;
		L(4,numdof*(i-1)+2)=0;
		L(4,numdof*(i-1)+3)=0;
		L(4,numdof*(i-1)+4)=0;
	end
end
