3-D plot with 3 variables resulting in 3 responses

Hello Everyone,

I am trying to create a 3-d plot with 3 independent variables that result in 3 responses in the three orthogonal directions of the structure. EF= scaling factor that I used amplifies the earthquake effect in the dynamic analysis. I want to see a 3-d plot concerning the change of EF in three orthogonal directions. Any help will be appreciated. Here is the code that I simplified and attached below.

Best,

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('Earthquake_Elcenctro_data.mat') % loading Elcentro data in g 
load('repsonses.mat')
Qx=Elctorsion(2,:)                  % rotation in radyan about x-direction

Qx = 1×2501

1.0e+00 * -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0001 -0.0001 -0.0000 -0.0000 -0.0000 0.0001 0.0001 0.0000 -0.0000 -0.0001 -0.0002 -0.0001 -0.0000 0.0001 0.0000 0.0000 0.0000 -0.0000

Qy=Qx;                              % rotation in radyan about y-direction
Qz=Qx;                              % rotation in radyan about z-direction
T=zeros(4,4,5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ux=max(abs(yout(:,9)));
Uy=max(abs(yout(:,18)));
Uz=0;
%%%% building geometric and earthquake data inputs
Tx=max(abs(Elcxyonu(2,:)));   % earthquake data x-translational component in g
Ty=max(abs(Elcyyonu(2,:)));   % earthquake data y-translational component in g
Tz=max(abs(Elczyonu(2,:)));   % earthquake data z-translational component in g
dx=max(abs(Elctorsion(2,:))); % earthquake data x-rotational component in radyan
dy=dx;                        % earthquake data y-rotational component in radyan
dz=dx;                        % earthquake data z-rotational component in radyan
Umax=[Ux Uy Uz]';             % top floor displacement vector over time without rotaional components of earthquake.
for i=1:5
    
    %%%% rotation matrix with Eular angles 
    R=[cos(Qy(i))*cos(Qz(i)),                                -cos(Qy(i))*sin(Qz(i)),                        sin(Qy(i));
        (sin(Qx(i))*sin(Qy(i))*cos(Qz(i))+cos(Qx(i)*sin(Qz(i)))),   sin(Qx(i))*sin(Qy(i))*sin(Qz(i))+cos(Qx(i))*cos(Qz(i)),   sin(Qx(i))*sin(Qy(i));
        sin(Qx(i))*sin(Qz(i))-cos(Qx(i))*sin(Qy(i))*cos(Qz(i)),   sin(Qx(i))*cos(Qz(i))-cos(Qx(i))*sin(Qy(i))*sin(Qz(i)) ,  cos(Qx(i))*cos(Qy(i))];
    %%%% Transformation matrix  
    T(:,:,i)=[R            Umax;
        zeros(1,3) ones(1,1)];
end
Ef=zeros(1,4);
t=1;
xyz=zeros(90405:4);
UT_max=zeros(4,1,5);
%%%% Final total translation response including rotational components
for i=0:0.5:10
    for j=0:0.5:10
        for k=0:0.5:20
            for z=1:5
                Ef=[i,j,k  ones(1,1)]';
                UT_max(:,:,t)=T(:,:,z)*Ef;
                xyz(t,:)=Ef;
                t=t+1;
            end
        end
    end
end
N=length(UT_max)

N = 90405

result=zeros(N,1);
for i=1:N
    resultx(i)=UT_max(1,1,i);
    resulty(i)=UT_max(2,1,i);
    resultz(i)=UT_max(3,1,i);
end
result=[resultx; resulty;resultz]';
figure(1)
scatter3(xyz(:,1),xyz(:,2),xyz(:,3),result);

Error using scatter3
Value must be a vector of positive numeric type or nan.

xlabel('Amplification factor (Ef_x)');
ylabel('Amplification factor (Ef_y)') 
zlabel('Amplification factor (Ef_z)')
figure(2)
pdeplot3D(xyz(:,1),xyz(:,2),xyz(:,3),"ColorMapData",result);
xlabel('Amplification factor (Ef_x)');
ylabel('Amplification factor (Ef_y)') 
zlabel('Amplification factor (Ef_z)')

Answer 1

[Cris LaPierre]

Some of your values in your variable result are <=0. They must be positive (can't have a marker with a size <=0)

[Osman AKYUREK]

Hi Cris,

Thank you for your quick reply. You are right. I can correct it by using the 'abs' command, but indeed this graph is not what I want. Do you know how to plot 3-d by interpolating with undefined variable intervals? I also want to see the color maps beside the graph. For example, EFx is changing between 0 to 10 with one incremental. What about 0.4 between 0 and 1 when EFx is employed? Can I use the "mesh grid" command in my case? if it is possible, how?

Best regards,

[Cris LaPierre]

You will need to do more than just abs, as some of your values in result are 0, which will result in the same error.

With that correction, right now your code is creating 3 markers at each point. I'm not sure how helpful that is, as I can't discern between them. Perhaps you could better describe what it is you are trying to see? Is there an example you can share?

You may be insterested in scatteredinterpolant or griddedinterpolant, depending what it is you are trying to do.

[Osman AKYUREK]

Hi Cris,

Let me try to explain what I want to do, EFx, EFy, and EFz are scaling (amplification) factors. For example, when Efx=1, Efy=2,Efz=5, the structure, which is exposed to a torsional earthquake with these scaling factors, gives responses in the x- direction, y-direction, and z direction. Therefore, EFx, Efy, Efz take values between 0 to 10 with arbitrary, they result in 3 times each combination of EFx, Efy, Efz. Actually, I want to have mesh grid graph in 3-d format but I couldn't. I hope this clarifies my case then you can help.

Thank you for your interest in this issue.

[Cris LaPierre]

A meshgrid creates a surface - one point in Z for every X,Y combination, while scatter3 creates a cube of points, one at every X, Y and Z coordinate. When you say you want meshgrid, does that mean you want a surface? If yes, how would I know which Z value to select for each X-Y pair?

My current confusion is this. If xyz contains your amplification factors (according to the labels of the scatter3 plot), then Efx=1, Efy=2,Efz=5 corresponds to a single point. So which variable contains the data you want to see as a meshgrid? What data is multiplied by these amplification factors to see a torsional earthquake?

[Osman AKYUREK]

Hi Criss,

You are right. I changed my code. There are only two direction I am considering right now which are EFx and EFy, and EFz is neglected. The results separately resultx and resulty depends on EFx and EFy. We can now mesh-grid (EFx, EFy) and get the surf plot according to your explanation. Thanks in advance. The updated code is below;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('Earthquake_Elcenctro_data.mat') % loading Elcentro data in g 
load('repsonses.mat')
Qx=Elctorsion(2,:)                  % rotation in radyan about x-direction
Qy=Qx;                              % rotation in radyan about y-direction
Qz=Qx;                              % rotation in radyan about z-direction
T=zeros(4,4,5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ux=max(abs(yout(:,9)));
Uy=max(abs(yout(:,18)));
Uz=0;
%%%% building geometric and earthquake data inputs
Tx=max(abs(Elcxyonu(2,:)));   % earthquake data x-translational component in g
Ty=max(abs(Elcyyonu(2,:)));   % earthquake data y-translational component in g
Tz=max(abs(Elczyonu(2,:)));   % earthquake data z-translational component in g
dx=max(abs(Elctorsion(2,:))); % earthquake data x-rotational component in radyan
dy=dx;                        % earthquake data y-rotational component in radyan
dz=dx;                        % earthquake data z-rotational component in radyan
Umax=[Ux Uy Uz]';             % top floor displacement vector over time without rotaional components of earthquake.
for i=1:5
    
    %%%% rotation matrix with Eular angles 
    R=[cos(Qy(i))*cos(Qz(i)),                                -cos(Qy(i))*sin(Qz(i)),                        sin(Qy(i));
        (sin(Qx(i))*sin(Qy(i))*cos(Qz(i))+cos(Qx(i)*sin(Qz(i)))),   sin(Qx(i))*sin(Qy(i))*sin(Qz(i))+cos(Qx(i))*cos(Qz(i)),   sin(Qx(i))*sin(Qy(i));
        sin(Qx(i))*sin(Qz(i))-cos(Qx(i))*sin(Qy(i))*cos(Qz(i)),   sin(Qx(i))*cos(Qz(i))-cos(Qx(i))*sin(Qy(i))*sin(Qz(i)) ,  cos(Qx(i))*cos(Qy(i))];
    %%%% Transformation matrix  
    T(:,:,i)=[R            Umax;
        zeros(1,3) ones(1,1)];
end
Ef=zeros(1,4);
t=1;
xyz=zeros(90405:4);
UT_max=zeros(4,1,5);
%%%% Final total translation response including rotational components
for i=0:0.5:10
    for j=0:0.5:10
        %         for k=0:0.5:20
        k=0
        for z=1:5
            Ef=[i,j,k  ones(1,1)]';
            UT_max(:,:,t)=T(:,:,z)*Ef;
            xyz(t,:)=Ef;
            t=t+1;
        end
        %         end
    end
end
N=length(UT_max)
result=zeros(N,1);
for i=1:N
    resultx(i)=abs(UT_max(1,1,i));
    resulty(i)=abs(UT_max(2,1,i));
    resultz(i)=abs(UT_max(3,1,i));
end
result=[resultx; resulty;resultz]';
figure(1)
surf(xyz(:,1),xyz(:,2),resultx);
xlabel('Amplification factor (Ef_x)');
ylabel('Amplification factor (Ef_y)') 
zlabel('Result in the x-direction')
figure(2)
pdeplot3D(xyz(:,1),xyz(:,2),resultx,"ColorMapData",result);
xlabel('Amplification factor (Ef_x)');
ylabel('Amplification factor (Ef_y)') 
zlabel('Result in the x-direction')

[Cris LaPierre]

In this case, I would just use reshape to turn the X, Y, and Z data into matrices

X = reshape(xyz(:,1),105,[]);

Y = reshape(xyz(:,2),105,[]);

Z = reshape(resultx,105,[]);

The 105 comes from the fact that the first 105 rows of xyz are all for the same X value, so I assume that is one of the dimensions of your grid. The '[]' allows reshape to figure out what the second dimension is automatically.

I'm not sure what you are trying to do with pdeplot, but you do not appear to be using it correctly (for example, you have not created nor passed in a PDEModel object), so I have removed it.

load('Earthquake_Elcenctro_data.mat') % loading Elcentro data in g 
load('repsonses.mat')
Qx=Elctorsion(2,:)                  % rotation in radyan about x-direction

Qx = 1×2501

1.0e+00 * -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0001 -0.0001 -0.0000 -0.0000 -0.0000 0.0001 0.0001 0.0000 -0.0000 -0.0001 -0.0002 -0.0001 -0.0000 0.0001 0.0000 0.0000 0.0000 -0.0000

Qy=Qx;                              % rotation in radyan about y-direction
Qz=Qx;                              % rotation in radyan about z-direction
T=zeros(4,4,5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ux=max(abs(yout(:,9)));
Uy=max(abs(yout(:,18)));
Uz=0;
%%%% building geometric and earthquake data inputs
Tx=max(abs(Elcxyonu(2,:)));   % earthquake data x-translational component in g
Ty=max(abs(Elcyyonu(2,:)));   % earthquake data y-translational component in g
Tz=max(abs(Elczyonu(2,:)));   % earthquake data z-translational component in g
dx=max(abs(Elctorsion(2,:))); % earthquake data x-rotational component in radyan
dy=dx;                        % earthquake data y-rotational component in radyan
dz=dx;                        % earthquake data z-rotational component in radyan
Umax=[Ux Uy Uz]';             % top floor displacement vector over time without rotaional components of earthquake.
for i=1:5
    
    %%%% rotation matrix with Eular angles 
    R=[cos(Qy(i))*cos(Qz(i)),                                -cos(Qy(i))*sin(Qz(i)),                        sin(Qy(i));
        (sin(Qx(i))*sin(Qy(i))*cos(Qz(i))+cos(Qx(i)*sin(Qz(i)))),   sin(Qx(i))*sin(Qy(i))*sin(Qz(i))+cos(Qx(i))*cos(Qz(i)),   sin(Qx(i))*sin(Qy(i));
        sin(Qx(i))*sin(Qz(i))-cos(Qx(i))*sin(Qy(i))*cos(Qz(i)),   sin(Qx(i))*cos(Qz(i))-cos(Qx(i))*sin(Qy(i))*sin(Qz(i)) ,  cos(Qx(i))*cos(Qy(i))];
    %%%% Transformation matrix  
    T(:,:,i)=[R            Umax;
        zeros(1,3) ones(1,1)];
end
Ef=zeros(1,4);
t=1;
xyz=zeros(90405:4);
UT_max=zeros(4,1,5);
%%%% Final total translation response including rotational components
for i=0:0.5:10
    for j=0:0.5:10
        %         for k=0:0.5:20
        k=0;
        for z=1:5
            Ef=[i,j,k  ones(1,1)]';
            UT_max(:,:,t)=T(:,:,z)*Ef;
            xyz(t,:)=Ef;
            t=t+1;
        end
        %         end
    end
end
N=length(UT_max)

N = 2205

result=zeros(N,1);
for i=1:N
    resultx(i)=abs(UT_max(1,1,i));
    resulty(i)=abs(UT_max(2,1,i));
    resultz(i)=abs(UT_max(3,1,i));
end
result=[resultx; resulty;resultz]';

X = reshape(xyz(:,1),105,[]);
Y = reshape(xyz(:,2),105,[]);
Z = reshape(resultx,105,[]);
figure(1)
surf(X,Y,Z);
xlabel('Amplification factor (Ef_x)');
ylabel('Amplification factor (Ef_y)') 
zlabel('Result in the x-direction')