%--------------------------------------
function [x,w]= create_1d_hp_quadrature(L,p)
%p = number of Gauss points per element; L = number of layers for quadrature
%creates quadrature formula for \int_0^1 f(x) with geometric refinement towards x = 0
%assumes: f will be evaluated at the Gauss point for each subinterval 
%
integral = 0; 
sigma = 0.15; 
[xx,ww]=gauleg(p); 
x = zeros((L+1)*p,1); 
w = zeros((L+1)*p,1); 
%first interval of the geometric mesh
xleft = 0; xright = sigma^L; m = 0.5*sigma^L; hh = sigma^L; 
x(1:p,1) = 0.5*hh*(1+xx); 
w(1:p,1) = 0.5*hh*ww; 
%
%do the remaining elements of the geometric mesh 
for l=1:L
  xleft = sigma^(L-l+1); xright = sigma^(L-l); hh = sigma^(L-l)*(1-sigma); 
x(l*p+1:(l+1)*p,1) = xleft+0.5*hh*(1+xx); 
w(l*p+1:(l+1)*p,1) = 0.5*hh*ww; 
end; 
return 
end
%---------------------------------------
function [x,w]=gauleg(n)

x=zeros(n,1);
w=zeros(n,1);
  m=(n+1)/2;
  xm=0.0;
  xl=1.0;
  for i=1:m
    z=cos(pi*(i-0.25)/(n+0.5));
    while 1
      p1=1.0;
      p2=0.0;
      for j=1:n
        p3=p2;
        p2=p1;
        p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
      end
      pp=n*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
      if (abs(z-z1)<eps), break, end
    end
    x(i)=xm-xl*z;
    x(n+1-i)=xm+xl*z;
    w(i)=2.0*xl/((1.0-z*z)*pp*pp);
    w(n+1-i)=w(i);
  end
return
end