function OneD(Order,NumE)
%1D Code for ELASTOSTICS/HEAT EQUATION  (using linear shape functions only)
%AUTHOR: S. Srivastava, University of Twente, NL
%Run as OneD(1,1), OneD(1,10), OneD(2,10) ..etc.
clear global;
format compact;     %CONTROLS NUMBER OF SIGNIFICANT DIGITS 

global KGlobal FGlobal Total_Elements Coord_Table L2GNodeTable L2GDofTable;  %GLOBAL VARIABLES
global Area Youngs Traction L DeltaX TSF;    %CODE PARAMETERS

L=1; %1m rod
Area=1e-4; %m^2
Youngs=70e9; %N/m^2     %MATERIAL AND GEOMETRIC PROPS.

Total_Elements=NumE;
Traction=1000; %N/m   %APPLIED LOAD

%Traction=Traction/(Youngs*Area);    %for solving -del(U)=f/(E*A)

%Discretize
counter=1;
DeltaX=L/(Total_Elements);

%Initaialize variables
L2GNodeTable=0;
L2GDofTable=0;
KGlobal=0;
FGlobal=0;

c=1; %A COUNTER
for i=1:Total_Elements
    for j=1:2    %ONLY 2Nodes  
        L2GNodeTable(i,j)=counter;
        Coord_Table(i,j)=(c-1)*DeltaX;
        counter=counter+1;
        c=c+1;
    end
    counter=counter-1;
    c=c-1;
end

counter=1;
for i=1:Total_Elements
    for j=1:2   %ONLY LINEAR(2) SHAPE ELEMENTS  
        L2GDofTable(i,j)=counter;
        counter=counter+1;
    end
    counter=counter-1;
end
counter=counter+1;
for i=1:Total_Elements
    for j=3:Order+1   %Higher order (bubble) shape functions
        L2GDofTable(i,j)=counter;
        counter=counter+1;
    end
end
%Generate Shape Funcitons
%SHAPE FUNCTIONS STORE POOLYNOMIAL COEFFICIENTS
for i=1:Order+1
    N(:,i)=putzero(N1D(i),Order+1); %fill zeros for zero coeffs of shap funcs
end
TSF=Order+1;   %TOTAL NUMBER OF SHAPE FUNCTIONS on each element
n=Order+1;
Gauss=GLQuads(n)      %GAUSS-LEGENDRE QUADRATURE ARRAY

%MAIN PROGRAM STARTS       
for E=1:Total_Elements   %LOOP OVER ELEMENTS
    KElement=0; %use zeros(..) for efficiency
    FElement=zeros(TSF,1);
    E
    for i=1:n           %LOOP OVER GAUSS POINTS
       J=Jacob(E,Gauss(1,i));       %EVALUATE JACOBIAN
       DJ=det(J);
       KElement=KElement+Gauss(2,i)*ButCBu(E,J,N,Gauss(1,i))/DJ;   %ELEMENT STIFFNESS MATRIX
    
       for j=1:TSF
            FElement(j,1)=FElement(j,1)+Gauss(2,i)*Evaluate(N(:,j),Gauss(1,i))*DJ;     %ELEMENT FORCE OR R.H.S VECTOR
       end
       
    end%END GAUSS
    
    FElement=Traction*FElement;
    AssembleK(E,KElement);
    AssembleF(E,FElement);      %ASSMEBLE ELEMENT STIFFNESS AND VECTOR INTO GLOBAL MATRICES
    KElement
end%END ELEMENTS
KGlobal
FGlobal'
RANK=rank(KGlobal)
SIZEK=size(KGlobal)

disp('Applying boundary conditions.....');
BNodes=[1];    %BOUNDARY CONDITION NODES
ApplyBC(BNodes); %Remove rows column corresponding to fixed DOF's

KGlobal
FGlobal'
RANK=rank(KGlobal)
SIZEK=size(KGlobal)
SIZEF=size(FGlobal)

Sol=KGlobal\FGlobal'; %SOLVE THE SYSTEM OF EQUATIONS USING MATLAB'S '\' OPERATOR

File=strcat('1D','.sol');
dlmwrite(File,Sol,',');
File=strcat('1D','.dof');        %STORE RESULTS ON DISK AS DELIMETED VALUES OF COEFFICIENTS
dlmwrite(File,L2GDofTable,',');

disp('Solution has been computed and stored on the disk, Use Postprocessing routines to view results-BYE.');
%CALL POSTPROCESSING ROUTINES.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%***END***%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%OTHER FUNCTIONS
function [Ans]=ButCBu(E,J,N,x)
%1D stiffness integrand
global Area Youngs TSF;
Bu=zeros(TSF,TSF);
J=inv(J);
for i=1:TSF
    Bu(1,i)=Evaluate(Derivative(N(:,i)),x);
end
Ans=Bu'*(Area*Youngs)*Bu;    %<<<<<< THIS QUANTITY IS A MATRIX!

function [answer]=Derivative(poly) 
%simply differentiates a given polynomial(analytically)
n=length(poly);
a(1)=0;
for i=1:n-1
    if poly(i)~=0
        answer(i+1)=poly(i)*(n-i);
    else answer(i+1)=0;
    end
end
answer=answer';

function [coff]=N1D(order)
%1D SHAPE FUNCITONS
if order==1
    coff=[-0.5,0.5];
    elseif order==2
            coff=[0.5,0.5];
    else coff=(1/sqrt(4*order-6))*(polyadd(Legendre_poly(order-1),-Legendre_poly(order-3)));         
end

function [Lcoff]=Legendre_poly(order)
%USED BY N1D, CREATES POLYNOMIALS (LEGENDRE)
x=[1,0];
if order==0
    Lcoff=1;
elseif order==1
        Lcoff=[1,0];
else Lcoff=(1.0/order)*polyadd((2*order-1)*conv(x,Legendre_poly(order-1)),(1-order)*Legendre_poly(order-2));
end

function [J]=Jacob(E,zeta)
%1D Jacobian
global Coord_Table;
J=0;
for j=1:2
   f=N1D(j);
   J=J+Coord_Table(E,j)*Evaluate(Derivative(f),zeta);
end

function [answer]=Evaluate(func,x)%evaluate shape function value at a point in master element(1D only)
if func==0
    answer=0;
else
    if x>=-1 & x<=1
        answer=polyval(func,x);
    else error('Error:CANNOT Evaluate Shape Function At this Point(x<=-1 >=1)..');
    end
end

function  [NKe]=ApplyBC(BNodes)
%IMPOSE ESSENTIAL BC(at x==0)
global KGlobal FGlobal Total_Elements L2GDofTable TSF;
for i=1:length(BNodes)
    Gi=L2GDofTable(1,1);
    KGlobal(Gi,:)=0;
    KGlobal(:,Gi)=0;
    KGlobal(Gi,Gi)=1;
    FGlobal(Gi)=0;
end

%ASSMEBLY ROUTINES
function AssembleK(E,Ke)
%assemble stiffness matrix add Ke to Kglobal
global KGlobal L2GDofTable;

Dof=length(Ke);
for i=1:Dof
    Gr=L2GDofTable(E,i);
    for j=1:Dof
        Gc=L2GDofTable(E,j);
        [R,C]=size(KGlobal);
        
        if Gr==0 | Gc==0
            Error('Zero GDof in AssembleK.....');    
        end
        if Gr>R | Gc>C
            KGlobal(Gr,Gc)=0;
        end
        KGlobal(Gr,Gc)=KGlobal(Gr,Gc)+Ke(i,j);
    end
end

function AssembleF(E,Fe)
global FGlobal L2GDofTable;
%Assemble load Vector
Dof=length(Fe);

for i=1:Dof
    Gi=L2GDofTable(E,i);
    if Gi==0
        error('Zero Gdof in AssembleF....');
    end
    l=length(FGlobal);
    if Gi>l
        FGlobal(Gi)=0;
    end
    FGlobal(Gi)=FGlobal(Gi)+Fe(i);
end

function[poly]=polyadd(poly1,poly2)
if length(poly1)<length(poly2)
  short=poly1;
  long=poly2;
else
  short=poly2;
  long=poly1;
end
mz=length(long)-length(short);
if mz>0
  poly=[zeros(1,mz),short]+long;
else
  poly=long+short;
end

function [p]=putzero(m,n)
l=length(m);
p=zeros(n,1);
for i=1:l
    p((n-l)+i)=m(i);
end

function [Gauss]=GLQuads(n)
%RETURNS Gauss-Legendre Integration POINTS AND WEIGHTS (QUADRATURE) (row 1= POINTS) & (row 2= WEIGHTS) 
switch n
  case 1
    Gauss(1,1)=0;
    Gauss(2,1)=2;
   
  case 2
      Gauss(1,1)=-0.57735;
      Gauss(1,2)=0.57735;
      %weight
      Gauss(2,1)=1; 
      Gauss(2,2)=1;
      
  case 3     
      Gauss(1,1)=-0.7746;
      Gauss(1,2)=0;
      Gauss(1,3)=0.7746;
      
      Gauss(2,1)=0.5555555;
      Gauss(2,2)=0.8888888;
      Gauss(2,3)=0.5555555;
      
  case 4
      Gauss(1,1)=-0.861136311594;
      Gauss(1,2)=-0.339981043585;
      Gauss(1,3)=0.339981043585;
      Gauss(1,4)=0.861136311594;
      
      Gauss(2,1)=0.347854845137;
      Gauss(2,2)=0.652145154863;
      Gauss(2,3)=0.652145154863;
      Gauss(2,4)=0.347854845137;
      
  case 5      
      Gauss(1,1)=-0.906179845939;
      Gauss(1,2)=-0.538469310106;
      Gauss(1,3)=0;
      Gauss(1,4)=0.538469310106;
      Gauss(1,5)=0.906179845939;
      
      Gauss(2,1)=0.236926885056;
      Gauss(2,2)=0.478628670499;
      Gauss(2,3)=0.568888888889;
      Gauss(2,4)=0.478628670499;
      Gauss(2,5)=0.236926885056;
      
  case 6    
      Gauss(1,1)=-0.932469514203;
      Gauss(1,2)=-0.661209386466;
      Gauss(1,3)=-0.238619186083;
      Gauss(1,4)=0.238619186083;
      Gauss(1,5)=0.661209386466;
      Gauss(1,6)=0.932469514203;
      
      Gauss(2,1)=0.171324492379;
      Gauss(2,2)=0.360761573048;
      Gauss(2,3)=0.467913934573;
      Gauss(2,4)=0.467913934573;
      Gauss(2,5)=0.360761573048;
      Gauss(2,6)=0.171324492379;
      
  case 7
      Gauss(1,1)=-0.949107912343;
      Gauss(1,2)=-0.741531185599;
      Gauss(1,3)=-0.405845151377;
      Gauss(1,4)=0;
      Gauss(1,5)=0.405845151377;
      Gauss(1,6)=0.741531185599;
      Gauss(1,7)=0.949107912343;
      
      Gauss(2,1)=0.129484966169;
      Gauss(2,2)=0.279705391489;
      Gauss(2,3)=0.381830050505;
      Gauss(2,4)=0.417959183673;
      Gauss(2,5)=0.381830050505;
      Gauss(2,6)=0.279705391489;
      Gauss(2,7)=0.129484966169;
      
  case 8
      Gauss(1,1)=-0.960289856498;
      Gauss(1,2)=-0.796666477414;
      Gauss(1,3)=-0.525532409916;
      Gauss(1,4)=-0.183434642496;
      Gauss(1,5)=0.183434642496;
      Gauss(1,6)=0.525532409916;
      Gauss(1,7)=0.796666477414;
      Gauss(1,8)=0.960289856498;
      
      Gauss(2,1)=0.10122853629;
      Gauss(2,2)=0.222381034453;
      Gauss(2,3)=0.313706645878;
      Gauss(2,4)=0.362683783378;
      Gauss(2,5)=0.362683783378;
      Gauss(2,6)=0.313706645878;
      Gauss(2,7)=0.222381034453;
      Gauss(2,8)=0.10122853629;
      
  case 9
      Gauss(1,1)=-0.968160239508;
      Gauss(1,2)=-0.836031107327;
      Gauss(1,3)=-0.613371432701;
      Gauss(1,4)=-0.324253423404;
      Gauss(1,5)=0;
      Gauss(1,6)=0.324253423404;
      Gauss(1,7)=0.613371432701;
      Gauss(1,8)=0.836031107327;
      Gauss(1,9)=0.968160239508;
      
      Gauss(2,1)=0.0812743883616;
      Gauss(2,2)=0.180648160695;
      Gauss(2,3)=0.260610696403;
      Gauss(2,4)=0.31234707704;
      Gauss(2,5)=0.330239355001;
      Gauss(2,6)=0.31234707704;
      Gauss(2,7)=0.260610696403;
      Gauss(2,8)=0.180648160695;
      Gauss(2,9)=0.0812743883616;
      
 case 10
      Gauss(1,1)=-0.973906528517;
      Gauss(1,2)=-0.865063366689;
      Gauss(1,3)=-0.679409568299;
      Gauss(1,4)=-0.433395394129;
      Gauss(1,5)=-0.148874338982;
      Gauss(1,6)=0.148874338982;
      Gauss(1,7)=0.433395394129;
      Gauss(1,8)=0.679409568299;
      Gauss(1,9)=0.865063366689;
      Gauss(1,10)=0.973906528517;
      
      Gauss(2,1)=0.0666713443087;
      Gauss(2,2)=0.149451349151;
      Gauss(2,3)=0.219086362516;
      Gauss(2,4)=0.26926671931;
      Gauss(2,5)=0.295524224715;
      Gauss(2,6)=0.295524224715;
      Gauss(2,7)=0.26926671931;
      Gauss(2,8)=0.219086362516;
      Gauss(2,9)=0.149451349151;
      Gauss(2,10)=0.0666713443087;
      
 case 11
      Gauss(1,1)=-0.978228658146;
      Gauss(1,2)=-0.887062599768;
      Gauss(1,3)=-0.730152005574;
      Gauss(1,4)=-0.519096129207;
      Gauss(1,5)=-.269543155952;
      Gauss(1,6)=0;
      Gauss(1,7)=0.269543155952;
      Gauss(1,8)=0.519096129207;
      Gauss(1,9)=0.730152005574;
      Gauss(1,10)=0.887062599768;
      Gauss(1,11)=0.978228658146;

      Gauss(2,1)=0.0556685671162;
      Gauss(2,2)=0.125580369465;
      Gauss(2,3)=0.186290210928;
      Gauss(2,4)=0.233193764592;
      Gauss(2,5)=0.26280454451;
      Gauss(2,6)=0.272925086778;
      Gauss(2,7)=0.26280454451;
      Gauss(2,8)=0.233193764592;
      Gauss(2,9)=0.186290210928;
      Gauss(2,10)=0.125580369465;
      Gauss(2,11)=0.0556685671162;
      
  case 12
      Gauss(1,1)=-0.981560634247;
      Gauss(1,2)=-0.90411725637;
      Gauss(1,3)=-0.769902674194;
      Gauss(1,4)=-0.587317954287;
      Gauss(1,5)=-0.367831498998;
      Gauss(1,6)=-0.125233408511;
      Gauss(1,7)=0.125233408511;
      Gauss(1,8)=0.367831498998;
      Gauss(1,9)=0.587317954287;
      Gauss(1,10)=0.769902674194;
      Gauss(1,11)=0.90411725637;
      Gauss(1,12)=0.981560634247;
      
      Gauss(2,1)=0.0471753363866;
      Gauss(2,2)=0.106939325995;
      Gauss(2,3)=0.160078328543;
      Gauss(2,4)=0.203167426723;
      Gauss(2,5)=0.233492536538;
      Gauss(2,6)=0.249147045813;
      Gauss(2,7)=0.249147045813;
      Gauss(2,8)=0.233492536538;
      Gauss(2,9)=0.203167426723;
      Gauss(2,10)=0.160078328543;
      Gauss(2,11)=0.106939325995;
      Gauss(2,12)=0.0471753363866;
      
  otherwise display('Error:No more Gauss Points..');
  end