Not-A-Number 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   

Here the value of C is 0 and value of S is 1 but C2 = C * C, S2 = S * S, CS = C * S are giving NaN as output. How can I prevent NaN value to get the actual values of S2, C2, CS obtained by the arithmetic operation? This is a finite element analysis code done in MATLAB.

Answer 1

Seems this is an error (as @Adiel said in comments)

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

at y2-y2, which should be y2-y1.