function [xp,w] = FEM4x4Gausspts
% [xp,w] = FEM4x4Gausspts (actually 4x4x4)
% Computed by table lookup from coefficients calculated with
% Mathematica using the following computation:
% 
% n = 4;
% f[x_]=LegendreP[n,x];
% g[x_]=D[f[x],x];
% z=Solve[f[x]==0,x];
% r=x/.z;
% w[x_]= 2/((1 - x^2) g[x]^2);
% pts=N[r,20];
% wei=N[w[r],20];
% See for example:
% Th. J.R. Hughes, The finite element method Sec 3.8 pp 143,
% For matlab:
% f = inline('3/8 - (15*x.^2)/4 + (35*x.^4)/8');
% g  = inline('(35*x.^3-15*x)/2');

  
  r =[-0.3399810435848562648;
      0.3399810435848562648;
      -0.86113631159405257522;
      0.86113631159405257522];
  
  wt =[  0.65214515486254614263;
	0.65214515486254614263;
	0.34785484513745385737;
	0.34785484513745385737];
  
  n=4;
  nnn=n^3;
  xp=zeros(nnn,3);
  w = zeros(nnn,1);
  
  kk=0;
  for i=1:n,
    for j=1:n,
      for k=1:n,
	kk=kk+1;
	xp(kk,:) = [r(i),r(j),r(k)];
	w(kk) = wt(i)*wt(j)*wt(k);
      end;
    end;
  end;