function [sol,xnod] = fem1d_2017 (a, b, xnel, myorder)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEM1D_2017
%% for MTH 655
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This function shows how to set up a finite element solution
%% of a two-point BVP 
%% -(au')' =f, u(a)=ua, u(b)=ub 
%% using Galerkin linear Finite Elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function sol = solfem (a,b,xnel)
%% a,f,ua,ub are provided as functions in code
%% numerical integration is used in computation of stiffness matrix and rhs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT:
%% <a,b> = endpoints of interval 
%% <xnel> =  the number or location of elements
%% if xnel is one number (scalar) then xnel = number of elements, 
%% if xnel is a vector, then it contains position of first node of each element
%% <myorder> =  degree of polynomials per element
%% if myorder is one number (scalar) then it is uniform over all elements, 
%% if myorder is a vector, then it must be as long as the vector of elements
%% OUTPUT:
%% <sol> is the computed numerical solution at nodes
%% <xnod> position of all nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 01/15/2008 
%% M. Peszynska, modified for MTH 655/9 Winter 2017
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% references: 
%% ACF software in Alberty, Carstensen, Funken, 1999
%% Becker, Carey, Oden "Introduction to Finite Elements" Vol I-VI. 
%% Class handouts for MTH 655/9, Winter 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% how to run:
%% uniform grid/element order 
%% fem1d_2017(0,1,30,1);
%% nonuniform grid: 
%% fem1d_2017(0,1,[0,0.18,0.82],[1 2 1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
set(0,'DefaultLineLineWidth',3);
%%%%%%%%%%%%%%%%%%%%%%%%% set up grid 
if size(xnel) == 1   %% uniform grid
    
   nels = xnel;
   hel = (b-a)/nels; 
   xel = a : hel : b-hel;   

elseif size(xnel,1) == 1
   nels = size(xnel,2);xel=xnel;
   hel=zeros(nels-1,1);
   for k=1:nels-1, hel(k) = xnel(k+1)-xnel(k); end;

else 

   error('wrong number of elements/grid');

end; 
%% set up the  order of elements 
%% = 1 + degree of polynomial = number of degrees of freedom 
%% for ex.: for linears, use two degrees of freedom for each element
ord = zeros(nels,1) + myorder' + 1;  
maxord = max(ord);

%% number  of nodes (including boundary values)
nnodes = sum(ord-1)+1;

%% derive the global indexing of nodes: 
%% nod(i,1) is the global number of j'th node in element i
nod = zeros(nels,maxord); myel = zeros(nnodes,2);
n = 1;
for k = 1:nels 
    for j = 1:ord(k)
        nod(k,j) = n;
        if j == 1, myel(n,2) = k;
        elseif j == ord(k), myel(n,1) = k;
        else myel(n,1) = k; myel(n,2) = k;
        end
        if j ~= ord(k), n = n+1; end
    end; 
end
%% uncomment to see the way the nodes and elements are numbered
%nod, pause
%myel,pause

%% xnod (i=1..nnodes): coordinates of node i 
xnod = zeros(nnodes,1);
for k=1:nels-1 
    %% h is the size of the element. hi is the size of subdivision
    h = xel(k+1)-xel(k); hi = h/(ord(k)-1);
    for j=1:ord(k), xnod (nod(k,j)) = xel(k) + hi*(j-1);end;
end;  
k = nels; h = b-xel(k); hi=h/(ord(k)-1);
for j=1:ord(k), xnod (nod(k,j)) = xel(k) + hi*(j-1);end;

hmax = max(hel);
%%%%%%%%%%%%%% set up numerical integration for stiffness matrix
%% set up quadrature parameters on the reference element (-1,1)
%% set up number of integration points nw, nodes xw, and weights w
if maxord == 1      %% exact for linears
      nw = 1;
      xw(1) = 0.;  
      w(1) = 2.;
elseif maxord == 2  %% exact for cubics
      nw = 2;
      xw(1) = -1/sqrt(3); xw(2) = -xw(1);
      w(1) = 1; w(2) = 1;
elseif maxord == 3  %% exact for polynomials of degree 5
      nw = 3;
      xw(1) = -sqrt(3./5.); xw(2)=0.; xw(3) =- xw(1);
      w(1) = 5./9.; w(2)=8./9.; w(3)=w(1);
else 
    error('Error: no code implemented for the desired order\n');
end;

%%%%%%%%%%%%% set up matrix and rhs of linear system 
%% set up can be easily modified to include higher order elements
%% as well as 
mat = sparse(nnodes,nnodes);  rhsf = zeros(nnodes,1); 

for el = 1:nels
    x1 = xnod(nod(el,1)); x2 = xnod(nod(el,ord(el))); %% left and right endpoint 
    dx = (x2-x1)/2.;                            %% Jacobian of transformation 

    %% uncomment to see the geometrical information about the nodes and elements
    %fprintf('Element %d with %d nodes\n',el,ord(el));
    %for k=1:ord(el)
    %    fprintf('.. %g ..',xnod(nod(el,k)));
    %end
    %fprintf('\n');

    %% compute element stiffness matrix and load vector
    m = zeros(ord(el),ord(el));                 %% element stiffness matrix
    f = zeros(ord(el),1);                       %% element load vector 

    for l = 1:nw
        x = x1 + (1 + xw(l))*dx;                %% x runs in true element, 
                                                %% xw runs in reference element
        [psi,dpsi] = shape(xw(l),ord(el));      %% calculations on ref.element
        %
        aval = feval(@afun,x);
        fval = feval(@rhsfun,x);
     
        f = f + fval * psi * w(l)*dx;  
        m = m + aval * (dpsi)*(dpsi')/dx/dx * w(l)*dx;
                
    end
    %% uncomment to see the local stiffness matrix and local load vector 
    %m, pause,
    %f, pause,
    
    %% add the computed element stiffness matrix and load vector to 
    %% the global matrix and vector

    rhsf (nod(el,1:ord(el))) = rhsf(nod(el,1:ord(el))) + f(1:ord(el));
    mat (nod(el,1:ord(el)),nod(el,1:ord(el))) = mat(nod(el,1:ord(el)),nod(el,1:ord(el))) + m;
    
end 

%% impose Dirichlet boundary conditions
%% eliminate known values from the system

freenodes = 2:nnodes-1;
sol = zeros(nnodes,1); 
sol (1) = exfun(a); sol (nnodes) = exfun(b);

rhsf = rhsf - mat*sol;

%%%%%%%%%%%%%%%%%%%%%%%%% solve linear system
sol (freenodes) = mat(freenodes,freenodes) \ rhsf (freenodes);

%%%%%%%%%% Now, postprocess: plot, calculate the error etc.
xplot = linspace(a,b,nels*10+1);
plot(xplot,exfun(xplot),'k');
hold on;

%%%% PLOTTTING element by element
for el = 1:nels
    x1 = xnod(nod(el,1)); x2 = xnod(nod(el,ord(el))); %% left and right endpoint 
    dx = (x2-x1)/2.;                            %% Jacobian of transformation 
    
    xpp = linspace(x1,x2,10);                   %% set up 10 points in the true element for plotting
    xww =(xpp-x1)/dx-1;                         %% xww runs in reference element
    ypp = 0*xpp;
    
    for j=1:length(xpp)
        [psi,dpsi] = shape(xww(j),ord(el));     %% evaluate all basis function in reference element
        uhval = 0;                              %% compute  
        for k=1:ord(el)
            mynum = nod(el,k);
            uhval = uhval + sol(mynum)*psi(k);
        end
        ypp(j)=uhval;
    end
    
    plot(xpp,ypp,'r-');
end;  
hold off;

legend('Exact solution',...
    sprintf('FE solution of max order =%d',maxord-1),...
    'Location','Best');

%%% 

%%% ERROR COMPUTATION
l2err = 0; h1err = 0;

for el = 1:nels
    x1 = xnod(nod(el,1)); x2 = xnod(nod(el,ord(el))); %% left and right endpoint 
    dx = (x2-x1)/2.;                            %% Jacobian of transformation 

 
    for l = 1:nw
        x = x1 + (1 + xw(l))*dx;                %% x runs in true element, 
                                                %% xw runs in reference element
        [psi,dpsi] = shape(xw(l),ord(el));          %% calculations on ref.element
     
        uval = feval(@exfun,x);
        duval = feval(@dexfun,x);
        
        uhval = 0; duhval = 0;
        for k=1:ord(el)
            mynum = nod(el,k);
            uhval = uhval + sol(mynum)*psi(k);
            duhval = duhval + sol(mynum)*dpsi(k)/dx;
        end
        l2err = l2err + (uval-uhval).^2*w(l)*dx;
        h1err = h1err + (duval-duhval).^2*w(l)*dx;
                      
    end
    
end 

h1err = h1err + l2err;
h1err = sqrt(h1err);
l2err = sqrt(l2err);
fprintf('hmax=%g l2err =%g h1err =%g\n',hmax,l2err,h1err);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% end of algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

function [y,dy] = shape(x,n)
%% shape function on reference element (-1,1)
%% n = 2: linear
%% n = 3: quadratic 
if n == 2
   y  (1,:) = .5.*(1.-x);
   y  (2,:) = .5.*(1.+x);
   dy (1,:) = -.5;
   dy (2,:) = .5;
elseif n == 3 
      error('You must code this first!');
%%    y(1,:) = YOU JUST CODE THIS
%%    y(2,:) = 1-x.^2;
%%    y(3,:) = YOU JUST CODE THIS
%%    dy(1,:) = YOU JUST CODE THIS
%%    dy(2,:) = -2.*x;
%%    dy(3,:) = YOU JUST CODE THIS
else
    error('Implementation for n=%d not provided\n',n);
end
end

function y = afun(x)    %%%%%%%%%%%%%% coefficient a(x)
    y = 1.0;
end
function y = rhsfun(x)  %%%%%%%%%%%%% rhs function
    y = pi*pi*sin(pi*x);
end
function y = exfun(x)   %%%%%%%%%%%%% exact solution: provides ua,ub
     y = sin(pi*x);
end  
function y = dexfun(x)   %%%%%%%%%%%%% derivative of the exact solution: provides ua,ub
     y = pi*cos(pi*x);
end

