function gqp = GaussQuadPoints()
%gqp = GaussQuadPoints()    Calculate the Gaussian quadrature points
%    for a triangular prism.  The result is a 3*6 matrix.

    onesixth = 1/6;
    twothirds = 2/3;
    sqrtthird = sqrt(1/3);

    gqp = [ ...
            [  onesixth;  onesixth; -sqrtthird ], ...
            [ twothirds;  onesixth; -sqrtthird ], ...
            [  onesixth; twothirds; -sqrtthird ], ...
            [  onesixth;  onesixth;  sqrtthird ], ...
            [ twothirds;  onesixth;  sqrtthird ], ...
            [  onesixth; twothirds;  sqrtthird ], ...
          ];
end
