Efficiency and alternatives for for loops

Good afternoon,

Recently I have been trying to plot multiple 2D plots of a yield surface equation, the goal was to be able to show the multiple parameters. The problem I am having is the enormous amount of for loops I am having to run, which im now finding our very expensive (Im relatively new). Could someone help me find an alternative? would declaring a matrix of data points rather then running a for loop be faster?

Heres a sample of the code, for refrence I would have to repeat this around 5 times for 1 surface.

Is there a way to have everything within the first 3 for loops? I have tried to puzzle my way around this but nothing seems to work.

Thank you for the time!

Kevin

clc; clear; close all;

tic;


    Yieldstr = 1000; %KPa 
    
f = 0.8; g = 0.6; h = 0.4;
f1 = 0.2; g1 = 0.4; h1 = 0.6;
f2 = 0.8; g2 = 0.4; h2 = 0.2;

l = 0.5; m = 0.6; n = 0.4;
l1 = 0.4; m1 = 0.2; n1 = 0.7;
l2 = 0.8; m2 = 0.1; n2 = 1.0;

PStress = zeros(134343,3);
PStress12 = zeros(134343,2);
PStress23 = zeros(134343,2);
Shear12 = zeros(42334,2);
Shear23 = zeros(42334,2);

minValue = -4*Yieldstr;
maxValue = 4*Yieldstr;
fidelity = Yieldstr/125;
i=1;
j=1;
k=1;
  
 for sigx = minValue:fidelity:maxValue
    for sigy = minValue:fidelity:maxValue
        for sigz = 0:0
            
            Hill12 = sqrt(((1/2).*((g+h).*sigx.^2)) + ((1/2).*((f+h).*sigy.^2)) - h.*sigx.*sigy);
         
             if Hill12 < Yieldstr
                if (sigx==sigy) 
                %stops code
            else
                PStress12(i,1) = sigx;
                PStress12(i,2) = sigy;
                
                i = i + 1;
                
                end
            end
        end
    end
 end
 
  for sigx = minValue:fidelity:maxValue
    for sigy = minValue:fidelity:maxValue
        for sigz = 0:0
            
            Hill121 = sqrt(((1/2).*((g1+h1).*sigx.^2)) + ((1/2).*((f1+h1).*sigy.^2)) - h1.*sigx.*sigy);
            Hill122 = sqrt(((1/2).*((g2+h2).*sigx.^2)) + ((1/2).*((f2+h2).*sigy.^2)) - h2.*sigx.*sigy);
         
             if Hill121 < Yieldstr
                if (sigx==sigy) 
                %stops code
            else
                PStress121(j,1) = sigx;
                PStress121(j,2) = sigy;
                
                j = j + 1;
                
                end
            end
        end
    end
 end
 
  for sigx = minValue:fidelity:maxValue
    for sigy = minValue:fidelity:maxValue
        for sigz = 0:0
            
            Hill122 = sqrt(((1/2).*((g2+h2).*sigx.^2)) + ((1/2).*((f2+h2).*sigy.^2)) - h2.*sigx.*sigy);
         
             if Hill122 < Yieldstr
                if (sigx==sigy) 
                %stops code
            else
                PStress122(k,1) = sigx;
                PStress122(k,2) = sigy;
                
                k = k + 1;
                
                end
            end
        end
    end
 end
 
 
 G = boundary(PStress12,0);
 plot(PStress12(G,1),PStress12(G,2),'LineWidth',1.75);
 hold on
 G1 = boundary(PStress121,0);
 plot(PStress121(G1,1),PStress121(G1,2),'LineWidth',1.75);
 G2 = boundary(PStress122,0);
 plot(PStress122(G2,1),PStress122(G2,2),'LineWidth',1.75);
 
hold off

axis('tight')
title('\sigma_y vs \sigma_x')
legend('f = 0.8, g = 0.6, h = 0.4','f = 0.2, g = 0.4, h = 0.6','f = 0.8, g = 0.4, h = 0.2','Location',"best")
xlabel('\sigma_x (KPa)');
ylabel('\sigma_y (KPa)');

Answer 1

[DGM]

I'm not sure if this is going to give you exactly what you want. The following simplified example will produce the same boundary curve, but the output arrays will be shorter, since they only include the points within the ROI. I'm only doing the first loop for sake of clarity, but the other should be very similar.

clc; clf; clearvars

Yieldstr = 1000; %KPa

f = 0.8; g = 0.6; h = 0.4;
f1 = 0.2; g1 = 0.4; h1 = 0.6;
f2 = 0.8; g2 = 0.4; h2 = 0.2;
l = 0.5; m = 0.6; n = 0.4;
l1 = 0.4; m1 = 0.2; n1 = 0.7;
l2 = 0.8; m2 = 0.1; n2 = 1.0;

kw = 4; % 4
kf = 125; % 125
minValue = -kw*Yieldstr;
maxValue = kw*Yieldstr;
fidelity = Yieldstr/kf;


% OLD METHOD 
tic
PStress12 = zeros(134343,2);

i=1;

for sigx = minValue:fidelity:maxValue
    for sigy = minValue:fidelity:maxValue
        for sigz = 0:0
            
            Hill12 = sqrt(((1/2).*((g+h).*sigx.^2)) + ((1/2).*((f+h).*sigy.^2)) - h.*sigx.*sigy);
            
            if Hill12 < Yieldstr
                if (sigx==sigy)
                    %stops code
                else
                    PStress12(i,1) = sigx;
                    PStress12(i,2) = sigy;
                    
                    i = i + 1;
                    
                end
            end
        end
    end
end
toc

Elapsed time is 0.040100 seconds.

G_orig = boundary(PStress12,0);


% ALT METHOD
tic
sigx = minValue:fidelity:maxValue;
sigy = sigx.'; % if sigx,sigy are never unique, this may simplify further
Hill12 = sqrt(((1/2).*((g+h).*sigx.^2)) + ((1/2).*((f+h).*sigy.^2)) - h.*sigx.*sigy);
[m,n,~] = find(Hill12 < Yieldstr);
PStress12 = [sigx(n).' sigy(m)]; % note we don't need to preallocate this array anymore
PStress12 = PStress12(PStress12(:,1)~=PStress12(:,2),:); % idk why this is necessary
toc

Elapsed time is 0.026506 seconds.

G = boundary(PStress12,0);
plot(PStress12(G,1),PStress12(G,2),'LineWidth',1.75);

immse(G_orig,G) % show that the two boundary curves are identical

ans = 0

I imagine there are yet ways to expedite things.

[Kevin Holly]

This is not related to the question, but I was just messing around with the data and drew some pictures. I thought that I would share.

Yieldstr = 1000; %KPa

f = 0.8; g = 0.6; h = 0.4;
f1 = 0.2; g1 = 0.4; h1 = 0.6;
f2 = 0.8; g2 = 0.4; h2 = 0.2;
l = 0.5; m = 0.6; n = 0.4;
l1 = 0.4; m1 = 0.2; n1 = 0.7;
l2 = 0.8; m2 = 0.1; n2 = 1.0;

kw = 4; % 4
kf = 125; % 125
minValue = -kw*Yieldstr;
maxValue = kw*Yieldstr;
fidelity = Yieldstr/kf;

sigx = minValue:fidelity:maxValue;
sigy = sigx.';
Hill12 = sqrt(((1/2).*((g+h).*sigx.^2)) + ((1/2).*((f+h).*sigy.^2)) - h.*sigx.*sigy);

[xg, yg] = meshgrid(sigx,sigx);

figure
subplot(2,2,1)
surf(xg',yg',Hill12)
lol = gca;
lol.Children.EdgeColor = 'none';

hold on
surf(xlim,ylim,Yieldstr*ones(2),'FaceColor','k','FaceAlpha',.5)

ZL = zlim;
subplot(2,2,2)
[c,h] = contourf(xg',yg',Hill12, [ZL(1) Yieldstr]);
subplot(2,2,3)
scatter(c(1,2:end),c(2,2:end),'.r')
subplot(2,2,4)
plot(c(1,2:end),c(2,2:end))

[Kevin Hanekom]

That contourf plot looks amazing, thank you both for these amazing ideas!! I think doing it your way @DGM will allow me to write each equation in "one loop" rather then the 3 for loops i needed for each equation resulting in a final 15 for loops (old way).

Thank you again!