function B=GetB_horiz(pentaelem,xyz_list,gauss_coord)
%GETB_HORIZ - compute the B matrix for the horizontal velocity of Pattyn's model
%
%   Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size numgrids*NDOF
%   For grid i, Bi can be expressed in the basic coordinate system by:
%   Bi_basic=[ dh/dx          0     ]
%            [   0          dh/dy   ]
%            [ 1/2*dh/dy  1/2*dh/dx ]
%            [ 1/2*dh/dz      0     ]
%            [   0        1/2*dh/dz ]
%   where h is the interpolation function for grid i.
%
%   Usage:
%      B=GetB_horiz(pentaelem,xyz_list,gauss_coord)
% 
%   See also GETBPRIME_HORIZ, CREATEKMATRIX

	num_grids=6;	
	NDOF2=2;

	dh1dh6_basic=GetNodalFunctionsDerivativesBasic(pentaelem,xyz_list,gauss_coord);
	%Build B:
	B=zeros(5,NDOF2*num_grids);
	for i=1:num_grids,
		B(1,NDOF2*(i-1)+1)=dh1dh6_basic(1,i);
		B(1,NDOF2*(i-1)+2)=0;
		B(2,NDOF2*(i-1)+1)=0;
		B(2,NDOF2*(i-1)+2)=dh1dh6_basic(2,i);
		B(3,NDOF2*(i-1)+1)=dh1dh6_basic(2,i)/2;
		B(3,NDOF2*(i-1)+2)=dh1dh6_basic(1,i)/2;	
		B(4,NDOF2*(i-1)+1)=dh1dh6_basic(3,i)/2;
		B(4,NDOF2*(i-1)+2)=0;
		B(5,NDOF2*(i-1)+1)=0;
		B(5,NDOF2*(i-1)+2)=dh1dh6_basic(3,i)/2;
	end
