function fem
clear all;clc;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Standard finite element method FEM (1D)                                %
%                                                                         %
%    - d/dx (k(x) du/dx) + c(x) du/dx + b(x) u(x) = f(x)                  %
%                                                                         %
%  with boundary condition at x =0, and x = l                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finite Element Method
%            "Q.Fang, T.Tsuchiya, T.Yamamoto, Finite difference, finite 
%            element and finite volume methods applied to two-point 
%            boundary value problems, Journal of Computational and Applied 
%            Mathematics 139 (2002) 9-19
%
%            p(x)=exp(1-x), f(x)=1+exp(1-x), u(0)=u(1)=0;
%            Exact solution: u(x)=x(1-exp(x-1))
%            Max norm: FDM 8.24e-5,8.31e-7,8.31e-9
%                      FEM 6.35e-5,6.36e-7,6.39e-9
%                      FVM 3.18e-5,3.18e-7,3.18e-9
global nnode gk gf nelem u nodes  kind x
setint;     %Gaussian points and weights
propset;    %Input data
formkf;
applybc;
u=gk\gf;    % Solve the linear system
%[u,flag,er,iter]=bicgstab(gk,gf,1e-14,20000);

%outpts=[-1,-0.5,0,0.5,1];
%for nel=1:nelem,nel
%    for npt=1:5,[xx,uh,ux,sigh,sigx]=eval(nel,outpts(npt));du=ux-uh;
%        dsig=sigx-sigh;[xx,uh,ux,du,sigh,sigx,dsig],end
%end

[ue,uee]=exact(x);[enorm,sqnorm]=enorms;
plot(x,u);hold;plot(x,ue,'r');
%load femexp02.mat;
load femjournal1exp01.mat;
maxnorm=[maxnorm max(abs(u-ue))];L2norm=[L2norm sqnorm];

%maxnorm=[maxnorm max(abs(u-ue))]
%L2norm=[L2norm sqnorm]
save femjournal1exp01.mat maxnorm L2norm -append;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exact Solution                                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u,sig]=exact(x)
u=(x.*(1-exp(x-1))).';sig=(1-exp(x-1)-x.*exp(x-1)).';
%u=(4*x.^2.5.*(1-x).^2.5).';sig=(10*x.^(3/2).*(1-x).^(5/2)-10*x.^(5/2).*(1-x).^(3/2)).';
%u=(x-sinh(x)/sinh(1)).';sig=(1-cosh(x)/sinh(1)).';

function [xk,xc,xq,xf] = getmat(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        function returns the material parameters of the two points       %
%        boundary value problems                                          %
%                                                                         %
%              -(k(x) u_x)_x + c(x) u_x + q(x) u(x) = f(x)                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%xk=1;      xc=0;       xq=0;   xf=6;
%xk=1+x;    %xc=cos(x); xq=x*x; xf=sin(x);
%xf = (1+x)*sin(x)-cos(x)+cos(x)*cos(x)+x*x*sin(x);
%xk=1;xc=0;xq=1;xf=x;
xk=exp(1-x);    xc=0;   xq=0;   xf=1+exp(1-x);
%xk=exp(1-x);    xc=0;   xq=0;   xf=exp(1-x)*(x^(1.5)*(1-x)^(1.5)*(60-20*x)-15*x^(0.5)*(1-x)^(0.5)*(1-2*x+2*x^2));

function propset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT DATA FOR THE PROBLEM                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global nnode nelem gk gf x kind nodes nint kbc vbc
a=0;b=1;         % Boundary points for a<x<b
nnode=129;       % Number of nodes
kind_e=1;       % 1 for linear, 2 for quadratic, 3 for cubic
nelem=(nnode-1)/kind_e;
gk=zeros(nnode,nnode);gf=zeros(nnode,1);
x=linspace(a,b,nnode);
% Basis functions
for i=1:nelem
    kind(i)=kind_e;	  
    nint(i)=4;  % order of integration
    for j=1:kind(i)+1,nodes(j,i)=j+kind(i)*(i-1);end
end
% boundary conditions
kbc(1) = 1; % BC at x = a, 1: Dirichlet, 2: Neumann BC; 3: Mixed.
kbc(2) = 1; % BC at x = b,

%ua=sin(a);ub=sin(b); % Dirichlet BC.
ua=0;ub=0;

%uxa=-(1+a)*cos(a);uxb=(1+ b)*cos(b);% Neumann BC.

%alf=3;beta=4;% Mixed BC
%gamma=alf*sin(a)+beta*cos(a);
%uxma=-alf*(1+a)/beta;
%uaa=gamma/alf;.

%gammab=alf*sin(b)+beta*cos(b);
%uxmb=alf*(1+b)/beta;
%ubb=gammab/alf;

switch kbc(1)
    case 1,vbc(1,1)=ua;vbc(2,1) = 0;
        case 2,vbc(1,1)=uxa;vbc(2,1)=0;
            case 3,vbc(1,1)=uxma;vbc(2,1)=uaa;
end

switch kbc(2)
    case 1,vbc(1,2)=ub; vbc(2,2)=0;
        case 2,vbc(1,2)=uxb;vbc(2,2)=0;
            case 3,vbc(1,2)=uxmb;vbc(2,2)=ubb;
end

function setint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function evaluate                                                  %
%       x(:,i) the Gaussian points                                        %
%       w(:,i) the weights of quadrature                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global xi w
for i=1:4,[xi(1:i,i),w(1:i,i)]=gnodewt(i,-1,1);end

function formkf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      kind(i): The order of basis elements                               %
%      ninit(i): The order of Gaussian quadrature                         %      
%      x:        The coordinates of the donal points                      %
%      xi, w:    The Gaussian points and weights                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global xi w nelem kind nodes x nint
for nel=1:nelem,
    n=kind(nel)+1;
    i1=nodes(1,nel);
    i2=nodes(n,nel);
    i3=nint(nel);
    xic=xi(:,i3);  wc=w(:,i3);
    [ek,ef]=elem(x(i1),x(i2),n,i3,xic,wc);
    assemb(ek,ef,nel,n);
end

function [ek,ef] = elem(x1,x2,n,nl,xi,w)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function elem evaluates the stiffness matrix gk and the right         %
%   hand side gf.                                                         %
%  x1,x2 : Coordinates x1 and x2 of the ends of the element.              %
%  n     : Number of nodal points (and shape functions) in the element.   %
%  nl    :  Order of the integration rule of Gaussian formula, 1,2,3,4    %
%  xi(l) :  The value of the global coordinate at an integration point.   %
%  w(l)  :  The integration weights.                                      %
%  dx    :  dx/d\xi = *x2-x1)/2                                           %
%  psi(i),dps(i): The value of shape functions and their derivatives      %
%  xk,xc,xb,xf: values of the material properties at x_l obtained bt      % 
%  calling getmat(x).                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dx=(x2-x1)/2;
% Initialize element arrays
ef=zeros(n,1);ek=zeros(n,n);
% Begin integration point loop
for l=1:nl
    x = x1 + (1.0 + xi(l))*dx;
    [xk,xc,xb,xf] = getmat(x);
    [psi,dpsi] = shape(xi(l),n);
    for i=1:n,
       ef(i) = ef(i) + psi(i)*xf*w(l)*dx;
       for j=1:n,
          ek(i,j)=ek(i,j)+(xk*dpsi(i)*dpsi(j)/(dx*dx) ...
               +xc*psi(i)*dpsi(j)/dx+xb*psi(i)*psi(j) )*w(l)*dx;
       end
    end
end

function [psi,dpsi]=shape(xi,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        Function ``shape'' evaluates the values of the basis functions   %
%        and their derivatives at a point xi.                             %
%                                                                         %
%        n: The basis function. n=2, linear, n=3, quadratic, n=4, cubic.  %
%        xi: The point where the base function is evaluated.              %
%        Output:                                                          %
%        psi:  The value of the base function at xi.                      %
%        dpsi: The derivative of the base function at xi.                 %
%        Reference: Finite element. An introduction y E.Becker, G.Carey,  %
%        and J.Oden, Vol.1., pp. 95-96.                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch n
case 2,
    % Linear base function
    psi(1)=(1-xi)/2;
    psi(2)=(1+xi)/2;
    dpsi(1)=-0.5;
    dpsi(2)=0.5;
    return
    case 3,
        % quadratic base function
        psi(1) = xi*(xi-1)*0.5;
        psi(2) = 1-xi*xi;
        psi(3) = xi*(xi+1)*0.5;
        dpsi(1) =  xi-0.5;
        dpsi(2) = -2*xi;
        dpsi(3) = xi + 0.5;
        return
        case 4,
            % cubic  base function
            psi(1) = 9*(1/9-xi*xi)*(xi-1)/16;
            psi(2) = 27*(1-xi*xi)*(1/3-xi)/16;
            psi(3) = 27*(1-xi*xi)*(1/3+xi)/16;
            psi(4) = -9*(1/9-xi*xi)*(1+xi)/16;
            dpsi(1) = -9*(3*xi*xi-2*xi-1/9)/16;
            dpsi(2) = 27*(3*xi*xi-2*xi/3-1)/16;
            dpsi(3) = 27*(-3*xi*xi-2*xi/3+1)/16;
            dpsi(4) = -9*(-3*xi*xi-2*xi+1/9)/16;
            return
end

function assemb(ek,ef,nel,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       function assemb assemble the stiffness matrix gk and the right    %
%       hand side gf by looping through the nel-th elements.              %     
%                                                                         %
%       nel: The nel-th element.                                          %
%       n:   number of nodal points in the element, e.g. linear n=2;      %
%            quadrtic, n=3; cubic: n=4.                                   %
%       nodes:  nodal point numbers of nodes in the nel-th element,       %
%            nodes(1,nel), nodes(2,nel), ..., nodes(n,nel).               % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gk gf nodes
for i=1:n,
  ig = nodes(i,nel);                      
  gf(ig) = gf(ig) + ef(i);% Assemble global vector gf
  for j=1:n,
     jg = nodes(j,nel);
     gk(ig,jg)=gk(ig,jg)+ek(i,j);%Assemble global stiffness matrix gk
  end
end

function applybc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        function aplybc applies the boundary condition to modifies the   %
%        linear system obtained from the finite element method.           %
%        For the meanings of the parameter, see propset.m.                %
%        Fot the treatment of the BC, see Finite element. An introduction %
%        by E.Becker, G.Carey, and J.Oden, Vol.1., pp. 80-84, 107-109.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global nnode gk gf kbc vbc
if kbc(1) == 1,
    drchlta(vbc);	            % Dirichlet BC at x=a.
else
    if kbc(1) == 2,	            % Neumann DC at x=a.
       gf(1) = gf(1) + vbc(1,1);
    else				        % Mixed BC at x=a.
       gk(1,1) =  gk(1,1) + vbc(1,1);
       gf(1) = gf(1) + vbc(1,1)*vbc(2,1); 
    end
end

if kbc(2) == 1,		            % Dirichlet BC at x=b.
   drchltb(vbc);
else
   if kbc(2) == 2,              % Neumann DC at x=b.
      gf(nnode) = gf(nnode) + vbc(1,2); 
   else				            % Mixed BC at x=b.
      gf(nnode) = gf(nnode) + vbc(1,2)*vbc(2,2);
      gk(nnode,nnode) =  gk(nnode,nnode) + vbc(1,2);
   end
end

function drchlta(vbc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        function incorporates the Dirichlet BC at x=a to the discrete    %
%        linear system. Reference: Finite element. An introduction by     %
%        E.Becker, G.Carey, and J.Oden, Vol.1., pp. 107- 108              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global nnode gk gf
for i=1:nnode,
   gf(i) = gf(i) - gk(i,1)*vbc(1,1);
   gk(i,1) = 0;
   gk(1,i) = 0;
end
gf(1) = vbc(1,1);
gk(1,1) = 1;
   
function drchltb(vbc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        function incorporates the Dirichlet BC at x=b to the discrete    %
%        linear system. Reference: Finite element. An introduction by     %
%        E.Becker, G.Carey, and J.Oden, Vol.1., pp. 107- 108              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global nnode gk gf
for i=1:nnode,
   gf(i) = gf(i) - gk(i,nnode)*vbc(1,2);
   gk(nnode,i) = 0;
   gk(i,nnode) = 0;
end
gk(nnode,nnode)=1;
gf(nnode) = vbc(1,2);

function [enorm,sqnorm]=enorms
global nelem nodes kind x
sqnorm=0;enorm=0;
for i=1:nelem

    n1=nodes(1,i);    n=kind(i)+1;    n2=nodes(n,i);    wt=(x(n2)-x(n1))/50;
 
    xi=-1;
    for k=1:51
        [xx,uh,ux,sigh,sigx]=eval(i,xi);
        [xk,xc,xb,xf]=getmat(xx);
        if k==1 || k==51,w=wt/2;else w=wt;end
        xi=xi+0.02;
        sqnorm=sqnorm+w*(ux-uh)^2;
        enorm=enorm+(sigx-sigh)^2*w/xk+(ux-uh)^2*w*xb;
    end
end
sqnorm=sqrt(sqnorm);enorm=sqrt(enorm);

function [xx,uh,ux,sigh,sigx]=eval(nel,xi)
global kind nodes u x
n=kind(nel)+1;
[psi,dpsi]=shape(xi,n);
uh=0;sigh=0;
for i=1:n
    i1=nodes(i,nel);
    if i==1,x1=x(i1);end
    if i==n,x2=x(i1);end
    uh=uh+psi(i)*u(i1);
    sigh=sigh+dpsi(i)*u(i1);
end
dx=(x2-x1)*0.5;
xx=x1+(1+xi)*dx;
[xk,xc,xb,xf]=getmat(xx);
sigh=sigh*xk/dx;
[ux,sigx]=exact(xx);

function [x,w]=gnodewt(n,a,b)
[mu,v]=glnodewt(n);
w=0.5*(b-a)*v;
x=a+0.5*(mu+1)*(b-a);

function [x,w]=glnodewt(n)
b=(1:n-1)./sqrt(4*(1:n-1).^2-1);
J=diag(b,-1)+diag(b,1);
[V,D]=eig(J);[x,ix]=sort(diag(D));w=2*V(1,ix)'.^2;