Issues with plotting/best method

Greetings,

I would like to plot the following, but I get a blank graph and I am not sure if it is the way I am inputting the values

any recommendations? or a different approach?

close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio
%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31
%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter
%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));
%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));
% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
    Q12,Q22,0;
    0,0,Q66];
%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];
%%Laminate rotations in radians
% theta=pi/180*[0;45;-45;90;30;60;-30;-60]';%%DEGREES TO RADIANS
% theta=[theta;theta(end:-1:1)];
z=(-length(theta)/2+(0:length(theta)))'*h;
% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(theta)
    m=cosd(theta(i));
    n=sind(theta(i));
    T=[m.^2,n.^2,2.*m.*n;
        n.^2,m.^2,-2.*m.*n;
        -m.*n,m.*n,m.^2-n.^2];
      
    % Global transformed reduced stiffness coefficients Q_bar
    Q_bar{i}=inv(T).*Q.*T;
end
Aij=0;
Bij=0;
Dij=0;
for j=1:length(theta)
    A_j=Q_bar{j}*(z(j+1)-z(j));
    Aj{j}=A_j;
    Aij=Aij+Aj{j};
    B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
    Bj{j}=B_j;
    Bij=Bij+Bj{j};
    D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
    Dj{j}=D_j;
    Dij=Dij+Dj{j};
end
A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;
%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];
%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];
%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];
%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];
%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];
%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];
test=Q_bar{1}*shearf+z(1,1)*Q_bar{1}*bendingf;
test1=Q_bar{1}.*shearf+z(2,1)*Q_bar{1}.*bendingf;
% plot(test(1,1), z(1,1));
for k=1:length(theta)
 stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
 stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
 stress1f{k}=stress1; %Bottom
 stress2f{k}=stress2; %Top
end

b0=stress1f{1};
b1=stress1f{2};
b2=stress1f{3};
b3=stress1f{4};
b4=stress1f{5};
b5=stress1f{6};
b6=stress1f{1};
b7=stress1f{7};
b8=stress1f{8};
b9=stress1f{9};
b10=stress1f{10};
b11=stress1f{11};
b12=stress1f{12};
b13=stress1f{13};
b14=stress1f{14};
b15=stress1f{15};
t0=stress2f{1};
t1=stress2f{2};
t2=stress2f{3};
t3=stress2f{4};
t4=stress2f{5};
t5=stress2f{6};
t6=stress2f{1};
t7=stress2f{7};
t8=stress2f{8};
t9=stress2f{9};
t10=stress2f{10};
t11=stress2f{11};
t12=stress2f{12};
t13=stress2f{13};
t14=stress2f{14};
t15=stress2f{15};
plot(b0(1,1), z(1,1))
keyboard
hold on
plot(t0(1,1), z(2,1))
plot(b1(1,1), z(2,1))
plot(t1(1,1), z(3,1))
plot(b2(1,1), z(3,1))
plot(t2(1,1), z(4,1))
plot(b3(1,1), z(4,1))
plot(t3(1,1), z(5,1))
plot(b4(1,1), z(5,1))
plot(t4(1,1), z(6,1))
plot(b5(1,1), z(6,1))
plot(t5(1,1), z(7,1))
plot(b6(1,1), z(7,1))
plot(t6(1,1), z(8,1))
plot(b7(1,1), z(8,1))
plot(t7(1,1), z(9,1))
plot(b8(1,1), z(9,1))
plot(t8(1,1), z(10,1))
plot(b9(1,1), z(10,1))
plot(t9(1,1), z(11,1))
plot(b10(1,1), z(11,1))
plot(t10(1,1), z(12,1))
plot(b11(1,1), z(12,1))
plot(t11(1,1), z(13,1))
plot(b12(1,1), z(13,1))
plot(t12(1,1), z(14,1))
plot(b13(1,1), z(14,1))
plot(t13(1,1), z(15,1))
plot(b14(1,1), z(15,1))
plot(t14(1,1), z(16,1))
plot(b15(1,1), z(16,1))
plot(t15(1,1), z(17,1))
hold off
keyboard

[Clayton Gotberg]

Check what line/marker style you are using in those plot commands. MATLAB uses lines by default and these will not be visible if only a single point is plotted in your command. Check out the LineSpec section of the plot command documentation for more detail.

[Eddy Ramirez]

i tried to add lines and as '-k' and i still dont get anything

[Clayton Gotberg]

The fact that MATLAB is using lines is the issue, as I said before. Use points or asterisks if you want to plot a cloud of points, or send a series of points to plot if you want to use lines.

[Eddy Ramirez]

I tried asterisks as well though and it doesnt work either

[Clayton Gotberg]

I've tried it myself and I can see the chart with the point marker:

plot(x_point,y_point,'k.')

[Eddy Ramirez]

hmm I was using '.k" no wonder is there a way to connect the dots?

I am sorry for my little ignorance with the coding mechanism my goal is to achieve something like the attached picture (top right)

[Eddy Ramirez]

should I run a stairs instead of a plot?

[Clayton Gotberg]

To connect the dots, just plot all of the points at once in the order you want them to be connected, using a line.

If I have three points (x_a, y_a), (x_b,y_b), and (x_c,y_c), I can plot them separately or together:

%Separately
plot(x_a, y_a,'k.')
plot(x_b, y_b,'k.')
plot(x_c, y_c,'k.')

%Together
plot([x_a x_b x_c],[y_a y_b y_c],'k-')
% Or
X = [x_a x_b x_c];
Y = [y_a y_b y_c];
plot(X,Y,'k-')

Answer 1

[LO]

short answer: you are plotting points singularly so you can't see them

short solution: pool them together in arrays then use the command to plot

arrays can be made using the line

[a,b,c,d,.....,z]

where letters are your values for bottom or top

I am not familiar with the sym data format but I think the issue might be just in the way you plot there. Also comment the keyboard otherwise it gets stuck in debugging.

see your code with scatter instead of plot (scatter can plot single points, plot can too but if you plot a single points one by one you will not see them, they have to be organized in arrays)

close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio
%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31
%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter
%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));
%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));
% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
    Q12,Q22,0;
    0,0,Q66];
%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];
%%Laminate rotations in radians
% theta=pi/180*[0;45;-45;90;30;60;-30;-60]';%%DEGREES TO RADIANS
% theta=[theta;theta(end:-1:1)];
z=(-length(theta)/2+(0:length(theta)))'*h;
% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(theta)
    m=cosd(theta(i));
    n=sind(theta(i));
    T=[m.^2,n.^2,2.*m.*n;
        n.^2,m.^2,-2.*m.*n;
        -m.*n,m.*n,m.^2-n.^2];
      
    % Global transformed reduced stiffness coefficients Q_bar
    Q_bar{i}=inv(T).*Q.*T;
end
Aij=0;
Bij=0;
Dij=0;
for j=1:length(theta)
    A_j=Q_bar{j}*(z(j+1)-z(j));
    Aj{j}=A_j;
    Aij=Aij+Aj{j};
    B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
    Bj{j}=B_j;
    Bij=Bij+Bj{j};
    D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
    Dj{j}=D_j;
    Dij=Dij+Dj{j};
end
A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;
%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];
%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];
%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];
%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];
%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];
%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];
test=Q_bar{1}*shearf+z(1,1)*Q_bar{1}*bendingf;
test1=Q_bar{1}.*shearf+z(2,1)*Q_bar{1}.*bendingf;
% scatter(test(1,1), z(1,1));
for k=1:length(theta)
 stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf; %%BOTTOM LAYER
 stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf; %%TOP LAYER
 stress1f{k}=stress1; %Bottom
 stress2f{k}=stress2; %Top
end

b0=stress1f{1};
b1=stress1f{2};
b2=stress1f{3};
b3=stress1f{4};
b4=stress1f{5};
b5=stress1f{6};
b6=stress1f{1};
b7=stress1f{7};
b8=stress1f{8};
b9=stress1f{9};
b10=stress1f{10};
b11=stress1f{11};
b12=stress1f{12};
b13=stress1f{13};
b14=stress1f{14};
b15=stress1f{15};
t0=stress2f{1};
t1=stress2f{2};
t2=stress2f{3};
t3=stress2f{4};
t4=stress2f{5};
t5=stress2f{6};
t6=stress2f{1};
t7=stress2f{7};
t8=stress2f{8};
t9=stress2f{9};
t10=stress2f{10};
t11=stress2f{11};
t12=stress2f{12};
t13=stress2f{13};
t14=stress2f{14};
t15=stress2f{15};
scatter(b0(1,1), z(1,1))
% keyboard
hold on
scatter(t0(1,1), z(2,1))
scatter(b1(1,1), z(2,1))
scatter(t1(1,1), z(3,1))
scatter(b2(1,1), z(3,1))
scatter(t2(1,1), z(4,1))
scatter(b3(1,1), z(4,1))
scatter(t3(1,1), z(5,1))
scatter(b4(1,1), z(5,1))
scatter(t4(1,1), z(6,1))
scatter(b5(1,1), z(6,1))
scatter(t5(1,1), z(7,1))
scatter(b6(1,1), z(7,1))
scatter(t6(1,1), z(8,1))
scatter(b7(1,1), z(8,1))
scatter(t7(1,1), z(9,1))
scatter(b8(1,1), z(9,1))
scatter(t8(1,1), z(10,1))
scatter(b9(1,1), z(10,1))
scatter(t9(1,1), z(11,1))
scatter(b10(1,1), z(11,1))
scatter(t10(1,1), z(12,1))
scatter(b11(1,1), z(12,1))
scatter(t11(1,1), z(13,1))
scatter(b12(1,1), z(13,1))
scatter(t12(1,1), z(14,1))
scatter(b13(1,1), z(14,1))
scatter(t13(1,1), z(15,1))
scatter(b14(1,1), z(15,1))
scatter(t14(1,1), z(16,1))
scatter(b15(1,1), z(16,1))
scatter(t15(1,1), z(17,1))
hold off
% keyboard

Answer 2

[per isakson]

I added 'd' to all your plot-commands

plot(b0(1,1), z(1,1), 'd' )
keyboard
hold on
plot(t0(1,1), z(2,1), 'd' )
plot(b1(1,1), z(2,1), 'd' )
plot(t1(1,1), z(3,1), 'd' )
plot(b2(1,1), z(3,1), 'd' )
etc.

Now the script outputs

[LO]

Using a series of plotting commands, as done, is highly inefficient.

Organizing data in arrays and then plotting once, is recommended.

[Eddy Ramirez]

thank you for the response, I was able to get the dots but I am looking for a stairs plot scenario (pic attached)

I tried to run something like this, but it is not even close

figure(1)
stairs(b0(1,1), z(1,1), 'LineWidth',2,'Marker','d','MarkerFaceColor','c')
hold on
stairs(t0(1,1), z(2,1),'LineWidth',2,'Marker','d','MarkerFaceColor','c')
stairs(b1(1,1), z(2,1), 'LineWidth',2,'Marker','d','MarkerFaceColor','c')
stairs(t1(1,1), z(3,1), 'LineWidth',2,'Marker','d','MarkerFaceColor','c')
hold off

[per isakson]

To draw a line you need a vector of at least two elements. That's true for both plot and stairs. I joined only a subset of your points.

x = [t0(1,1),b1(1,1),t1(1,1),b2(1,1),t2(1,1),b3(1,1),t3(1,1),b4(1,1),t4(1,1),b5(1,1),t5(1,1),b6(1,1),t6(1,1),b7(1,1),t7(1,1)];
y = [z(2,1),z(2,1),z(3,1),z(3,1),z(4,1),z(4,1),z(5,1),z(5,1),z(6,1),z(6,1),z(7,1),z(7,1),z(8,1),z(8,1),z(9,1)];
stairs( x, y, '-dk' )

outputs

and replace '-dk' by '-k' to remove the markers

Not all of the lines of unnamed.png are vertical or horizontal. plot( x, y, '-dk' ) produces