function [sol,xnod] = FEM1D (a, b, xnel)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEM1D
%% 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
%% OUTPUT:
%% <sol> is the computed numerical solution at nodes
%% <xnod> position of all nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 01/15/2008
%% © 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. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
   for i=1:nels-1 hel(i) = xnel(i+1)-xnel(i); end;

else 

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

end; 


%% set up uniform order of elements = 1 + degree of polynomial =
%% = number of degrees of freedom 

ord = zeros(nels,1) + 2;  %% type of elements: linear 
maxord = max(ord);

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

%% derive 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 i = 1:nels 
    for j = 1:ord(i)
        nod(i,j) = n;
        if j == 1 myel(n,2) = i;
        elseif j == ord(i) myel(n,1) = i;
        else myel(n,1) = i; myel(n,2) = i;
        end;
        if j ~= ord(i) n = n+1; end
    end; 
end;

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

%%%%%%%%%%%%%%%%%%%%%%%%% set up numerical integration

%% 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);
end;

%%%%%%%%%%%%%%%%%%%%%%%%% 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,2)); %% left and right endpoint 
    dx = (x2-x1)/2.;                            %% Jacobian of transformation 
  

    %% 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);          %% 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

    %% add the computed element stiffness matrix and load vector to 
    %% the global matrix and vector

    rhsf (nod(el,:)) = rhsf(nod(el,:)) + f(:);
    mat (nod(el,:),nod(el,:)) = mat(nod(el,:),nod(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);

%%%% compute FD solution for demonstration
xin_node = xnod(2:end-1);
fd_rhs = h*rhsfun(xin_node);
fd_sol = mat(freenodes,freenodes) \fd_rhs;
%%%%

xplot = linspace(a,b,xnel*10);
plot(xnod,exfun(xnod),'k',xnod,sol,'r-o');
legend('Exact solution','Numerical solution (linear FE)');
pause,
plot(xplot,exfun(xplot),'k',xnod,sol,'r-o');
legend('Exact solution','Numerical solution (linear FE)');
%%% display FD solution
pause,
plot(xplot,exfun(xplot),'k',xnod,sol,'r-o',xin_node,fd_sol,'b*');
legend('Exact solution','Numerical solution (linear FE)','Numerical solution (FD)');
%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% end of algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y,dy] = shape(x,n)
%% shape function on reference element (-1,1)
%% n = 2: linear
%% n = 3: quadratic (must be coded)

if n == 2

   y  (1,:) = .5.*(1.-x);
   y  (2,:) = .5.*(1.+x);
   dy (1,:) = -.5;
   dy (2,:) = .5;

end

function y = afun(x)    %%%%%%%%%%%%%% coefficient a(x)
y = 1.0;

function y = rhsfun(x)  %%%%%%%%%%%%% rhs function
%    y = 2.0*ones(size(x));
    y = pi*pi*sin(pi*x);

function y = exfun(x)   %%%%%%%%%%%%% exact solution: provides ua,ub
%    y=x.*(1-x);
     y = sin(pi*x);
