function [ngp,w,gp] = gauss(type,order)
% Get gauss points coordinates in the parametric space of element and corresponding weights for a given quadrature type and order.
% Input arguments:
%  type:  gauss quadrature type (either line, triangle or quadrilateral)
%  order: quadrature order
% Output arguments:
%  ngp:   number of gauss points
%  w:     array of gauss point weights
%  gp:    array of gauss point parametric coordinates

 mlock;                   % prevent clearing this function from memory

 include_gblrefs;

 persistent qpts;
 persistent qwgt;

 if isempty(qpts)
  qpts = [.000000000000000 -.577350269189626 -.774596669241483 -.861136311594053
          .000000000000000  .577350269189626  .000000000000000 -.339981043584856
          .000000000000000  .000000000000000  .774596669241483  .339981043584856
          .000000000000000  .000000000000000  .000000000000000  .861136311594053];
 end

 if isempty(qwgt)
  qwgt = [2.00000000000000 1.000000000000000 0.555555555555555 0.347854845137454
          0.00000000000000 1.000000000000000 0.888888888888888 0.652145154862546
          0.00000000000000 0.000000000000000 0.555555555555555 0.652145154862546
          0.00000000000000 0.000000000000000 0.000000000000000 0.347854845137454];
 end

 % Initialize arrays of gauss point parametric coordenates and weights
 if(type == LINE_QUADRATURE)
  ngp = order;
 elseif(type == TRIA_QUADRATURE)
  if( order == 1 )
   ngp = 1;
  elseif( order == 2 || order == 3 )
   ngp = 3;
  elseif( order == 4 )
   ngp = 4;
  else
   ngp = 0;
  end
 elseif(type == QUAD_QUADRATURE)
  ngp = order*order;
 else
  ngp = 0;
 end

 if(ngp > 0)
  qp = zeros(2,ngp);
  w = zeros(1,ngp);
 end

 % Assemble arrays of gauss point coordinates and weights
 if(type == LINE_QUADRATURE)

  for i = 1:order
   gp(1,i) = qpts(i,order);
   gp(2,i) = 0.0;
   w(i)    = qwgt(i,order);
  end

 elseif(type == TRIA_QUADRATURE)

  if( order == 1 )
   %%% COMPLETE HERE 10 %%%
  elseif( order == 2 || order == 3 )
   %%% COMPLETE HERE 11 %%%
  elseif( order == 4 )
   %%% COMPLETE HERE 12 %%%
  end

 elseif(type == QUAD_QUADRATURE)

   %%% COMPLETE HERE 13 %%%

 end

end
