Need some clarity on parfor loop.

I am try to create a FEM code in matlab using parallel processing. The following is my code.

%% bar Problem
% E - Youngs Modulus; l - length of the bar; A - Cross sectional area
% k local stiffness matrix; and K - Global stiffness matrix

clear();
clc;
tic;
syms x 
l = 300;                 % mm
E = 210*10^3;            % Mpa
A = 100 ;                % sqmm
p = 100 ;                 % N/mm
ne = 3;
nn = 4;
n_dof = 1;
ne_dof = 2;
nd = nn*n_dof;
le = l/ne;
N = [(le - x)/le , x/le];
Nx = diff(N,x);
k =double(int(E*A*Nx.'*Nx,x,0,le));
f = double(int(p*N.',x,0,le));
% initalize global stiffness matrixs and Force matrix to zero


% connectivity matrix 
c = [1,2;2,3;3,4;4,5;5,6];
K =  zeros(nd); 
F = zeros(nd,1);

parfor e = 1:ne
    for i =1:ne_dof
        for j = 1:ne_dof
            
            K(c(e,i),c(e,j)) =  K(c(e,i),c(e,j)) + k(i,j);
            
        end
        F(c(e,i)) =  F(c(e,i)) + f(i);
    end
end

% apply boundary condtion
% when x = 0, u1 = 0

K1 = K(2:end,2:end);
F1 = F(2:end);

X = linsolve(K1,F1);
toc;

I am geting "valid indices for 'K' are restricted in parfor" error on the statement inside the inner most and the outer one. I have tried multiple way to rectify it, but not so successful.

Kindly guide me to sort this out.

Regards

kavi

Answer 1

[Ruchika]

Hi, the error you are encountering, "valid indices for 'K' are restricted in parfor," occurs because the 'parfor' loop in MATLAB requires specific conditions to be met in order to execute in parallel. One of these conditions is that the loop indices must be known at the start of the loop, which is not the case when you are updating K and F inside the nested loops.

To resolve this issue, you can use temporary variables to accumulate the values in each iteration of the 'parfor' loop and then update K and F outside the loop. Here's an updated version of your code that should work:

%% Bar Problem
% E - Young's Modulus; l - length of the bar; A - Cross-sectional area
% k local stiffness matrix; and K - Global stiffness matrix
clear();
clc();
tic();

syms x;
l = 300;                 % mm
E = 210*10^3;            % MPa
A = 100;                 % sqmm
p = 100;                 % N/mm
ne = 3;
nn = 4;
n_dof = 1;
ne_dof = 2;
nd = nn * n_dof;
le = l / ne;
N = [(le - x) / le, x / le];
Nx = diff(N, x);
k = double(int(E * A * Nx.' * Nx, x, 0, le));
f = double(int(p * N.', x, 0, le));

% Initialize global stiffness matrix and force vector to zero
K = zeros(nd);
F = zeros(nd, 1);

% Connectivity matrix
c = [1, 2; 2, 3; 3, 4; 4, 5; 5, 6];

% Temporary matrices for accumulation
temp_K = zeros(nd);
temp_F = zeros(nd, 1);

parfor e = 1:ne
    temp_K_e = zeros(nd);
    temp_F_e = zeros(nd, 1);
    
    for i = 1:ne_dof
        for j = 1:ne_dof
            temp_K_e(c(e, i), c(e, j)) = temp_K_e(c(e, i), c(e, j)) + k(i, j);
        end
        temp_F_e(c(e, i)) = temp_F_e(c(e, i)) + f(i);
    end
    
    % Accumulate values in temporary matrices
    temp_K = temp_K + temp_K_e;
    temp_F = temp_F + temp_F_e;
end

Starting parallel pool (parpool) using the 'Processes' profile ...
Parallel pool using the 'Processes' profile is shutting down.

% Update K and F outside the parfor loop
K = temp_K;
F = temp_F;

% Apply boundary conditions
K1 = K(2:end, 2:end);
F1 = F(2:end);
X = linsolve(K1, F1);

toc();

Elapsed time is 22.088309 seconds.

By using temporary variables, you avoid the issue of updating K and F inside the 'parfor' loop, allowing the code to execute in parallel without any indexing errors.