Index in position 1 exceeds array bounds (must not exceed 8).

I am trying to find the global stiffness matrix of 4 Q4 elements, not sure why im getting this error can someone help, kinda very noob at coding

this is my code:

%(a)
% Material properties
E = 207e9;  % Young's modulus in Pa
nu = 0.32;  % Poisson's ratio

% Pressure applied on the plate
P = 39 * 1e6;  % Pa

% Element properties
t = 20*10^-3;  % Plate thickness in meters

% Gauss quadrature points and weights
dof = 2;  % Number of integration points in each direction
xiVec = [-1/sqrt(3), 1/sqrt(3)];  % Gauss quadrature points in the xi direction
etaVec = xiVec;  % Gauss quadrature points in the eta direction
wVec = [1, 1];  % Gauss weights in both directions

% Initialize the element stiffness matrices
Ke1 = zeros(8, 8);
Ke2 = zeros(8, 8);
Ke3 = zeros(8, 8);
Ke4 = zeros(8, 8);

% Loop over the integration points
for i = 1:dof
    for j = 1:dof
        xi = xiVec(i);
        eta = etaVec(j);
        
        % Shape functions and derivatives
        N = 1/4 * [(1 - xi) * (1 - eta);
            (1 + xi) * (1 - eta);
            (1 + xi) * (1 + eta);
            (1 - xi) * (1 + eta)];
        dN_dxi = 1/4 * [-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)];
        dN_deta = 1/4 * [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi];
        
        % Jacobian matrix
        J = 10^-3 *[dN_dxi; dN_deta] * [0, -100; 100, -100; 100, 0; 0, 0];  % x1, y1, etc. are the nodal coordinates
        
        % Determinant of the Jacobian matrix
        detJ = det(J);
        
        % Inverse of the Jacobian matrix
        invJ = inv(J);
        
        % B matrix (strain-displacement matrix)
        B = zeros(3, 8);
        for k = 1:4
            B(1, 2*k-1) = dN_dxi(k);
            B(2, 2*k) = dN_deta(k);
            B(3, 2*k-1) = dN_deta(k);
            B(3, 2*k) = dN_dxi(k);
        end
        
        % Element stiffness matrix contribution
        C = E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2];
        Ke = t * B' * C * B * detJ;
        
        % Add the contribution to the respective element stiffness matrix
        if i == 1 && j == 1
            Ke1 = Ke1 + wVec(i) * wVec(j) * Ke;
        elseif i == 2 && j == 1
            Ke2 = Ke2 + wVec(i) * wVec(j) * Ke;
        elseif i == 2 && j == 2
            Ke3 = Ke3 + wVec(i) * wVec(j) * Ke;
        elseif i == 1 && j == 2
            Ke4 = Ke4 + wVec(i) * wVec(j) * Ke;
        end
    end
end

% Display the element stiffness matrices
disp('Element Stiffness Matrix for Element 1:');

Element Stiffness Matrix for Element 1:

disp(Ke1);

   1.0e+06 *

    2.4027    1.1834   -1.6297   -0.4559   -0.6438   -0.3171   -0.1292   -0.4104
    1.1834    2.4027   -0.4104   -0.1292   -0.3171   -0.6438   -0.4559   -1.6297
   -1.6297   -0.4104    1.8368   -0.3171    0.4367    0.1100   -0.6438    0.6175
   -0.4559   -0.1292   -0.3171    0.7384    0.1222    0.0346    0.6508   -0.6438
   -0.6438   -0.3171    0.4367    0.1222    0.1725    0.0850    0.0346    0.1100
   -0.3171   -0.6438    0.1100    0.0346    0.0850    0.1725    0.1222    0.4367
   -0.1292   -0.4559   -0.6438    0.6508    0.0346    0.1222    0.7384   -0.3171
   -0.4104   -1.6297    0.6175   -0.6438    0.1100    0.4367   -0.3171    1.8368

disp('Element Stiffness Matrix for Element 2:');

Element Stiffness Matrix for Element 2:

disp(Ke2);

   1.0e+06 *

    1.8368    0.3171   -1.6297    0.4104   -0.6438   -0.6175    0.4367   -0.1100
    0.3171    0.7384    0.4559   -0.1292   -0.6508   -0.6438   -0.1222    0.0346
   -1.6297    0.4559    2.4027   -1.1834   -0.1292    0.4104   -0.6438    0.3171
    0.4104   -0.1292   -1.1834    2.4027    0.4559   -1.6297    0.3171   -0.6438
   -0.6438   -0.6508   -0.1292    0.4559    0.7384    0.3171    0.0346   -0.1222
   -0.6175   -0.6438    0.4104   -1.6297    0.3171    1.8368   -0.1100    0.4367
    0.4367   -0.1222   -0.6438    0.3171    0.0346   -0.1100    0.1725   -0.0850
   -0.1100    0.0346    0.3171   -0.6438   -0.1222    0.4367   -0.0850    0.1725

disp('Element Stiffness Matrix for Element 3:');

Element Stiffness Matrix for Element 3:

disp(Ke3);

   1.0e+06 *

    0.1725    0.0850    0.0346    0.1100   -0.6438   -0.3171    0.4367    0.1222
    0.0850    0.1725    0.1222    0.4367   -0.3171   -0.6438    0.1100    0.0346
    0.0346    0.1222    0.7384   -0.3171   -0.1292   -0.4559   -0.6438    0.6508
    0.1100    0.4367   -0.3171    1.8368   -0.4104   -1.6297    0.6175   -0.6438
   -0.6438   -0.3171   -0.1292   -0.4104    2.4027    1.1834   -1.6297   -0.4559
   -0.3171   -0.6438   -0.4559   -1.6297    1.1834    2.4027   -0.4104   -0.1292
    0.4367    0.1100   -0.6438    0.6175   -1.6297   -0.4104    1.8368   -0.3171
    0.1222    0.0346    0.6508   -0.6438   -0.4559   -0.1292   -0.3171    0.7384

disp('Element Stiffness Matrix for Element 4:');

Element Stiffness Matrix for Element 4:

disp(Ke4);

   1.0e+06 *

    0.7384    0.3171    0.0346   -0.1222   -0.6438   -0.6508   -0.1292    0.4559
    0.3171    1.8368   -0.1100    0.4367   -0.6175   -0.6438    0.4104   -1.6297
    0.0346   -0.1100    0.1725   -0.0850    0.4367   -0.1222   -0.6438    0.3171
   -0.1222    0.4367   -0.0850    0.1725   -0.1100    0.0346    0.3171   -0.6438
   -0.6438   -0.6175    0.4367   -0.1100    1.8368    0.3171   -1.6297    0.4104
   -0.6508   -0.6438   -0.1222    0.0346    0.3171    0.7384    0.4559   -0.1292
   -0.1292    0.4104   -0.6438    0.3171   -1.6297    0.4559    2.4027   -1.1834
    0.4559   -1.6297    0.3171   -0.6438    0.4104   -0.1292   -1.1834    2.4027

% (b) AssembleThe code for assembling the global stiffness matrix is as follows:

% Total number of nodes and elements
numNodes = 9;  % Assuming 9 nodes in total
numElements = 4;  % Number of elements

% Degree of freedom (DOF) per node
DOFperNode = 2;

% Total number of DOFs
numDOFs = numNodes * DOFperNode;

% Initialize the global stiffness matrix
KGlob = zeros(18);

% Element connectivity matrix
Elconn = [1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8];  % Assuming the nodes are numbered in a specific order

% Assemble the global stiffness matrix
for element = 1:numElements
    % Get the element stiffness matrix for the current element
    if element == 1
        Ke = Ke1;
    elseif element == 2
        Ke = Ke2;
    elseif element == 3
        Ke = Ke3;
    else
        Ke = Ke4;
    end
    
    % Loop over the element nodes
    for i = 1:4
        % Get the global node number for the current element node
        node = Elconn(element, i);
        
        % Loop over the degrees of freedom (DOFs) of the current node
        for j = 1:DOFperNode
            % Get the global DOF number
            DOF = (node - 1) * DOFperNode + j
            
            % Add the contribution to the global stiffness matrix
            for k = 1:DOFperNode
                DOF_global = (node - 1) * DOFperNode + k;
                KGlob(DOF, DOF_global) = KGlob(DOF, DOF_global) + Ke(DOF, DOF_global);
            end
        end
    end
end

DOF = 1

DOF = 2

DOF = 3

DOF = 4

DOF = 9

Index in position 1 exceeds array bounds. Index must not exceed 8.

% Display the global stiffness matrix
disp('Global Stiffness Matrix:');
disp(KGlob);

% (c) Initialize the element force vectors
fe1 = zeros(8, 1);
fe2 = zeros(8, 1);
fe3 = zeros(8, 1);
fe4 = zeros(8, 1);

% Loop over the integration points
for i = 1:dof
    for j = 1:dof
        xi = xiVec(i);
        eta = etaVec(j);
        
        % Shape functions and derivatives
        N = 1/4 * [(1 - xi) * (1 - eta);
            (1 + xi) * (1 - eta);
            (1 + xi) * (1 + eta);
            (1 - xi) * (1 + eta)];
        
        % Jacobian matrix
        J = [dN_dxi; dN_deta] * [x1, y1; x2, y2; x3, y3; x4, y4];  % x1, y1, etc. are the nodal coordinates
        
        % Determinant of the Jacobian matrix
        detJ = det(J);
        
        % Element force vector contribution
        fe = P * N * detJ;
        
        % Add the contribution to the respective element force vector
        if i == 1 && j == 1
            fe1 = fe1 + wVec(i) * wVec(j) * fe;
        elseif i == 2 && j == 1
            fe2 = fe2 + wVec(i) * wVec(j) * fe;
        elseif i == 2 && j == 2
            fe3 = fe3 + wVec(i) * wVec(j) * fe;
        elseif i == 1 && j == 2
            fe4 = fe4 + wVec(i) * wVec(j) * fe;
        end
    end
end

% Display the element force vectors
disp('Element Force Vector for Element 1:');
disp(fe1);

disp('Element Force Vector for Element 2:');
disp(fe2);

disp('Element Force Vector for Element 3:');
disp(fe3);

disp('Element Force Vector for Element 4:');
disp(fe4);

% (D) Initialize the global force vector
F_global = zeros(numDOFs, 1);

% Assemble the global force vector
for element = 1:numElements
    % Get the element force vector for the current element
    if element == 1
        fe = fe1;
    elseif element == 2
        fe = fe2;
    elseif element == 3
        fe = fe3;
    else
        fe = fe4;
    end
    
    % Loop over the element nodes
    for i = 1:4
        % Get the global node number for the current element node
        node = Elconn(element, i);
        
        % Loop over the degrees of freedom (DOFs) of the current node
        for j = 1:DOFperNode
            % Get the global DOF number
            DOF = (node - 1) * DOFperNode + j;
            
            % Add the contribution to the global force vector
            F_global(DOF) = F_global(DOF) + fe(DOF);
        end
    end
end

% Display the global force vector
disp('Global Force Vector:');
disp(F_global);

% (e) Apply boundary conditions
% Assuming node 1 is fixed in both x and y directions
fixedDOFs = [1, 2];  % DOFs corresponding to node 1

% Set the corresponding rows and columns in the global stiffness matrix to zero, except for the diagonal entries
KGlob(fixedDOFs, :) = 0;
KGlob(:, fixedDOFs) = 0;
KGlob(fixedDOFs, fixedDOFs) = eye(length(fixedDOFs));

% Set the corresponding entries in the global force vector to zero
F_global(fixedDOFs) = 0;

% (f) Solve for the unknown displacements
U = KGlob \ F_global;

% Display the unknown displacements
disp('Unknown Displacements:');
disp(U);

% (g) Calculate the plane stress components in the center of each element
stress_center = zeros(numElements, 3);  % [sigma_xx, sigma_yy, tau_xy]

for element = 1:numElements
    % Get the element stiffness matrix for the current element
    if element == 1
        Ke = Ke1;
    elseif element == 2
        Ke = Ke2;
    elseif element == 3
        Ke = Ke3;
    else
        Ke = Ke4;
    end
    
    % Get the element displacements for the current element
    Ue = U(Elconn(element, :));
    
    % Calculate the strains using the strain-displacement relation
    strains = B * Ue;
    
    % Calculate the stress components using Hooke's law
    stresses = C * strains;
    
    % Store the stress components in the center of the element
    stress_center(element, :) = stresses(:, 1)';
end

% Display theplane stress components in the center of each element:
disp('Plane Stress Components in the Center of Each Element:');
disp(stress_center);

[Dyuman Joshi]

See the above edit, you should be able to spot where the error occurs, which is due to the value of variable DOF.

Answer 1

[James Tursa]

Here is how to debug these types of errors:

Type the following at the command prompt:

dbstop if error

Then run your code. When the error occurs, the program will pause with all variables intact. Examine the variables on the line with the error to determine which ones are causing the error. Then backtrack in your code to figure out why they are not the values and/or sizes you expected.