expanded matrix wont be bigger than 8 by 8

why wont K_global be bigger than 8 by 8?

this is my code:

clc; clear; clear all;

% Topology and element matrix
topologydegree = [1 2 3 4; 3 4 5 6; 5 6 7 8;
    7 8 9 10; 9 10 11 12; 11 12 13 14;
    13 14 15 16; 1 2 15 16; 3 4 15 16;
    3 4 13 14; 5 6 13 14; 7 8 13 14; 7 8 11 12]

topologydegree = 13×4

1 2 3 4 3 4 5 6 5 6 7 8 7 8 9 10 9 10 11 12 11 12 13 14 13 14 15 16 1 2 15 16 3 4 15 16 3 4 13 14

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

element = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8;
    1 8; 2 8; 2 7; 3 7; 4 7; 4 6]

element = 13×2

1 2 2 3 3 4 4 5 5 6 6 7 7 8 1 8 2 8 2 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Number of nodes
unique_elements = unique(element);
nodemax = numel(unique_elements)

nodemax = 8

% Number of elements
n = size(topologydegree, 1);

% Given variables
E = 24.17;
A = 1;
L = [300*sqrt(2); 300; 300; 300*sqrt(2);
    300; 300; 300; 300; 300; 300*sqrt(2);
    300; 300*sqrt(2); 300]

L = 13×1

424.2641 300.0000 300.0000 424.2641 300.0000 300.0000 300.0000 300.0000 300.0000 424.2641

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

co = [0, 0; 300, 300; 600, 300; 
    900, 300; 1200, 0; 900, 0; 600, 0; 300, 0]

co = 8×2

0 0 300 300 600 300 900 300 1200 0 900 0 600 0 300 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Initialize the global stiffness matrix
K_global = zeros(nodemax, nodemax);

for i = 1:n
    % Stiffness matrix for the current element
    Ke(:,:,i) = (E * A) / L(i) * [1, -1; -1, 1]
    
    % Calculate nxx and nyx
    x1 = co(element(i,1), 1);
    x2 = co(element(i,2), 1);
    y1 = co(element(i,1), 2);
    y2 = co(element(i,2), 2);
    
    nxx = (x2 - x1) / L(i);
    nyx = (y2 - y1) / L(i);
    
    % Geometry matrix for the current element
    G(:,:,i) = [nxx, nyx, 0, 0; 0, 0, nxx, nyx];
    G_(:,:,i) = G(:,:,i)';
    
    % Stiffness matrix in global coordinates
    Kea(:,:,i) = G_(:,:,i) * Ke(:,:,i) * G(:,:,i)
    
    % Assemble the global stiffness matrix
    for ii = 1:size(topologydegree, 2)
        for jj = 1:size(topologydegree, 2)
            row = topologydegree(i, ii);
            col = topologydegree(i, jj);
            K_global(row, col) = K_global(row, col) + Kea(ii, jj);
        end
    end
end

Ke = 2×2

0.0570 -0.0570 -0.0570 0.0570

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kea = 4×4

0.0285 0.0285 -0.0285 -0.0285 0.0285 0.0285 -0.0285 -0.0285 -0.0285 -0.0285 0.0285 0.0285 -0.0285 -0.0285 0.0285 0.0285

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Ke = 
Ke(:,:,1) =

    0.0570   -0.0570
   -0.0570    0.0570


Ke(:,:,2) =

    0.0806   -0.0806
   -0.0806    0.0806

Kea = 
Kea(:,:,1) =

    0.0285    0.0285   -0.0285   -0.0285
    0.0285    0.0285   -0.0285   -0.0285
   -0.0285   -0.0285    0.0285    0.0285
   -0.0285   -0.0285    0.0285    0.0285


Kea(:,:,2) =

    0.0806         0   -0.0806         0
         0         0         0         0
   -0.0806         0    0.0806         0
         0         0         0         0

Ke = 
Ke(:,:,1) =

    0.0570   -0.0570
   -0.0570    0.0570


Ke(:,:,2) =

    0.0806   -0.0806
   -0.0806    0.0806


Ke(:,:,3) =

    0.0806   -0.0806
   -0.0806    0.0806

Kea = 
Kea(:,:,1) =

    0.0285    0.0285   -0.0285   -0.0285
    0.0285    0.0285   -0.0285   -0.0285
   -0.0285   -0.0285    0.0285    0.0285
   -0.0285   -0.0285    0.0285    0.0285


Kea(:,:,2) =

    0.0806         0   -0.0806         0
         0         0         0         0
   -0.0806         0    0.0806         0
         0         0         0         0


Kea(:,:,3) =

    0.0806         0   -0.0806         0
         0         0         0         0
   -0.0806         0    0.0806         0
         0         0         0         0

Ke = 
Ke(:,:,1) =

    0.0570   -0.0570
   -0.0570    0.0570


Ke(:,:,2) =

    0.0806   -0.0806
   -0.0806    0.0806


Ke(:,:,3) =

    0.0806   -0.0806
   -0.0806    0.0806


Ke(:,:,4) =

    0.0570   -0.0570
   -0.0570    0.0570

Kea = 
Kea(:,:,1) =

    0.0285    0.0285   -0.0285   -0.0285
    0.0285    0.0285   -0.0285   -0.0285
   -0.0285   -0.0285    0.0285    0.0285
   -0.0285   -0.0285    0.0285    0.0285


Kea(:,:,2) =

    0.0806         0   -0.0806         0
         0         0         0         0
   -0.0806         0    0.0806         0
         0         0         0         0


Kea(:,:,3) =

    0.0806         0   -0.0806         0
         0         0         0         0
   -0.0806         0    0.0806         0
         0         0         0         0


Kea(:,:,4) =

    0.0285   -0.0285   -0.0285    0.0285
   -0.0285    0.0285    0.0285   -0.0285
   -0.0285    0.0285    0.0285   -0.0285
    0.0285   -0.0285   -0.0285    0.0285

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

disp('Global Stiffness Matrix K:')
disp(K_global)

Answer 1

[Walter Roberson]

unique_elements = unique(element);
nodemax = numel(unique_elements)

nodemax is 8.

K_global = zeros(nodemax, nodemax);

You initialize K_global as 8 x 8

            row = topologydegree(i, ii);
            col = topologydegree(i, jj);
            K_global(row, col) = K_global(row, col) + Kea(ii, jj);

Although row and col can each be larger than 8, you are adding a value to K_global(row,col) which is a problem when row or col exceeds 8 as K_global does not exist for row or col greater than 8.

Answer 2

[Swastik Sarkar]

Hi @Olivier,

In the provided code, K_global is being initialized to zeros based on the maximum of unique values from the element array, as seen in the following lines:

% Number of nodes
unique_elements = unique(element);
nodemax = numel(unique_elements)
K_global = zeros(nodemax, nodemax);

However, K_global is accessed using values from topologydegree instead of element:

row = topologydegree(i, ii);
col = topologydegree(i, jj);
K_global(row, col) = K_global(row, col) + Kea(ii, jj);

To align with the requirements, nodemax should be set to max(topologydegree) instead. This adjustment ensures that K_global is sized appropriately for the indices used.

Hope this helps.