function Bs = BatGaussPoints(hh,nGrid)
h=1.0/nGrid;
X = h/2:h:1-h/2;
X = X.*hh;
[x,y,z]=meshgrid(X);

Bs=LocalBs(x,y,z,0,0,0,hh);
% intB = sum(Bs,[3 4 5]).*hh^3/nGrid^3;
% intB = zeros(6,24);
% for ix=1:nGrid
%     for iy=1:nGrid
%         for iz=1:nGrid
%             intB=intB+Bs(:,:,ix,iy,iz);
%         end
%     end
% end
% intB=intB.*hh^3/nGrid^3;
end