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