I face a problem with the dimensions of Z as it should be a 2x2 matrix.

%------Creating a strain matrix-----------%
E = zeros(2*(size(K1,1)-1),2*(size(K1,1)-1));
for j = 1: (size(K1,1)-1)
    for i = 1: (size(K1,1)-1)
        v = [zpk1(j,i);zpk2(j,i);zpk1(j,i+1);zpk2(j,i+1);zpk1(j+1,i+1);zpk2(j+1,i+1);zpk1(j+1,i);zpk2(j+1,i)];
        e1= strain1(v);
        e2= strain2(v);
        e3= strain3(v);
        e4= strain4(v);
        E((2*j-1),(2*i-1)) = e1(1);
        E((2*j-1),2*i) = e2(1);
        E(2*j,(2*i-1)) = e4(1);
        E(2*j,2*i) = e3(1);
    end
end


su = round(4*d/g)+1;
r = E((((su-1)- round(d/g)):-1:1),(su-1));
p = E((((su-1)- round(d/g)):-1:1),(su));
st = (p+r)./(1.5*(10)^(-4));
m = g*(size(r)-1);
H = 0:g:(m);
h = H./10;

writematrix(st)
type 'st.txt'
writematrix(h)
type 'h.txt'


plot(h,st)
xlim([0 1.5])
xlabel("X2/d",'fontsize',14)
ylabel("Normalised Strain along vertical path", 'fontsize',13)
legend("d/w = 0.5")
saveas(gcf,'plot.png')
contourf(st)
colorbar
saveas(gcf,'contour.png')

The final plot of the above code should resemble something like this:

I have double chekced, the formula is correct and the value of st at a given r and p match with the experimental results. However, I face a problem with the dimensions of Z as it should be a 2x2 matrix.

Thank you in advance.

[Arnav]

Hi @Mohd Babar,

Can you provide executable code with data that reproduces the issue?

[Mohd Babar]

clc;
clear all;

%-----Geometry, Number of simulations------%
d = 10;
N = 1;
a11 = 250 * (d / 10);
a12 = 100 * (d / 10);

%-----Mesh grid boundary definition--------%
xms = a11 / 2 - 2 * d;
xmf = a11 / 2 + 2 * d;
yms = a12 / 2 - 2 * d;
ymf = a12 / 2 + 2 * d;

%----------------Grid size----------------%
g = 0.085;
xq1 = xms:g:xmf;
xq = xq1';
yq1 = yms:g:ymf;
yq = yq1';

%-----Mesh grid around the hole----------%
[K1, K2] = meshgrid(xq, yq);
a = size(xq1, 2);
b = size(yq1, 2);
for k = 1:b
    for j = 1:a
        e = ((k * g) - 2 * d)^2 + ((j * g) - 2 * d)^2 - (0.5 * d)^2;
        if e < 0
            K1(j, k) = NaN;
            K2(j, k) = NaN;
        end
    end
end

zz1 = 0;
zz2 = 0;

%----Importing files from Abaqus---------%
for i = 1:N
    fx = ['d:/babar/Abaqus/S2000/aax' num2str(i) '.txt'];
    fy = ['d:/babar/Abaqus/S2000/aay' num2str(i) '.txt'];
    fu = ['d:/babar/Abaqus/S2000/aau1' num2str(i) '.txt'];
    fv = ['d:/babar/Abaqus/S2000/aau2' num2str(i) '.txt'];
    
    x = textfile(fx);
    y = textfile(fy);
    z1 = textfile(fu);
    z2 = textfile(fv);
    U1 = griddata(x, y, z1, K1, K2);
    U2 = griddata(x, y, z2, K1, K2);
    
    zz1 = U1 + zz1;
    zz2 = U2 + zz2;
end

%--------Averaging the displacements-------%
zpk1 = zz1 / N;
zpk2 = zz2 / N;

xlswrite('zpk1.xlsx', zpk1)
xlswrite('zpk2.xlsx', zpk2)

%------Creating a strain matrix-----------%
E = zeros(2 * (size(K1, 1) - 1), 2 * (size(K1, 1) - 1));
for j = 1:(size(K1, 1) - 1)
    for i = 1:(size(K1, 1) - 1)
        v = [zpk1(j, i); zpk2(j, i); zpk1(j, i + 1); zpk2(j, i + 1); zpk1(j + 1, i + 1); zpk2(j + 1, i + 1); zpk1(j + 1, i); zpk2(j + 1, i)];
        e1 = strain1(v);
        e2 = strain2(v);
        e3 = strain3(v);
        e4 = strain4(v);
        E((2 * j - 1), (2 * i - 1)) = e1(1);
        E((2 * j - 1), 2 * i) = e2(1);
        E(2 * j, (2 * i - 1)) = e4(1);
        E(2 * j, 2 * i) = e3(1);
    end
end

su = round(4 * d / g) + 1;
r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);

st = (p + r) / (1.5 * 10^(-4));
m = g * (size(r, 1) - 1);
H = 0:g:m;
h = H / 10;

writematrix(st, 'st.txt')
writematrix(h, 'h.txt')

figure;
plot(h, st)
xlim([0 1.5])
xlabel("X2/d", 'fontsize', 14)
ylabel("Normalized Strain along vertical path", 'fontsize', 13)
legend("d/w = 0.5")
saveas(gcf, 'plot.png')

figure;
contourf(st)
colorbar
saveas(gcf, 'contour.png')


function E1 = strain1(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain2(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain3(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain4(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function B21 = B2()
B21 = [0.05 / 0.0025 0 0 0; 0 0.05 / 0.0025 0 0; 0 0 0.05 / 0.0025 0; 0 0 0 0.05 / 0.0025];
end

function B31 = B3(x1, eta1)
B31 = [...
    -(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25 0; ...
    -(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25 0; ...
    0 -(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25; ...
    0 -(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25];
end

function value = textfile(filename)
fileID = fopen(filename, 'r');
if fileID == -1
    error('Failed to open file: %s', filename);
end
value = fscanf(fileID, '%f');
fclose(fileID);
end



%"These are the extract data from abaqus and
su = round(4 * d / g) + 1;
r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);
[r, p] = meshgrid(r,p);

st = (p + r) / (1.5 * 10^(-4));
m = g * (size(r, 1) - 1);
H = 0:g:m;
h = H / 10; 

i have tried with meshgrid option, after using meshgrid i rid of the error but results are not correct."

Thanks in advance

[Mohd Babar]

Here Z is "st"

Answer 1

[Walter Roberson]

d = 10;
g = 0.085;

d and g are scalars

su = round(4 * d / g) + 1;

su is formed from scalar d and g, so su is scalar

r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);

su-1 and su are scalars, so you are extracting exactly one column from E into r and p, so r and p are vectors (not 2D)

st = (p + r) / (1.5 * 10^(-4));

vector plus vector (in the same direction) is vector, so st is a vector.

contourf(st)

You are trying to contourf() a vector.

[Mohd Babar]

So, How i debug this error

[Walter Roberson]

I already went through the debugging steps and showed exactly what the cause of the problem is.

Your r and p are vectors, so your st is a vector, and you cannot contour a vector.

Your su needs to be a scalar for the colon operator of the indexing of r and p to work properly.

What we are left with is the possibility that you want non-scalar second dimension for r and p, or else that you want to contourf() something different than st.

[Mohd Babar]

[r, p] = meshgrid(r, p);

mehsgrid function is use for converte the vector to a matrix, after defining r and p we can use this it will resolve the error but plot is still incurrect

[Walter Roberson]

Using meshgrid like that is not going to work. You have

st = (p + r) / (1.5 * 10^(-4));

and making p and r into grids cannot possibly give you the kind of curved output that you hope for.

clc;
clear all;

%-----Geometry, Number of simulations------%
d = 10;
N = 1;
a11 = 250 * (d / 10);
a12 = 100 * (d / 10);

%-----Mesh grid boundary definition--------%
xms = a11 / 2 - 2 * d;
xmf = a11 / 2 + 2 * d;
yms = a12 / 2 - 2 * d;
ymf = a12 / 2 + 2 * d;

%----------------Grid size----------------%
g = 0.085;
xq1 = xms:g:xmf;
xq = xq1';
yq1 = yms:g:ymf;
yq = yq1';

%-----Mesh grid around the hole----------%
[K1, K2] = meshgrid(xq, yq);
a = size(xq1, 2);
b = size(yq1, 2);
for k = 1:b
    for j = 1:a
        e = ((k * g) - 2 * d)^2 + ((j * g) - 2 * d)^2 - (0.5 * d)^2;
        if e < 0
            K1(j, k) = NaN;
            K2(j, k) = NaN;
        end
    end
end

zz1 = 0;
zz2 = 0;

%----Importing files from Abaqus---------%
for i = 1:N
    fx = ['aax' num2str(i) '.txt'];
    fy = ['aay' num2str(i) '.txt'];
    fu = ['aau1' num2str(i) '.txt'];
    fv = ['aau2' num2str(i) '.txt'];
    
    x = textfile(fx);
    y = textfile(fy);
    z1 = textfile(fu);
    z2 = textfile(fv);
    U1 = griddata(x, y, z1, K1, K2);
    U2 = griddata(x, y, z2, K1, K2);
    
    zz1 = U1 + zz1;
    zz2 = U2 + zz2;
end

%--------Averaging the displacements-------%
zpk1 = zz1 / N;
zpk2 = zz2 / N;

xlswrite('zpk1.xlsx', zpk1)

Warning: Unable to write to Excel format, attempting to write file to CSV format.

xlswrite('zpk2.xlsx', zpk2)

Warning: Unable to write to Excel format, attempting to write file to CSV format.

%------Creating a strain matrix-----------%
E = zeros(2 * (size(K1, 1) - 1), 2 * (size(K1, 1) - 1));
for j = 1:(size(K1, 1) - 1)
    for i = 1:(size(K1, 1) - 1)
        v = [zpk1(j, i); zpk2(j, i); zpk1(j, i + 1); zpk2(j, i + 1); zpk1(j + 1, i + 1); zpk2(j + 1, i + 1); zpk1(j + 1, i); zpk2(j + 1, i)];
        e1 = strain1(v);
        e2 = strain2(v);
        e3 = strain3(v);
        e4 = strain4(v);
        E((2 * j - 1), (2 * i - 1)) = e1(1);
        E((2 * j - 1), 2 * i) = e2(1);
        E(2 * j, (2 * i - 1)) = e4(1);
        E(2 * j, 2 * i) = e3(1);
    end
end

su = round(4 * d / g) + 1;
r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);
[r, p] = meshgrid(r,p);

st = (p + r) / (1.5 * 10^(-4));
m = g * (size(r, 1) - 1);
H = 0:g:m;
h = H / 10;

writematrix(st, 'st.txt')
writematrix(h, 'h.txt')

figure;
plot(h, st)
xlim([0 1.5])
xlabel("X2/d", 'fontsize', 14)
ylabel("Normalized Strain along vertical path", 'fontsize', 13)
legend("d/w = 0.5")
saveas(gcf, 'plot.png')

figure;
contourf(st)
colorbar
saveas(gcf, 'contour.png')

function E1 = strain1(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain2(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain3(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function E1 = strain4(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end

function B21 = B2()
B21 = [0.05 / 0.0025 0 0 0; 0 0.05 / 0.0025 0 0; 0 0 0.05 / 0.0025 0; 0 0 0 0.05 / 0.0025];
end

function B31 = B3(x1, eta1)
B31 = [...
    -(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25 0; ...
    -(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25 0; ...
    0 -(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25; ...
    0 -(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25];
end

function value = textfile(filename)
fileID = fopen(filename, 'r');
if fileID == -1
    error('Failed to open file: %s', filename);
end
value = fscanf(fileID, '%f');
fclose(fileID);
end

[Walter Roberson]

Maybe you want to use r and p and some other variable as x, y, z vectors together with scatteredInterpolant() or gridddata()