# Can someone rearrange the code to run

3 views (last 30 days)
MINATI on 25 Jun 2020
Commented: MINATI on 25 Jun 2020
%subdivisions space /time
ht=0.01; Tmax=1.2; nx=33; hx=1/(nx-1); x=[0:hx:1]';
%matrices
K=stiff(1/pi^2,hx,nx); M=mass(1/ht,hx,nx); A=M+K;
%boundary conditions by deleting rows
A(nx,:)=[]; A(:,nx)=[]; A(1,:)=[]; A(:,1)=[];
%creation
R=chol(A);
%initial condition
u=cos(pi*(x-1/2));
%time step loop
k=0;
while (k*ht<Tmax)
k=k+1;
b=M*u;
b(nx)=[];b(1)=[];
%solve
u=zeros(nx,1);
u(2:nx-1)=R\(R'\b);
end
function R=stiff(nu,h,n)
%
[m1,m2]=size(nu);
if (length(nu)==1)
ee=nu*ones(n-1,1)/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=.5*(nu(1:n-1)+nu(2:n))/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
R=spdiags([-e1 e -e2],-1:1,n,n);
return
function M=mass(alpha,h,n)
[m1,m2]=size(alpha);
if(length(alpha)==1)
ee=alpha*h*ones(n-1,1)/6;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=h*(alpha(1:n-1)+alpha(2:n))/12;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
M=spdiags([e1 2*e e2],-1:1,n,n);
return
% Set the matrices to zero
k=zeros(6,6,2); K=zeros(6,6,2); Gamma=zeros(6,6,2);
% Enter parameter values, in units: mm^2, mm^4, and MPa(10^6 N/m)
A=6000; II=200*10^6; EE=200000;
% Convert units into meters and kN
A=A/10^6; II = II/10^12; EE =EE*1000;
% Element 1
i=[0,0]; j=[7.416,3];
[k(:,:,1),K(:,:,1),Gamma(:,:,1)]=stiff(EE,II,A,i,j);
% Element 2
i=j; j=[15.416,3];
[k(:,:,2),K(:,:,2),Gamma(:,:,2)]=stiff(EE,II,A,i,j);
% Define the elementary stiffness matrix
ID=[-4 1 -7;-5 2 -8;-6 3 -9];
% Define the connection matrix
LM=[-4 -5 -6 1 2 3; 1 2 3 -7 -8 -9];
% Assemble the augmented stiffness matrix
Kaug = zeros(9);
for elem=1:2
for r=1:6
lr=abs(LM(elem,r));
for c=1:6
lc=abs(LM(elem,c));
Kaug(lr,lc)=Kaug(lr,lc)+K(r,c,elem);
end
end
end
% Extract the stiffness matrix
Ktt=Kaug(1:3,1:3);
% Determine the reactions at the nodes in local coordinates
fea(1:6,1)=0;
fea(1:6,2)=[0,8*4/2,4*8^2/12,0,8*4/2,-4*8^2/12]';
% Determine the reactions in the global coordinate system
FEA(1:6,1)=Gamma(:,:,1)'*fea(1:6,1);
FEA(1:6,2)=Gamma(:,:,2)'*fea(1:6,2);
% FEA_Rest for constrained nodes
FEA_Rest=[0,0,0,FEA(4:6,2)'];
% Assemble the right-hand side for non-constrained nodes
P(1)=50*3/8; P(2)=-50*7.416/8-FEA(2,2); P(3)=-FEA(3,2);
% Solve to find the displacements in meters and in radians
Displacements=inv(Ktt)*P'
% Extract Kut
Kut = Kaug(4:9,1:3);
% Compute the reactions and introduce boundary conditions
Reactions=Kut*Displacements+FEA_Rest'
% Solve to find the internal forces, excluding fixed points
dis_global(:,:,1)=[0,0,0,Displacements(1:3)'];
dis_global(:,:,2)=[Displacements(1:3)',0,0,0];
for elem=1:2
dis_local= Gamma(:,:,elem)*dis_global(:,:,elem)';
int_forces= k(:,:,elem)*dis_local+fea(1:6,elem)
end
%The above script calls the stiff function, which can be implemented as follows:
function [k,K,Gamma] = stiff( EE,II,A,i,j )
% Find the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);
% Compute the angle theta
if(j(1)-i(1))~=2
alpha=atan((j(2)-i(2))/(j(1)-i(1)))
else
alpha=-pi/2;
end
% Form the rotation matrix Gamma
Gamma =[cos(alpha) sin(alpha) 0 0 0 0;
-sin(alpha) cos(alpha) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(alpha) sin(alpha) 0;
0 0 0 -sin(alpha) cos(alpha) 0;
0 0 0 0 0 1];
% Form the elementary stiffness matrix in local coordinates
EI=EE*II; EA=EE*A;
k=[EA/L, 0, 0, -EA/L, 0, 0;
0, 12*EI/L^3, 6*EI/L^2, 0, -12*EI/L^3,6*EI/L^2;
0, 6*EI/L^2, 4*EI/L, 0 -6*EI/L^2, 2*EI/L;
-EA/L, 0 ,0 , EA/L, 0, 0;
0, -12*EI/L^3, -6*EI/L^2, 0, 12*EI/L^3, -6*EI/L^2;
0, 6*EI/L^2, 2*EI/L, 0, -6*EI/L^2, 4*EI/L];
% Elementary matrix in global coordinates
K=Gamma'*k*Gamma;
end
end
end
%%%% I got this code from a book to solve heat eqn but unable to arrange this.
%% Please someone arrange this so that it could be of use

Rafael Hernandez-Walls on 25 Jun 2020
You just have to define the function mass
MINATI on 25 Jun 2020
Dear Rafael
Like stiff, mass is also defined in line 31 (I guess) or, please modify if possible

Rafael Hernandez-Walls on 25 Jun 2020
%subdivisions space /time
ht=0.01; Tmax=1.2; nx=33; hx=1/(nx-1); x=[0:hx:1]';
%matrices
K=stiff2(1/pi^2,hx,nx); M=mass(1/ht,hx,nx);
A=M+K;
%boundary conditions by deleting rows
A(nx,:)=[]; A(:,nx)=[]; A(1,:)=[]; A(:,1)=[];
%creation
R=chol(A);
%initial condition
u=cos(pi*(x-1/2));
%time step loop
k=0;
while (k*ht<Tmax)
k=k+1;
b=M*u;
b(nx)=[];b(1)=[];
%solve
u=zeros(nx,1);
u(2:nx-1)=R\(R'\b);
end
% Set the matrices to zero
k=zeros(6,6,2); K=zeros(6,6,2); Gamma=zeros(6,6,2);
% Enter parameter values, in units: mm^2, mm^4, and MPa(10^6 N/m)
A=6000; II=200*10^6; EE=200000;
% Convert units into meters and kN
A=A/10^6; II = II/10^12; EE =EE*1000;
% Element 1
i=[0,0]; j=[7.416,3];
[k(:,:,1),K(:,:,1),Gamma(:,:,1)]=stiff(EE,II,A,i,j);
% Element 2
i=j; j=[15.416,3];
[k(:,:,2),K(:,:,2),Gamma(:,:,2)]=stiff(EE,II,A,i,j);
% Define the elementary stiffness matrix
ID=[-4 1 -7;-5 2 -8;-6 3 -9];
% Define the connection matrix
LM=[-4 -5 -6 1 2 3; 1 2 3 -7 -8 -9];
% Assemble the augmented stiffness matrix
Kaug = zeros(9);
for elem=1:2
for r=1:6
lr=abs(LM(elem,r));
for c=1:6
lc=abs(LM(elem,c));
Kaug(lr,lc)=Kaug(lr,lc)+K(r,c,elem);
end
end
end
% Extract the stiffness matrix
Ktt=Kaug(1:3,1:3);
% Determine the reactions at the nodes in local coordinates
fea(1:6,1)=0;
fea(1:6,2)=[0,8*4/2,4*8^2/12,0,8*4/2,-4*8^2/12]';
% Determine the reactions in the global coordinate system
FEA(1:6,1)=Gamma(:,:,1)'*fea(1:6,1);
FEA(1:6,2)=Gamma(:,:,2)'*fea(1:6,2);
% FEA_Rest for constrained nodes
FEA_Rest=[0,0,0,FEA(4:6,2)'];
% Assemble the right-hand side for non-constrained nodes
P(1)=50*3/8; P(2)=-50*7.416/8-FEA(2,2); P(3)=-FEA(3,2);
% Solve to find the displacements in meters and in radians
Displacements=inv(Ktt)*P'
% Extract Kut
Kut = Kaug(4:9,1:3);
% Compute the reactions and introduce boundary conditions
Reactions=Kut*Displacements+FEA_Rest'
% Solve to find the internal forces, excluding fixed points
dis_global(:,:,1)=[0,0,0,Displacements(1:3)'];
dis_global(:,:,2)=[Displacements(1:3)',0,0,0];
for elem=1:2
dis_local= Gamma(:,:,elem)*dis_global(:,:,elem)';
int_forces= k(:,:,elem)*dis_local+fea(1:6,elem)
end
%The above script calls the stiff function, which can be implemented as follows:
function [k,K,Gamma] = stiff( EE,II,A,i,j )
% Find the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);
% Compute the angle theta
if(j(1)-i(1))~=2
alpha=atan((j(2)-i(2))/(j(1)-i(1)))
else
alpha=-pi/2;
end
% Form the rotation matrix Gamma
Gamma =[cos(alpha) sin(alpha) 0 0 0 0;
-sin(alpha) cos(alpha) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(alpha) sin(alpha) 0;
0 0 0 -sin(alpha) cos(alpha) 0;
0 0 0 0 0 1];
% Form the elementary stiffness matrix in local coordinates
EI=EE*II; EA=EE*A;
k=[EA/L, 0, 0, -EA/L, 0, 0;
0, 12*EI/L^3, 6*EI/L^2, 0, -12*EI/L^3,6*EI/L^2;
0, 6*EI/L^2, 4*EI/L, 0 -6*EI/L^2, 2*EI/L;
-EA/L, 0 ,0 , EA/L, 0, 0;
0, -12*EI/L^3, -6*EI/L^2, 0, 12*EI/L^3, -6*EI/L^2;
0, 6*EI/L^2, 2*EI/L, 0, -6*EI/L^2, 4*EI/L];
% Elementary matrix in global coordinates
K=Gamma'*k*Gamma;
end
function R=stiff2(nu,h,n)
%
[m1,m2]=size(nu);
if (length(nu)==1)
ee=nu*ones(n-1,1)/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=.5*(nu(1:n-1)+nu(2:n))/h; e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
R=spdiags([-e1 e -e2],-1:1,n,n);
return
end
function M=mass(alpha,h,n)
[m1,m2]=size(alpha);
if(length(alpha)==1)
ee=alpha*h*ones(n-1,1)/6;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
else
ee=h*(alpha(1:n-1)+alpha(2:n))/12;
e1=[ee;0]; e2=[0;ee]; e=e1+e2;
end
M=spdiags([e1 2*e e2],-1:1,n,n);
return
end

#### 1 Comment

MINATI on 25 Jun 2020
Dear Rafael
it worked
Thank you