Specifying an element and gaussian point for a script

I was wondering how exactly I would specify an individual element and gaussian point for the following:

function res = FEMbmatrix(L,n) %inputs n[2,2]=4 i.e. number of plates. L is length (2,2] squr

nx = n(1) + 1; % number of nodes = number of elements +1 (x) ;num. of nodes given num. element input established
ny = n(2) + 1; % number of nodes = number of elements +1 ? ; n will be a 1*2 matrix; n = [x,y]
nNodes = nx * ny; % total number of nodes
nE = n(1) * n(2); % total number of elements; e.g n[2,2]=4 elements
dx = L(1)/n(1); % length of x/number of elements in x
dy = L(2)/n(2); % height/number of element in y direction

% Gaussian points
Zeta = [-(1/3)^0.5; -(1/3)^0.5; (1/3)^0.5; (1/3)^0.5]; % local zeta; bottom left quadrant to bottom right quadrant clockwise.
Eta = [-(1/3)^0.5; (1/3)^0.5; (1/3)^0.5; -(1/3)^0.5]; % local eta;

eID = 1; %loop over all elemets
for j = 1:n(2)
    for i = 1:n(1)
        %node ID
        e(eID).n(1) = (j-1)*nx+i;
        e(eID).n(2) = j*nx+i;
        e(eID).n(3) = j*nx+i+1;
        e(eID).n(4) = (j-1)*nx+i+1;
        %global node position
        e(eID).x(1,1) = (i-1)*dx;
        e(eID).y(1,1) = (j-1)*dy;
        e(eID).x(2,1) = (i-1)*dx;
        e(eID).y(2,1) = j*dy;
        e(eID).x(3,1) = i*dx;
        e(eID).y(3,1) = j*dy;
        e(eID).x(4,1) = i*dx;
        e(eID).y(4,1) = (j-1)*dy;
        eID=eID+1;
    end
end


for eID = 1:nE
   
    for gp=1:4                 %loop over all gaussian points
        zeta = Zeta(gp,1);
        eta = Eta(gp,1);
        dNdzeta(1) = -0.25*(1-eta);
        dNdzeta(2) = -0.25*(1+eta);
        dNdzeta(3) = 0.25*(1+eta);
        dNdzeta(4) = 0.25*(1-eta);
        dNdeta(1) = -0.25*(1-zeta);
        dNdeta(2) = 0.25*(1-zeta);
        dNdeta(3) = 0.25*(1+zeta);
        dNdeta(4) = -0.25*(1+zeta);
        dNdlocal = [dNdzeta; dNdeta];
        %Jacobian
        J(1,1) = dNdzeta*e(eID).x;
        J(1,2) = dNdzeta*e(eID).y;
        J(2,1) = dNdeta*e(eID).x;
        J(2,2) = dNdeta*e(eID).y;
        dNdglobal = J\dNdlocal; %compacted matrix 
        %define 3*8 b matrix
        B(1,1) = dNdglobal(1,1);
        B(1,2) = 0;
        B(1,3) = dNdglobal(1,2);
        B(1,4) = 0;
        B(1,5) = dNdglobal(1,3);
        B(1,6) = 0;
        B(1,7) = dNdglobal(1,3);
        B(1,8) = 0;
        B(2,1) = 0;
        B(2,2) = dNdglobal(2,1);
        B(2,3) = 0;
        B(2,4) = dNdglobal(2,2);
        B(2,5) = 0;
        B(2,6) = dNdglobal(2,3);
        B(2,7) = 0;
        B(3,1) = dNdglobal(2,1);
        B(3,2) = dNdglobal(1,1);
        B(3,3) = dNdglobal(2,2);
        B(3,4) = dNdglobal(1,2);
        B(3,5) = dNdglobal(2,3);
        B(3,6) = dNdglobal(1,3);
        B(3,7) = dNdglobal(2,4);
        B(3,8) = dNdglobal(1,4);
        B    
    
    end
   
end

I have written this script for the B matrix for each point. Ideally, I would want to store each iteration of the gaussian point for each element somewhere.

Thanks!

Answer 1

[Vedant Shah]

Hi @Thomas Stokes,

To store the ‘B’ matrix corresponding to each Gaussian point within each finite element, a cell array can be effectively utilized. This approach allows for organized and efficient access to the ‘B’ matrices.

In this context, the ‘B’ matrix should be stored for every combination of element ID ‘eID’ and Gaussian point ‘gp’. To implement this, the following modifications can be made to the code:

Preallocate the cell array before entering the element and Gaussian point loops:

% Preallocate B matrix storage 
Bmatrices = cell(nE, 4); 

Store the computed B matrix within the inner loop, after its calculation:

% Store B matrix 
Bmatrices{eID, gp} = B; 

This ensures that each B matrix is retained and can be referenced later using its corresponding element and Gaussian point indices.

For more information, refer to the following documentation:

• https://www.mathworks.com/help/matlab/ref/cell.html