Non-A Numéro Matlab

Question

clear 

% input problem mesh 
node = [ 0 0; 
     20 0; 
     20 15; 
      0 15]*12; 

conn = [1 2; 
     1 3; 
     2 4; 
     2 3; 
     4 3; 
     1 4]; 

% properties 
A = 1.0; 
E = 10e6; 

% boundary conditions 
P = 1000; 
f=zeros(8,1); 
f(6)=P; 

isol = [3 4 5 6 ]; %unconstrained dof 

%-----------------END OF INPUT---------------------- 

nn = size(node,1); %number of nodes 
ne = size(conn,1); %number of elements 
ndof = 2*nn; % number of dof 

K = zeros(ndof,ndof); %global stiffness matrix 
d = zeros(ndof, 1); %displacement vector 

for e=1:ne %loop through all the elements 

     n1=conn(e,1); 
     n2=conn(e,2); 

     x1=node(n1,1); y1=node(n1,2); % x and y coordinate for the first node 
     x2=node(n2,1); y2=node(n2,2); %x and y coordinate for the second node 

     L =sqrt( (x2-x1)^2 + (y2-y2)^2 ); %Length of the element 

     C = (x2-x1)/L; %Cosine 
     S = (y2-y1)/L; %Sine 
     C2 = C*C; % Cosine square 
     S2 = S*S; % Sine square 
     CS = C*S; 

     ke =(A*E/L)*[ C2 CS -C2 -CS; 
        CS S2 -CS -S2; 
        -C2 -CS C2 CS; 
        -CS -S2 CS S2 ]; %Stiffness matrix for element e 

     sctr = [2*n1-1 2*n1 2*n2-1 2*n2]; %Location where ke is to scatter 
              %in the global stiffness matrix 
              %scatter array 


     K(sctr,sctr) = K(sctr,sctr) + ke; 
end 

%solve for the nodal displacement, d 
d(isol) = K(isol,isol)\f(isol); 
fprintf('\n N o d a l D i s p l a c e m e n t\n') 
fprintf(' NID  X-Force Y-Force\n') 
for i=1:nn 
     fprintf('%5d %8.3e  %8.3e\n',i,d(2*i-1),d(2*i)) 
end 

%compute the reaction forces 
f = K*d; 

fprintf('\n E x t e r n a l F o r c e s\n') 
fprintf(' NID  X-FORCE Y-FORCE\n') 
for i=1:nn 
     fprintf('%5d %8.3e  %8.3e\n',i,f(2*i-1),f(2*i)) 
end 

%compute the element stress 
fprintf('\n E l e m e n t S t r e s s\n') 
fprintf('EID  STRAIN  STRESS\n') 
for e = 1:ne 

     n1 = conn(e,1); 
     n2 = conn(e,2); 

     x1 = node(n1,1); y1 = node(n1,2); 
     x2 = node(n2,1); y2 = node(n2,2); 

     L = sqrt((x2-x1)^2 + (y2-y1)^2); 
     C =(x2-x1)/L; 
     S =(y2-y1)/L; 

     B = 1/L*[-C -S C S]; %the element b matrix 

     sctr = [2*n1-1 2*n1 2*n2-1 2*n2 ]; 
     strain = B*d(sctr); 
     stress = E*strain; 

     fprintf('% 5d  %8.3e %8.3e\n',i,strain,stress) 

end 

Ici, la valeur de C est 0 et la valeur de S est 1 mais C2 = C * C, S2 = S * S, CS = C * S sont donnant NaN en sortie. Comment puis-je empêcher la valeur NaN d'obtenir les valeurs réelles de S2, C2, CS obtenues par l'opération arithmétique? C'est un code d'analyse par éléments finis fait dans MATLAB.

Answer 1

semble que ce soit une erreur (comme @Adiel dit dans les commentaires)

L =sqrt( (x2-x1)^2 + (y2-y2)^2 ); %Length of the element

à y2-y2, qui devrait être y2-y1.