function Plotdisp(Order)
%1D Postprocess routine for OneD FEM CODE, will read output files from 1DFEM
%AUTHOR: S. Srivastava, University of Twente, NL
%Run as Plot1D(1), Plot1D(2) ..etc.
global L Area Youngs Traction Total_Elements TSF_U TSF_Sigma Coord_Table L2GDofTable ;

for i=1:Order+1
    N(:,i)=putzero(N1D(i),Order+1);
end
N

display('PostProcessing Solution for Displacement in 1D (Input files: 1d.dof 1d.sol).......');

File=strcat('1D','.sol');
sol1=dlmread(File,',');
File=strcat('1D','.dof');
L2GDofTable=dlmread(File,',');

figure(1);
clf;

subplot(211); %plot two curves in one window
Dplot(sol1,N,Order);  %plot displacement

title(strcat('Displacement Ux ',''));
xlabel('X Axis(Along Bar Length) (m)');
ylabel('Ux (m)');

subplot(212);
Splot(sol1,N,Order); %plot stress

title(strcat('Stress Sigmax ',''));
xlabel('X Axis(Along Bar Length) (m)');
ylabel('Sigma (N/m^2)');

display ('Disp Plot Over....................');


function Dplot(Coeff,N,order)
global L TSF Coord_Table Total_Elements;

i=0;
j=0;
num=0;
coeff=0;
X=-1:0.1:1;
i=length(X);
Z=0;
J=0;
counter=1;

TSF=order+1;
i=length(X);
for num=1:Total_Elements
    num
    coeff=LocalCoeff(num,Coeff);
    for il=1:i
        Z(counter)=FindDisp(coeff,num,N,X(il));
        counter=counter+1;  
    end
    
end

grid on;
colormap('default');
X=0:(1/(length(Z)-1)):L;
plot(X,Z);


function Splot(Coeff,N,order)
global L TSF Coord_Table Total_Elements;
i=0;
j=0;
num=0;
coeff=0;
X=-1:0.1:1;
i=length(X);
Z=0;
J=0;
counter=1;
%Generate Shape Funcitons
TSF=order+1;
i=length(X);
for num=1:Total_Elements
    num
    coeff=LocalCoeff(num,Coeff);
    for il=1:i
        Z(counter)=FindStress(coeff,num,N,X(il));
        counter=counter+1;  
    end
    
end

grid on;
colormap('default');
X=0:(1/(length(Z)-1)):L;
plot(X,Z);

function [U]=FindDisp(coeff,E,ShapeFunc,x)
%evaluates stress at a point x,y,z of type material and return the stress matrix
U=0;
[c,TSF]=size(ShapeFunc);

%Evaluate Ux Uy Uz using FEM basic interpolation
if length(coeff)~=TSF
    error('Error in FindDisp...check  here');
end
for i=1:TSF
    U=U+coeff(i)*Evaluate(ShapeFunc(:,i),x);
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 [S]=FindStress(coeff,E,ShapeFunc,x)
global Youngs;
%evaluates stress at a point x,y,z of type material and return the stress matrix
S=0;
[c,TSF]=size(ShapeFunc);

%Evaluate Ux Uy Uz using FEM basic interpolation
if length(coeff)~=TSF
    error('Error in FindDisp...check  here');
end
for i=1:TSF
    S=S+coeff(i)*Evaluate(Derivative(ShapeFunc(:,i)),x);
end
S=Youngs*S;

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 [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 [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 [answer]=LocalCoeff(E,Coeff)
global TSF L2GDofTable;

for i=1:TSF
    answer(i)=Coeff(L2GDofTable(E,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