function hpFEM_1D
%
%Beispiel einer Implementierung von hp-FEM in 1D
%
pmax = 4; 
mmax = 7; 
%create convergence plots for different values of p and uniform meshes 
exact_energy = pi; 
energy = zeros(mmax,pmax); 
h      = zeros(mmax,pmax); 
%
for p=1:pmax
 for m=1:mmax
  M = 2^m; 
  energy(m,p) = fem_1D(M,p)
  h(m,p) = 1.0/M; h
 end; 
end; 
error = sqrt(abs(exact_energy - energy)); 
figure(1)
loglog(h(:,1),error(:,1),'*-k',...
       h(:,2),error(:,2),'o--b',...
       h(:,3),error(:,3),'d:r',...
       h(:,4),error(:,4),'+-.m','LineWidth',2)
set(gca,'FontSize',14)
axis([1/200 1 10^(-7) 1])
xlabel('mesh size h')
ylabel('energy norm error')
title('-u^{\prime\prime} + u =f; solution u(x) = sin(x); uniform meshes')
legend('p=1','p=2','p=3','p=4','Location','SouthEast')
%
%check out sparsity pattern 
%
m = 4; 
mesh = [0:pi/m:pi]; 
p=1; 
[B,l] = assemble_neumann_problem(mesh,p); 
figure(2)
spy(B)
title('p=1; 4 elements')
p=2; 
[B,l] = assemble_neumann_problem(mesh,p); 
figure(3)
spy(B)
title('p=2; 4 elements')
p=5; 
[B,l] = assemble_neumann_problem(mesh,p); 
figure(4)
spy(B)
title('p=5 4 elements')
%
p=10; 
[B,l] = assemble_neumann_problem(mesh,p); 
figure(5)
spy(B)
title('p=10; 4 elements')

return
end
%-------------------------------------------------------------------------------
function energy = fem_1D(m,p)
%m = parameter for the mesh (here: uniform mesh)
%p = polynomial degree 
%
%numbering of the unknowns is such that the first and the 
%last unknown correspond to the left and right endpts 
%

g = -1;                   %value of the Neumann bc at the right endpt 

%define the mesh
xleft = 0; xright = pi; 
h = (xright-xleft)/m; 
mesh =[xleft:h:xright];   %uniform mesh

%compute the stiffness matrix and the load vector without boundary conditions 
%for the model problem - u^{\prime\prime} + u  = f
[B,l] = assemble_neumann_problem(mesh,p); 
N = size(B,1);            % get size of B
l(N) = l(N)+g;            % add the  Neumann boundary conditions 
                          % at the right endpoint

BD = B([2:N],[2:N]);      %drop the first row and column to enforce 
                          %homogeneous Dirichlet b.c. at left end point
lD = l([2:N]);            %also drop the first entry of load vector
u = BD\lD;                %solve 
energy = u'*lD;           %compute u^\top B u
return
end
%-------------------------------------------------------------------------------
function [B,l] = assemble_neumann_problem(mesh,p); 
%input: a mesh, i.e., a (sorted) list of nodes x1, x2,...,xM
%input: polynomial degree p
%
%
[T,N] = local_to_global_maps(mesh,p); 
%
B = sparse(N,N); 
l = zeros(N,1); 
for K = 1:length(mesh)-1
 xleft = mesh(K); xright = mesh(K+1); 
 vertex_coords = [xleft; xright]; 
 BK = element_stiffness_matrix(vertex_coords,p); 
 lK = element_load_vector(vertex_coords,p); 
 for i=1:p+1
  for j=1:p+1
    B(T(K,i),T(K,j)) = B(T(K,i),T(K,j))+BK(i,j); 
  end;
  l(T(K,i)) = l(T(K,i))+lK(i); 
 end;
end;  
return
end
%-------------------------------------------------------------------------------
function [T,N] = local_to_global_maps(mesh,p) 
%input: the mesh and the polynomial degree 
%output: matrix T, where T(K,i) gives the number of the global DOF 
%        corresponding to the function N_i on element K
%        N = total number of DOF 
%        recall: N_1, N_2 are linears, N_i for i \ge 3 are bubbles 
%        global numbering: for each element, first the left endpt, 
%                          then all bubbles, then the right endpt 
%
T = zeros(length(mesh)-1,p+1); 
N = 1; 
for K=1:length(mesh)-1 
  T(K,1) = N; 
  N = N +1; 
  for i=3:p+1
    T(K,i) = N ; 
    N = N +1; 
  end; 
  T(K,2) = N ; 
end;
return
end
%-------------------------------------------------------------------------------
function B = element_stiffness_matrix(vertex_coords,p) 
% vertex_coords(1:2) contains the two endpoints of the element
% p >= 1 polynomial degree to be employed 
%
h = vertex_coords(2) - vertex_coords(1); 
K = zeros(p+1,p+1); 
M = zeros(p+1,p+1); 
q = p+1;           %number of Gaussian quadrature points to be employed 
[x,w] = gauleg(q); %x = quadrature points on (-1,1); w = quadrature weights
x = x'; w = w';    %gauleg returns column vectors, we want row vectors 
PHI = N(x,p+1);           %evaluate the shape fcts N_i at the quadrature points 
GRAD_PHI = grad_N(x,p+1); %evaluate N_i' at the quadrature points 


   K = GRAD_PHI*diag(w)*GRAD_PHI'/(h/2); %realizes the numerical integration of 
                                         %\int_{-1}^1 N_i^\prime N_j^\prime/(h/2)
   M = PHI*diag(w)*PHI'*h/2;             %numerical integration of 
                                         %\int_{-1}^1 N_i N_j*(h/2)
B = K+M;                                                           
%
%---comment
% the quadratures to compute K,M are written very compactly. multiplying out, 
% they amount to  setting  M(i,j) = 
% \sum_{k=1}^{length(x)} w(k) PHI(i,k)PHI(j,k)*h/2 for all i,j =1,\ldots,p+1
% similarly for the matrix K
%---endcomment
return
end
%-------------------------------------------------------------------------------
function l = element_load_vector(vertex_coords,p) 
% vertex_coords(1:2) contains the two endpoints of the element
% p >= 1 polynomial degree to be employed 
%
h = vertex_coords(2) - vertex_coords(1); 
l = zeros(p+1,1); 
q = p+1;           %number of Gaussian quadrature points to be employed 
[x,w] = gauleg(q); %x = quadrature points on (-1,1); w = quadrature weights
x = x'; w = w';    %gauleg returns column vectors, we want row vectors 
PHI = N(x,p+1);                     %evaluate the shape fcts N_i at the quadrature points 
F = f(h/2*(x+1)+vertex_coords(1));  %evaluate the right-hand side at the quadrature points 

l = PHI*diag(w)*F'*h/2;             %\int_{-1}^1 f N_i*h/2
%
%---comment
% the quadrature to compute l,M is written very compactly. Multiplying out, the definition of l 
% amounts to  setting l(i) = \sum_{k=1}^{length(x)} w(k) PHI(i,k) * F(k) *h/2    for all i
% similarly for the matrix K
%---endcomment
return
end
%-------------------------------------------------------------------------------
function [U] =N(x,n)

%
%  Calculates the first n trial functions on [-1,1]; 
%  we assume implicitly that the vector x is a row vector
%
%  N_1 = (1-x)/2
%  N_2 = (1+x)/2
%  N_n = \sqrt{(2n-3)/2}\int_{-1}^x L_{n-2}(t) dt     if n \ge 3
%
  U(1,:) = (1-x)/2;
  U(2,:) = (1+x)/2;
if n >2 
[P] = legtable(x,n);
  for i = 3:n
      ip=i-2;
        U(i,:) = sqrt((2*i-3)/2)*1/(2*ip+1)*(P(ip+2,:)-P(ip,:));
  end;
end
return
end
%-------------------------------------------------------------------------------
function [U] =grad_N(x,n)

%
%  Calculates the derivatives of the first n trial functions on [-1,1]; 
%  we assume implicitly that the vector x is a row vector
%
%  N_1 = (1-x)/2
%  N_2 = (1+x)/2
%  N_n = \sqrt{(2n-3)/2}\int_{-1}^x L_{n-2}(t) dt     if n \ge 3
%
  U(1,:) = -1/2*ones(size(x));
  U(2,:) =  1/2*ones(size(x));
if n >2 
[P] = legtable(x,n);
  for i = 3:n
      ip=i-2;
        U(i,:) = sqrt((2*i-3)/2)*P(ip+1,:);
  end;
end
return
end
%-------------------------------------------------------------------------------
function [P] = legtable(x,m)
%
% input:  row vector x \subset (-1,1) of points
%         m the maximal order of the Legendre polynomials
% output: a matrix P=P(m+1,length(x)) containing the values of the
%         Legendre polynomials at the points x
%
  l=length(x);
  P=ones(m+1,l);
  P(2,:)=x;
  for i=2:m
    P(i+1,:)=((2*i-1)*x.*P(i,:)-(i-1)*P(i-1,:))/i;
  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
%-------------------------------------------------------------------------------
function y = f(x)
%right-hand side function f 
y = 2*sin(x); 
%exact solution: u(x) = sin(x) on (0,pi);    hence:  -u''+u = 2*sin(x)
%exact energy: pi
return
end
%-------------------------------------------------------------------------------