Why am I getting a error 'Dimensions of arrays being concatenated are not consistent'

Why am i getting an error message. I thought i was producing a 12x12 zeros matrix, then adding the matrix calculcated to the zeroes. This works when I use a 4x4. why doesnt 12x12 work. LINE 117 is the error.

%
% BarFE.m
%
% The purpose of BarFE.m is to introduce you to coding your own
% Finite element software. Herein, we will use bar elements in 2D for
% simplicity.
% 
% Use this file as a template to expand it to 3D, then further add
% DoF like shaft and beams.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Written by:     Terence Macquart.
%   Last Updated:   07/04/2021.
%   Updated by:     
%   Change log:
%
%
%   Contact details:
%   BCI (ACCIS) - Bristol Composite Institute
%   Department of Aerospace Engineering - University of Bristol
%   Queen's Building - University Walk
%   Bristol - BS8 1TR
%   U.K.
%   E-mail: [e-mail]
%   Phone:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% swap folder and add any subfolders/files to your matlab path
if(~isdeployed),  cd(fileparts(which('BarFE_Example')));end

Error using cd
Unable to change current folder to '' (Name is nonexistent or not a folder).

clear; close all; clc
maindir = pwd;
maindir(maindir=='\')='/'; % for linux
addpath (genpath(strcat(maindir)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Step 1/2 geometry 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% node and element connectivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% node coordinates (set 3rd column to 0 for 2D)
nodes = [...
    0, 0, 0 ;
    1, 1, 0;
    0, 1, 0;
    1, 0, 0;
    2, 1, 0;
    2, 0, 0;
    3, 1, 0];

% element connectivity matrix
eConn = [...
    1, 2;
    3, 2;
    1, 4;
    4, 5;
    4, 6;
    2, 5;
    6, 7;
    5, 7;
    4, 2;
    6, 5];

ne = size(eConn,1); % number of bar elements

% plot the structure nodes
figure('name','Nodes and Elements')
hold on
plot(nodes(:,1), nodes(:,2), 'blue o', 'markerfacecolor', 'blue')
xlabel('x (m)')
ylabel('y (m)')
axis equal
grid on
for i=1:size(nodes,1)
    text(nodes(i,1)-0.05,nodes(i,2)-0.05,num2str(i),'color','red')
end


% compute theta and length for each element
Le    = zeros(ne,1);    % element length
theta = zeros(ne,1);    % element orientation(rad)
for i=1:ne
    elemtNode1 = nodes(eConn(i,1),:)';
    elemtNode2 = nodes(eConn(i,2),:)';
    opp      = elemtNode2(2) - elemtNode1(2);
    adj      = elemtNode2(1) - elemtNode1(1);
    theta(i) = atan2(opp,adj);
    Le(i)    = sqrt( sum((elemtNode2-elemtNode1).^2) ); % = norm(elemtNode2-elemtNode1)
    
    % plot the structure elements
    plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'black')
    text(mean([elemtNode1(1), elemtNode2(1)]),mean([elemtNode1(2), elemtNode2(2)]),num2str(i))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Step 3/4/5 Element stiffness matrices (bars only in this example)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% element properties and stiffness matrix
%%%%%%% all element are assumed to have the same area and Young's modulus
%%%%%%% Ke is the 2D element stiffness matrix in the element ref. frame
%%%%%%% Kg is the 2D element stiffness matrix in the global  ref. frame
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 0.01;     % bar cross-sectional area (m^2)
E = 190*1e9;  % Young's modulus (Pa)
Izz = 5000;
Iyy = 5000;
G = 5000;
J = 5000;

Ke = zeros(12,12,ne);
Kg = zeros(12,12,ne);
for i=1:ne
    Ke(:,:,i) = [A*E/(Le(i)), 0, 0, 0, 0, 0, -A*E/(Le(i)), 0, 0, 0, 0, 0; 
        0, 12*E*Izz/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Izz/Le(i)^3, 0, -12*E*Izz/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Izz/Le(i)^3; 
        0, 0, 12*E*Iyy/Le(i)^3, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 0, 0, -12*E*Iyy/Le(i)^3, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0; 
        0, 0, 0, G*J/Le(i), 0, 0, 0, 0, 0, -G*J/Le(i), 0, 0; 
        0, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 4*Le(i)^2*E*Iyy/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 2*Le(i)^2*E*Iyy/Le(i)^3, 0; 
        6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3; 
        A*E/Le(i), 0, 0, 0, 0, 0, A*E/Le(i), 0, 0, 0, 0, 0; 
        0, -12*E*Izz/Le(i)^3, 0 ,0, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 12*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3; 
        0, 0, -12*E*Iyy/Le(i)^3, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 0, 0, 12*E*Iyy/Le(i)^3, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0; 
        0, 0, 0, -G*J/Le(i), 0, 0, 0, 0, 0, G*J/Le(i), 0, 0; 
        0, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 2*Le(i)^2*E*Iyy/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 4*Le(i)^2*E*Iyy/Le(i)^3, 0; 
        0, 6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3];
    
    R = [cos(theta(i)) -sin(theta(i));
        sin(theta(i)) cos(theta(i))];
    R = blkdiag(R,R);
    
    Kg(:,:,i) = R*Ke(:,:,i)*R';
end



%% Step 6 Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% element assembly
%%%%%%% the global stiffness matrix is 14x14 since we have 2 degrees of
%%%%%%% freedom per node (translations in x and y), and 7 nodes
%%%%%%%
%%%%%%% here I chose to number the degrees of freedom as follows
%%%%%%% Dof = description
%%%%%%% 1   = node 1 translation in x
%%%%%%% 2   = node 1 translation in y
%%%%%%% 3   = node 2 translation in x
%%%%%%% 4   = node 3 translation in y
%%%%%%% ...
%%%%%%% 14  = node 7 translation in y
%%%%%%%
%%%%%%% this choice is arbitrary and you could chose another one.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K = zeros(14,14); % global stiffness matrix
for i=1:ne
    nodeID     = eConn(i,:);                  % index of nodes linked to element i
    id         = sort([nodeID*2-1 nodeID*2]); % index of dofs  linked to element i
    K(id,id)   = K(id,id) + Kg(:,:,i);        % add stiffnes of element i to global stiffness matrix
end


%% Step 7/8 Boundary conditions & Solve 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Apply Boundary conditions & loads, then solve
%%%%%%% Clamped at node 1 and 3
%%%%%%% vertical downward force at node 7 Fy = -1e4 Newtons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u       = zeros(14,1); 
F       = zeros(14,1);
F(end)  = -4e6;                  % external force (N) applied at the last dof


% remove the boundary conditions and solve the remaing system of equations
idBC        = [1 2 5 6];         % dofs set to zero due to clamps
idFree      = [3 4 7:14];        % remaining free dofs 
Kf          = K(idFree,idFree);
Ff          = F(idFree);
u(idFree)   = Kf^-1*Ff;          % displacement solution caused by F



%% Step 9 Post-processing

% format the output and visualise the results
u_xy            = reshape(u,2,[])';
DeformedNode    = nodes(:,1:2)+u_xy;
 

L = zeros(ne,1); % deformed element length
for i=1:ne
   elemtNode1 = DeformedNode(eConn(i,1),:)';
   elemtNode2 = DeformedNode(eConn(i,2),:)';
   L(i)       = sqrt( sum((elemtNode2-elemtNode1).^2) );
   
   if L(i)>Le(i)
      % member in tension
      plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'blue','linewidth',2)
   else
      % member in compression
      plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'red','linewidth',2)
   end
end



% strains due to axial deformation (from element summary table for bar)
eps_xx = zeros(ne,1);
for i=1:ne 
   % for each element find the nodes and DoFs indices 
   idNodes    = eConn(i,:)';
   Dof_Node1  = [idNodes(1)*2-1, idNodes(1)*2];
   Dof_Node2  = [idNodes(2)*2-1, idNodes(2)*2];
   
   
   ug = u([Dof_Node1,Dof_Node2]);    % global displacement of the nodes
   
   
   % convert ug into displacement of the nodes expressed in the element coordinate system
   R  = [cos(theta(i)) -sin(theta(i));
         sin(theta(i)) cos(theta(i))];
   R  = blkdiag(R,R);
   ue = R'*ug; 
   
   % evaluate the element axial strain 
   % you can check that eps_xx = (L(i)-Le(i))/Le(i) for small displacements
   eps_xx(i) = [-1/Le(i) 1/Le(i)]*ue([1 3]); 
end




%% Step 10 Validate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Validation against Abaqus results
%       Node Label            U.U1            U.U2
%                           @Loc 1          @Loc 1
% -------------------------------------------------
%                1    -6.31579E-03    -37.1723E-03
%                2    -4.21053E-03    -14.3756E-03
%                3     6.31579E-03    -12.2704E-03
%                4     10.5263E-03    -35.0671E-03
%                5    -12.0000E-30    -4.00000E-30
%                6     12.6316E-03    -62.0743E-03
%                7     12.0000E-30              0.
%
% abaqus node numbering is different than our so we need to make sure to
% compare the correct nodes together --> this is why we use idmap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

AbaqusData = [ 1    -6.31579E-03    -37.1723E-03
               2    -4.21053E-03    -14.3756E-03
               3     6.31579E-03    -12.2704E-03
               4     10.5263E-03    -35.0671E-03
               5    -12.0000E-30    -4.00000E-30
               6     12.6316E-03    -62.0743E-03
               7     12.0000E-30              0.];
           
idmap = [5 3 7 2 4 1 6]; 
error = abs(u_xy-AbaqusData(idmap,2:3));
assert(all(error(:)<1e-6),'benchmark AbaqusData not maching the results')

Answer 1

[Voss]

This line has only 11 elements:

6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3;

And this line has only 10 elements:

0, -12*E*Izz/Le(i)^3, 0 ,0, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 12*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3;

The other 10 lines each have 12 elements.