function [W, HP] = GaussianPoints3D(NPoints)
% This program is for getting the Weight and Points of Gaussian Integration L1 and L2
if(~isnumeric(NPoints))
    error('Input order of Gaussian NPOINTS must be numeric and integer');
elseif(~ismember(NPoints, [1,2,3,4,5]))
    error('Input order of Gaussian NPOINTS must be 1, 2, 3, 4 or 5 ')
end


switch NPoints
    case 1
        WW = 2;
        HH = 0;
    case 2
        WW = [1,1];
        HH = [-1/sqrt(3), 1/sqrt(3)];
    case 3
        p = sqrt(3/5);
        HH = [-p, 0, p];
        WW = [5/9, 8/9, 5/9];
    case 4
        p = 3/7;
        q = 2/7*sqrt(6/5);
        r = sqrt(p+q);
        s = sqrt(p-q);
        HH = [-r, -s, s, r];
        p = (18+sqrt(30))/36;
        q = (18-sqrt(30))/36;
        WW = [q, p, p, q];
    case 5
        p = 5;
        q = 2*sqrt(10/7);
        r = sqrt(p+q)/3;
        s = sqrt(p-q)/3;
        HH = [-r, -s, 0, s, r];
        p = (322+13*sqrt(70))/900;
        q = (322-13*sqrt(70))/900;
        r = 128/225;
        WW = [q, p, r, p, q];
    otherwise
        error('No Gaussian points of order %d in 3d parametric element!', NPoints)
end
[Hx, Hy, Hz] = meshgrid(HH);
[Wx, Wy, Wz] = meshgrid(WW);
HP = [Hx(:), Hy(:), Hz(:)];
W =  Wx(:) .* Wy(:) .* Wz(:);
end