Can someone rearrange the code to run

Question

MINATI 2020-6-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/554314-can-someone-rearrange-the-code-to-run

评论： MINATI 2020-6-25

采纳的回答： Rafael Hernandez-Walls

%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

% 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

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

2 个评论
显示无隐藏无

Rafael Hernandez-Walls 2020-6-25

You just have to define the function mass

MINATI 2020-6-25

Dear Rafael

Like stiff, mass is also defined in line 31 (I guess) or, please modify if possible

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Rafael Hernandez-Walls 2020-6-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/554314-can-someone-rearrange-the-code-to-run#answer_456745

在 MATLAB Online 中打开

%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