function B=GetBStokes(pentaelem,xyz_list,gauss_coord)
%GETBSTOKES - computes the B matrix for Stokes' diagnostic
%
%   Compute B  matrix. B=[B1 B2 B3 B4 B5 B6 Bb] 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       0  ]
%            [   0           dh/dy           0       0  ]
%            [   0             0           dh/dy     0  ]
%            [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
%            [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
%            [   0          1/2*dh/dz    1/2*dh/dy   0  ]
%            [   0             0             0       h  ]
%            [ dh/dx         dh/dy         dh/dz     0  ]
%   where h is the interpolation function for grid i.
%   Same thing for Bb except the last column that does not exist.
% 
%   Usage:
%      B=GetBStokes(pentaelem,xyz_list,gauss_coord)
%
%   See also GETBPRIMESTOKES, CREATEKMATRIX

	num_grids=6;	 % the bubble in not included
	NDOF4=4;

	l1l6=GetNodalFunctions(pentaelem, gauss_coord);
	dh1dh7_basic=GetNodalFunctionsDerivativesBasicStokes(pentaelem,xyz_list,gauss_coord);

	%Build B:
	B=zeros(8,NDOF4*num_grids+3);
	for i=1:num_grids+1,
		B(1,NDOF4*(i-1)+1)=dh1dh7_basic(1,i);
		B(1,NDOF4*(i-1)+2)=0;
		B(1,NDOF4*(i-1)+3)=0;

		B(2,NDOF4*(i-1)+1)=0;
		B(2,NDOF4*(i-1)+2)=dh1dh7_basic(2,i);
		B(2,NDOF4*(i-1)+3)=0;

		B(3,NDOF4*(i-1)+1)=0;
		B(3,NDOF4*(i-1)+2)=0;
		B(3,NDOF4*(i-1)+3)=dh1dh7_basic(3,i);

		B(4,NDOF4*(i-1)+1)=dh1dh7_basic(2,i)/2;
		B(4,NDOF4*(i-1)+2)=dh1dh7_basic(1,i)/2;	
		B(4,NDOF4*(i-1)+3)=0;

		B(5,NDOF4*(i-1)+1)=dh1dh7_basic(3,i)/2;
		B(5,NDOF4*(i-1)+2)=0;
		B(5,NDOF4*(i-1)+3)=dh1dh7_basic(1,i)/2;

		B(6,NDOF4*(i-1)+1)=0;
		B(6,NDOF4*(i-1)+2)=dh1dh7_basic(3,i)/2;
		B(6,NDOF4*(i-1)+3)=dh1dh7_basic(2,i)/2;

		B(7,NDOF4*(i-1)+1)=0;
		B(7,NDOF4*(i-1)+2)=0;
		B(7,NDOF4*(i-1)+3)=0;

		B(8,NDOF4*(i-1)+1)=dh1dh7_basic(1,i);
		B(8,NDOF4*(i-1)+2)=dh1dh7_basic(2,i);
		B(8,NDOF4*(i-1)+3)=dh1dh7_basic(3,i);
	end
	
	%Add the columns corresponding to pressure (not for the bubble)
	for i=1:num_grids,
		B(7,NDOF4*(i-1)+4)=l1l6(i);
	end
