for loop return empty struct

syms x y z

thickness = 0.1:0.01:0.30;
unitcellsize = 1:0.1:6;
num = max(size(unitcellsize));

a = 0:step:unitcellsize;
b = 0:step:unitcellsize;
c = 0:step:unitcellsize;
[X,Y,Z] = meshgrid(a,b,c);

f = @(x,y,z,unitcellsize) 2.*(cos((2.*pi.*x)./unitcellsize).* cos((2.*pi.*y)./unitcellsize) + cos((2.*pi.*y)./unitcellsize).* cos((2.*pi.*z)./unitcellsize) + cos((2.*pi.*z)./unitcellsize).*cos((2.*pi.*x)./unitcellsize)) - (cos((4.*pi.*x)./unitcellsize)+cos((4.*pi.*y)./unitcellsize)+cos((4.*pi.*z)./unitcellsize));

for ii = 1:num;
    
    V{ii} = f(X,Y,Z,unitcellsize(ii));

end

for i = 1:num;
    
    co(i) = isosurface(a,b,c,V{i},0) ;

end

It turns out that 'co' not return all the loop it stoped some how. I don't know why ?

Answer 1

[Walter Roberson]

syms x y z

thickness = 0.1:0.01:0.30;
unitcellsize = 1:0.1:6;
num = max(size(unitcellsize));
step = 0.5;
a = 0:step:unitcellsize;
b = 0:step:unitcellsize;
c = 0:step:unitcellsize;
[X,Y,Z] = meshgrid(a,b,c);

f = @(x,y,z,unitcellsize) 2.*(cos((2.*pi.*x)./unitcellsize).* cos((2.*pi.*y)./unitcellsize) + cos((2.*pi.*y)./unitcellsize).* cos((2.*pi.*z)./unitcellsize) + cos((2.*pi.*z)./unitcellsize).*cos((2.*pi.*x)./unitcellsize)) - (cos((4.*pi.*x)./unitcellsize)+cos((4.*pi.*y)./unitcellsize)+cos((4.*pi.*z)./unitcellsize));

for ii = 1:num
    V{ii} = f(X,Y,Z,unitcellsize(ii));
end

for i = 1:num
    thisco = isosurface(a,b,c,V{i},0);
    co(i) = thisco;
    if isempty(thisco) || isempty(thisco.vertices)
        minV = min(V{i},[],'all');
        maxV = max(V{i},[],'all');
        fprintf('iteration #%d, empty isosurface, V range is %f to %f\n', i, minV, maxV);
        break
    end
end

iteration #26, empty isosurface, V range is 0.010885 to 3.000000

Which is to say that you are getting empty structure entries because you are asking for isosurface at level 0 for data that does not include 0.

You are dividing by larger and larger numbers as you go, so you should expect that your f values will at some point become only slightly wavy around some constant offset.

[Teerapong Poltue]

Well, how can I do then to find surface area of each 3d plot as unitcellsize changed.

Previously I extract vertices from isosurface then using

A = co(j).vertices(co(j).faces(:, 2), :) - co(j).vertices(co(j).faces(:, 1), :);
B = co(j).vertices(co(j).faces(:, 3), :) - co(j).vertices(co(j).faces(:, 1), :);
C = cross(A, B, 2);
area(j) = 1/2 * sum(sqrt(sum(C.^2, 2)));

for finding surface area.

Then how can I do to find surface area of each 3d plot as unitcellsize changed.

[Walter Roberson]

a = 0:step:unitcellsize;
b = 0:step:unitcellsize;
c = 0:step:unitcellsize;

but unitcellsize is a vector, and when you use a vector as a bound, the first element of the vector is used. So a and b and c are not increasing as unitcellsize increases: they are stuck at 1. Is that what you want?

[Teerapong Poltue]

I want unit cell to be able to change the size so I can find different surface area as unit cell size change.