% MEAM 599- 02 / 1998 Winter
% FEM for Side Frame Analysis using plane beam elements
% with possibly flexible joints
%                                                     N. Kikuchi
% ______________________________________________________________
%
% Pre-Processing / set up an analysis model
%
nfile=input('Is data in the datafile ?  [y/n] = ','s');
if nfile~='y' 
% read nodal coordinates ( Z, Y )
nx=input('number of nodes of the whole structure = ');
for i=1:nx
	i
	Z(i)=input('Z coordinate of node = ');
	Y(i)=input('Y coordinate of node = ');
end
% plot the nodes of the plane beam structure
plot(Z,Y,'+')
% read element connectivity and section type of beam elements
%     ijk(1,nel)=node i of beam element nel
%     ijk(2,nel)=node j of beam element nel
%     ijk(3,nel)=section type of beam element nel
nelx=input('number of beam elements = ');
for nel=1:nelx
	nel
	ijk(1,nel)=input('node i = ');
	ijk(2,nel)=input('node j = ');
	ijk(3,nel)=input('section type = ');
end
% plot the beam elements
for nel=1:nelx
	Ze(2*nel-1)=Z(ijk(1,nel));
	Ze(2*nel)=Z(ijk(2,nel));
	Ye(2*nel-1)=Y(ijk(1,nel));
	Ye(2*nel)=Y(ijk(2,nel));
end
plot(Z,Y,'+',Ze,Ye)
% read section type ( properties )
nsecx=input('total number of section type = ');
for i=1:nsecx
	i
	E(i)=input('Young,s modulus / E = ');
	A(i)=input('crossectional area / A = ');
	Iyy(i)=input('moment of inertia about the x axis / Iyy = ');
	Fy(i)=input('shear constant / Fy = ');
end
% read data for flexible joints
nfjx=input('number of flexible joints = ');
if nfjx>0
	for i=1:nfjx
		i
		fjoint(1,i)=input('node i of the flexible joint = ');
		fjoint(2,i)=input('node j of the flexible joint = ');
		kZ(i)=input('stiffness percent in the Z direction / kZ = ');
		kY(i)=input('stiffness percent in the Y direction / kY = ');
		kqX(i)=input('stiffness percent about the X axis rotation / kqX = ');
	end
end
% plot the flexible joints
if nfjx>0
	for i=1:nfjx
		Zfj(i)=Z(fjoint(1,i));
		Yfj(i)=Y(fjoint(1,i));
	end
plot(Z,Y,Ze,Ye,'+',Zfj,Yfj,'o')
end
% displacement constraints
spc=[];
nspc=input('number of single point constraints = ');
for i=1:nspc
	i
	spc(1,i)=input('node number = ');
	spc(2,i)=input('degree of freedom for spc = ');
	spc(3,i)=input('constrained value = ');
end	
% applied forces and moments at the nodes
afm=[];
nafm=input('number of applied forces and moment = ');
for i=1:nafm
	i
	afm(1,i)=input('node number = ');
	afm(2,i)=input('degrees of freedom = ');
	afm(3,i)=input('applied forces or moment = ');
end
%
save input.mat nx Y Z nelx ijk nsecx E A Iyy Fy nfjx fjoint kY kZ kqX nspc spc nafm afm
end
if nfile=='y'
	load input.mat
%	load input.mat nx Y Z nelx ijk nsecx E A Iyy Fy nfjx fjoint kY kZ kqX nspc spc nafm afm
end
%
% FE-Processing / forming the global stiffness matrix
%
% beam elements
%
sk=zeros(3*nx);
f=zeros(3*nx,1);
sksize=3*nx;
%
for nel=1:nelx
%
	Zji=Z(ijk(2,nel))-Z(ijk(1,nel));
	Yji=Y(ijk(2,nel))-Y(ijk(1,nel));
	Lji=sqrt(Zji^2+Yji^2);
	cji=Zji/Lji;
	sji=Yji/Lji;
	TG=zeros(6);
	TG(1,1)= cji;
	TG(1,2)=-sji;
	TG(2,1)= sji;
	TG(2,2)= cji;
	TG(3,3)= 1;
	TG(4,4)= cji;
	TG(4,5)=-sji;
	TG(5,4)= sji;
	TG(5,5)= cji;
	TG(6,6)= 1;
%
	skel=zeros(6);
	ijk3=ijk(3,nel);
	EIyy=E(ijk3)*Iyy(ijk3);
	Fy1=(1+Fy(ijk3));
	skel(1,1)=12*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji^3);
	skel(2,2)=E(ijk3)*A(ijk3)/Lji;
	skel(3,1)=-6*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji^2);
	skel(3,3)=(4+Fy(ijk3))*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji);
	skel(4,1)=-skel(1,1);
	skel(4,3)=-skel(3,1);
	skel(4,4)= skel(1,1);
	skel(5,2)=-skel(2,2);
	skel(5,5)= skel(2,2);
	skel(6,1)= skel(3,1);
	skel(6,3)=(2-Fy(ijk3))*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji);
	skel(6,4)=-skel(3,1);
	skel(6,6)= skel(3,3);
	for i=1:5
		for j=1+1:6
			skel(i,j)=skel(j,i);
		end
	end
%
	skeg=TG'*skel*TG;
%
	ijk1=ijk(1,nel);
	ijk2=ijk(2,nel);
	ndg=[3*ijk1-2,3*ijk1-1,3*ijk1,3*ijk2-2,3*ijk2-1,3*ijk2];
	for i=1:6
		for j=1:6
			sk(ndg(i),ndg(j))=sk(ndg(i),ndg(j))+skeg(i,j);
		end
	end
%
end
%
% flexible joint elements
%
if nfjx>0
	for nfj=1:nfjx
		skeg=zeros(6);
		ijk1=fjoint(1,nfj);
		ijk2=fjoint(2,nfj);
		kYb=(sk(3*ijk1-2,3*ijk1-2)+sk(3*ijk2-2,3*ijk2-2))/2;
		kZb=(sk(3*ijk1-1,3*ijk1-1)+sk(3*ijk2-1,3*ijk2-1))/2;
		kqXb=(sk(3*ijk1,3*ijk1)+sk(3*ijk2,3*ijk2))/2;
		refrencestiffness=[nfj,kYb,kZb,kqXb]
		skeg(1,1)= kY(nfj)*kYb;
		skeg(4,4)= kY(nfj)*kYb;
		skeg(1,4)=-kY(nfj)*kYb;
		skeg(4,1)=-kY(nfj)*kYb;
		skeg(2,2)= kZ(nfj)*kZb;
		skeg(5,5)= kZ(nfj)*kZb;
		skeg(2,5)=-kZ(nfj)*kZb;
		skeg(5,2)=-kZ(nfj)*kZb;
		skeg(3,3)= kqX(nfj)*kqXb;
		skeg(6,6)= kqX(nfj)*kqXb;
		skeg(3,6)=-kqX(nfj)*kqXb;
		skeg(6,3)=-kqX(nfj)*kqXb;
		ndg=[3*ijk1-2,3*ijk1-1,3*ijk1,3*ijk2-2,3*ijk2-1,3*ijk2];
		for i=1:6
			for j=1:6
				sk(ndg(i),ndg(j))=sk(ndg(i),ndg(j))+skeg(i,j);
			end
		end
	end
end
%
% single point constraints
%
for i=1:nspc
	spc1=spc(1,i);
	spc2=spc(2,i);
	spc3=spc(3,i);
	ndof=3*(spc1-1)+spc2;
	f=f-sk(:,ndof)*spc3;
	sk(:,ndof)=zeros(sksize,1);
	sk(ndof,:)=zeros(1,sksize);
	sk(ndof,ndof)=1;
end
%
% applied for#es and moments
%	
for i=1:nafm
	afm1=afm(1,i);
	afm2=afm(2,i);
	afm3=afm(3,i);
	ndof=3*(afm1-1)+afm2;
	f(ndof)=f(ndof)+afm3;
end
%
% solving the matrix equation
%
d=sk\f;
d'
%
% Post-Processing of the computed results
%
strainenergy=[];
Zp=[];
Yp=[];
uZ=[];
uY=[];
for nel=1:nelx
	Zji=Z(ijk(2,nel))-Z(ijk(1,nel));
	Yji=Y(ijk(2,nel))-Y(ijk(1,nel));
	Lji=sqrt(Zji^2+Yji^2);
	cji=Zji/Lji;
	sji=Yji/Lji;
	TG=zeros(6);
	TG(1,1)= cji;
	TG(1,2)=-sji;
	TG(2,1)= sji;
	TG(2,2)= cji;
	TG(3,3)= 1;
	TG(4,4)= cji;
	TG(4,5)=-sji;
	TG(5,4)= sji;
	TG(5,5)= cji;
	TG(6,6)= 1;
	deg=zeros(6,1);
	for i=1:2
		ijki=ijk(i,nel);
		for j=1:3
			deg(3*(i-1)+j)=d(3*(ijki-1)+j);
		end
	end
	del=TG*deg;
	skel=zeros(6);
	ijk3=ijk(3,nel);
	EIyy=E(ijk3)*Iyy(ijk3);
	Fy1=(1+Fy(ijk3));
	skel(1,1)=12*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji^3);
	skel(2,2)=E(ijk3)*A(ijk3)/Lji;
	skel(3,1)=-6*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji^2);
	skel(3,3)=(4+Fy(ijk3))*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji);
	skel(4,1)=-skel(1,1);
	skel(4,3)=-skel(3,1);
	skel(4,4)= skel(1,1);
	skel(5,2)=-skel(2,2);
	skel(5,5)= skel(2,2);
	skel(6,1)= skel(3,1);
	skel(6,3)=(2-Fy(ijk3))*E(ijk3)*Iyy(ijk3)/((1+Fy(ijk3))*Lji);
	skel(6,4)=-skel(3,1);
	skel(6,6)= skel(3,3);
	for i=1:5
		for j=1+1:6
			skel(i,j)=skel(j,i);
		end
	end
	nel
	strainenergynel=(del'*skel*del)/2
	strainenergy(nel)=strainenergynel;
% axial strain  ez = ez0 + y ( ezi*(1-z/l) + ezj*(z/l) )
	ez0=(del(5)-del(2))/Lji;
	ezi= 6*(del(1)-del(4))/Lji^2-(4*del(3)+2*del(6))/Lji;
	ezj=-6*(del(1)-del(4))/Lji^2+(2*del(3)+4*del(6))/Lji;
	axialstrain=[ez0,ezi,ezj]
% axial stress  sz = sz0 + y ( szi*(1-z/l) + szj*(z/l) )
	sz0=E(ijk3)*ez0;
	szi=E(ijk3)*ezi;
	szj=E(ijk3)*ezj;
	axialstress=[sz0,szi,szj]
% displacement of the beam axis
	ipx=11;
	for ip=1:ipx
		zpi=(ip-1)/(ipx-1);
		nv=[1-3*zpi^2+2*zpi^3,Lji*(zpi-2*zpi^2+zpi^3),3*zpi^2-2*zpi^3,Lji*(-zpi^2+zpi^3)];
		nw=[1-zpi,zpi];
		uyi=nv*[del(1),-del(3),del(4),-del(6)]';
		uzi=nw*[del(2),del(5)]';
		ipnel=ipx*(nel-1)+ip;
		Yp(ipnel)=nw*[Y(ijk(1,nel)),Y(ijk(2,nel))]';
		Zp(ipnel)=nw*[Z(ijk(1,nel)),Z(ijk(2,nel))]';
		uY(ipnel)= uyi*cji+uzi*sji;
		uZ(ipnel)=-uyi*sji+uzi*cji;
	end
end
if nfjx>0
	for nfj=1:nfjx
		skeg=zeros(6);
		ijk1=fjoint(1,nfj);
		ijk2=fjoint(2,nfj);
		kYb=(sk(3*ijk1-2,3*ijk1-2)+sk(3*ijk2-2,3*ijk2-2))/2;
		kZb=(sk(3*ijk1-1,3*ijk1-1)+sk(3*ijk2-1,3*ijk2-1))/2;
		kqXb=(sk(3*ijk1,3*ijk1)+sk(3*ijk2,3*ijk2))/2;
		skeg(1,1)= kY(nfj)*kYb;
		skeg(4,4)= kY(nfj)*kYb;
		skeg(1,4)=-kY(nfj)*kYb;
		skeg(4,1)=-kY(nfj)*kYb;
		skeg(2,2)= kZ(nfj)*kZb;
		skeg(5,5)= kZ(nfj)*kZb;
		skeg(2,5)=-kZ(nfj)*kZb;
		skeg(5,2)=-kZ(nfj)*kZb;
		skeg(3,3)= kqX(nfj)*kqXb;
		skeg(6,6)= kqX(nfj)*kqXb;
		skeg(3,6)=-kqX(nfj)*kqXb;
		skeg(6,3)=-kqX(nfj)*kqXb;
		for i=1:2
			ijki=fjoint(i,nfj);
			for j=1:3
				deg(3*(i-1)+j)=d(3*(ijki-1)+j);
			end
		end
		nfj
		strainenergynfj=(deg'*skeg*deg)/2
		strainenergy(nelx+nfj)=strainenergynfj;
	end
end
bar(strainenergy)
xlabel('beam elements and flexible joints')
ylabel('strain energy')
title('strain energy distribution')
% plot the deformed configuration
mag=input('magnification factor of the fisplacement = ');
Zup=Zp+mag*uZ;
Yup=Yp+mag*uY;
plot(Zp,Yp,Zup,Yup)
title('deformed configuration')
xlabel('Z')
ylabel('Y')

