Plotting a 3D mode shape by revolving 2D datapoints (in the r-z-plane) of an axisymmetric geometry around the z-axis

Dear Matlab community,

I would like to generate a 3D graph of a vibrational mode shape using 2D data of a solid of revolution (the cross-section of which lies in the r-z plane of a cylindrical coordinate system), by revolving the 2D data around the z-axis in the circumferential direction θ by 2π radians. I have a point cloud of data points of the meshed geometry where each corresponds to a node in a finite element model. From that model I also have the the three displacement components at each node, u (in radial direction r), v (in circumferential direction θ) and w(in axial direction z).

The undeformed geometry can be seen in the following graph (which is a simplification of the actual geometry I am solving for, but essentially includes similar geometric features, but us not as symmetric)

The red dots correspond to the nodes in the finite element model and the blue line corresponds to the extreme edge(s).

The displaced geometry (,,) of a mode at the plane where θ=0 looks as follows (after scaling it):

where the cyan line shows the edges of the undeformed geoemetry. I neglected the circumferential displacement here but it needs to be included in the 3D figure.

The displacements vary sinusoidally around the circumference, which needs to be included in the final 3D plot as:

Ideally, I would like to plot the 3D mode shape as shown in the figure below. It is taken directly from a finite element software, but I can't use it, since I apply further processing outside the software. But after finishing my calculations, I again end up with 3 displacement components at the same nodes and could plot the new ones. The colourmap c represents the magnitude of the overall displacement at each point as .

In the appendix I provide the data of the wheel in the data.mat file. The file is required to run the code below. It includes the nodal coordinates r and z of the above-shown undeformed 2D geometry and the three displacement components for the single mode. It would be sufficient to revolve the outer edges (blue lines in the first/second figure), instead of revolving all points. A good starting point is probably to try and plot the undeformed shape first, before going any further. Afterwards, the displacements could be added to each point and maybe it is possible to add a colourmap. I add a small piece of code below, which only includes my start. I have, however, already had problems coming up with something to plot the undeformed geoemtry, other than creating a 3D point cloud. I am very sorry, I can't provide you more. I'm not very experienced in plotting something in 3D in MATLAB and I am not sure how to tackle this problem, to be honest.

I am really hoping anyone could provide me some adivce on how to plot a 3D shape from the 2D data, as shown above. I am grateful for anything, I hope my description makes sense. This is the first time I am using the forum, I hope I was precise enough for you, to understand my problem. If not, please don't hesitate to tell me and I will do my best providing all neccesary information. Thank you very much in advance.

Best wishes,

Christopher

%Variables from data.mat are:
%r: radial coordinate of each node
%z: axial coordinate of each node
%u: displacement in radial direction of each node
%v: displacement in circumferential direction of each node
%w: displacement in axial direction of each node
%Connect: Connectivity matrix (links nodes to individual finite elements)
clc;clear variables

%load data from data.mat
load('data.mat')

%angle theta around circumference in steps of pi/N
N = 30;
theta = 0:pi/N:2*pi; 
%number of nodal diameters
n = 2; 

%displacement magnitude or colourmap
c = sqrt(u.^2+v.^2+w.^2); 

%re-scale displacements to max displacement
u = u./max(abs([u', v', w']));
v = v./max(abs([u', v', w']));
w = w./max(abs([u', v', w']));

%r,z,u,v,w on boundary or extreme edge nodes
rEdg = r(pEdge);
zEdg = z(pEdge);
uEdg = u(pEdge);
vEdg = v(pEdge);
wEdg = w(pEdge);

%rescale mode shape in plotby scale factor to look proper 
scale = 0.1 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%   plots   %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% create 3D point cloud of undeformed geometry for each node/value of theta
figure(1)
clf

%create arrays for each value of theta
%all nodes
R0 = r.*ones(1,length(theta));
Z0 = z.*ones(1,length(theta));
%edge nodes (comment out if all nodes shall be used)
% R0 = rEdg.*ones(1,length(theta))
% Z0 = zEdg.*ones(1,length(theta))

Th0 = theta.*ones(length(R0),1);

%transform cyl polar to cartesian coords
[X0,Y0,Z0] = pol2cart(Th0,R0,Z0);

plot3(X0,Y0,Z0,'k.')
axis equal
axis padded

%% create 3D point cloud of deformed geometry for each node/value of theta
figure(2)
clf

%create arrays for each value of theta
%all nodes
% R = r+u*scale.*ones(1,length(theta)).*cos(n*theta);
% Z = z+w*scale.*ones(1,length(theta)).*cos(n*theta);
% Th = theta+v*scale.*sin(n*theta);

%edge nodes (comment out if all nodes shall used)
R = rEdg+uEdg*scale.*ones(1,length(theta)).*cos(n*theta);
Z = zEdg+wEdg*scale.*ones(1,length(theta)).*cos(n*theta);
Th = theta+vEdg*scale.*sin(n*theta);

%transform cyl polar to cartesian coords
[X,Y,Z] = pol2cart(Th,R,Z);

plot3(X,Y,Z,'k.')
axis equal
axis padded

%% plot 3D as slices at angle theta increment

figure(3)
clf

for ith = 1:length(theta)
u_theta = u*scale.*cos(n*theta(ith));
v_theta = v*scale.*sin(n*theta(ith));
w_theta = w*scale.*cos(n*theta(ith));

for i = 1:size(Connect,1)
    Conn = (Connect(i,:));
    re = r(Conn)+u_theta(Conn);
    ze = z(Conn)+w_theta(Conn);
    the = theta(ith)+v_theta(Conn); 
    [Xe,Ye,Ze] = pol2cart(the,re,ze);
    patch(Xe,Ze,Ye,[0.5 0.5 0.5])
    hold on
end
view(-45,45)
axis equal
axis([-0.6 0.6 -0.1 0.1 -0.6 0.6])   
 
end

Answer 1

[Matt J]

I can't follow most of what you posted, but if you just want to generate a surface/solid of revolution, you can use cylinder()

https://www.mathworks.com/help/matlab/ref/cylinder.html

[Christopher Knuth]

Dear Matt, thank you for your answer. I am sorry if I was not very clear, I will try reformulating the question later and edit it. A short description follows below. I looked at the cylinder function already, but haven't figuried out yet, how I could use it for my case. I want to generate a surface of revolution, but the cross-section is varying with the angle or revolution, while revolving around one axis.

Essentially, I have data in a 2D grid (x,y), which is shown as red dots in the second figure of my original post. The coordinates of each point are made of a spatial coordinate (r and z) of the undeformed solid (as in the first figure of my origial post), plus a small displacement (u and w). Consequently, the coordinates of the 2D plot are and . This data represents the cross-section of a 3D structure which can be generated by revolution of the 2D data points around the y-axis. Additionally, the displacements vary sinusoidally around the circumference of the full 3D solid. For my purpose, it would be sufficient to only rotate the red dots which lie on the outer edge (conicide with blue line), since I only want to show the deformed boundary surface of the solid, as seen in the third figure of the original post. I can generate the cross-sectional data at each angle of θ between 0 and 2π but I don't know, if it is possible using that data to generate a 3D plot that looks like the third figure, i.e. one that shows the "enveloping surface" only. To each data point I can assign a value which could generate the colourmap. I hope that helps clarifying. Best wishes, Christopher

PS: I tried to put it into one figure, see below. Revolving the red data points on the outer edge of the solid in the x-y plane around the y axis. I show the slices at steps of to see how the cross-section varies with θ. The blue curve in 2D would generate the "enveloping" surface in 3D. I hope this is helpful.

[Christopher Knuth]

Thanks for trying to help me @Matt J. After some thoughts I was able to resolve the problem on my own. I did not use cylinder() as you suggested, I wasn't quite clear how to use in my case. But I found another way, which I describe in an answer below, with added code.

I never used the forum, can/shall I close the question by accepting my own answer? That seems slightly odd to me.

Best Chris

Answer 2

[Christopher Knuth]

I have managed to plot the 3D figure from my 2D data by using the patch function. First, I create a 3D boundary mesh from my 2D data (only blue line) by revolving the points (red dots) in the 2D plane step-by-step by an angle increment. I preallocate an array with the corner points of each boundary element. Then each element can be plotted with the patch() function.

The raw 2D data must look like:

[in my case each point has a radial and axial coordinate in (r,z) + sinusoidally varying mechanical displacements (u,v,w) in 3D in (r,θ,z)]

The resulting figure for above-shown arbitrary mode of a disc is depicted below.

My code looks like the following, in case someone is interested in doing something similar. The data is specific to finite element vibration analysis, but the code may be adopted for similar problems for plotting a figure of a solid of revolution. It can be used for any type of cross-section, if one has the connectivity matrix that links the nodes (red dots) to the edge elements (blue lines), though I am showing a disc which might be easier implemented, as well. It might not be most efficient, but it works and does the job.

Still, if anyone knows a function which does that job directly or another way, I am happy for suggestions to improvements, although with my current code speed is not really a problem for the number elements/angles I require.

clc;clear variables;

%% numerical example data

%radial coordinate (undeformed) (all nodes in 2D plane)
r = [0.10,0.10,0.20,0.10,0.20,0.30,0.20,0.30,0.40,0.30,0.40,0.50,0.40,0.50,0.60,0.50,0.60,0.60]';
%axial coordinate (undeformed) (all nodes in 2D plane)
z = [-0.050,0,-0.050,0.050,4.16e-18,-0.050,0.050,0,-0.050,0.050,0,-0.050,0.050,-1.04e-18,-0.050,0.050,0,0.050]';
%circumferential coordinate (undeformed) for extension to 3D
Nangle = 40;
th = linspace(0,2*pi,Nangle);
%displacement in radial direction (all nodes in 2D plane)
u = [0;0;0.0140;0;0;0.00300;-0.0140;0;-0.0170;-0.00300;0;-0.0300;0.0170;0;-0.0330;0.0300;0;0.0330];
%displacement in circumferential direction (all nodes in 2D plane)
v = [0;0;-0.020;0;0;-0.019;0.020;0;-0.010;0.019;0;0.0020;0.010;0;0.0080;-0.0020;0;-0.0080];
%displacement in axial direction (all nodes in 2D plane)
w = [0;0;0.0420;0;0.0430;0.0700;0.0420;0.0730;0.0540;0.0700;0.0570;0;0.0540;0.00100;-0.0680;0;-0.0670;-0.0680];
%Connectivity matrix that links node numbers to edge elements (not the element inside domain)
Connect = [2,1;1,3;4,2;3,6;7,4;6,9;10,7;9,12;13,10;12,15;16,13;15,17;18,16;17,18];

%scaling factor for displacements in illustration
scale = 0.1;
%number of nodal diameters of circumferential variation (zero lines)
n = 1;

%maximum displacement (for scaling)
Dmax = max(abs([u',v',w']));
%max displacement magnitude (for colourmap)
Cmax = max(sqrt(u.^2+v.^2+w.^2));

%number of elements and nodes/element
NElem = size(Connect,1);
NNode = size(Connect,2);

%% preallocate matrix including node coordinates of each patch
for iElem = 1:NElem % sweep over each edge element

    %connectivity for i-th element
    Conn = Connect(iElem,:);
    
    %elemental coordinates (re,ze) and displacements (ue,ve,we)
    re = r(Conn);
    ze = z(Conn);
    ue = u(Conn);
    ve = v(Conn);
    we = w(Conn);
    
    %rotate element in circle around circumferential coordinate
    for ith = 1:Nangle-1 %sweep over each angle increment(except last because is at 2pi or 0)
        
        %sinusiodally varying displacements of the nodes at angle i and angle i+1
        ueth = [ue*cos(n*th(ith));flipud(ue)*cos(n*th(ith+1))];
        veth = [ve*sin(n*th(ith));flipud(ve)*sin(n*th(ith+1))];
        weth = [we*cos(n*th(ith));flipud(we)*cos(n*th(ith+1))];
            
        %colour map/nodal displacement magnitude for boundary element and at each angle
        ceth(:,ith) = sqrt(ueth.^2+veth.^2+weth.^2);
        %nodal coordinate+displacement(re-scaled) for boundary element and at each angle
        reth(:,ith) = [re;flipud(re)]+ueth/Dmax*scale;
        theth(:,ith) = [th(ith)*ones(NNode,1);th(ith+1)*ones(NNode,1)]+veth/Dmax*scale;
        zeth(:,ith) = [ze;flipud(ze)]+weth/Dmax*scale;
        
    end
   
    %transform cylindrical to cartesian coordinates
    [Xe,Ye,Ze] = pol2cart(theth,reth,zeth);
    %store solution for each boundary element and each angle in 3D array
    Ctmp(:,iElem,1:size(ceth,2)) = ceth;
    Xtmp(:,iElem,1:size(ceth,2)) = Xe;
    Ytmp(:,iElem,1:size(ceth,2)) = Ye;
    Ztmp(:,iElem,1:size(ceth,2)) = Ze;
    
end

%reshape 3D arrays to 2D (row: nodes of patch ,column: patches)
C = reshape(Ctmp,2*NNode,[]);
X = reshape(Xtmp,2*NNode,[]);
Y = reshape(Ytmp,2*NNode,[]);
Z = reshape(Ztmp,2*NNode,[]);

%% plot 3D figure

figure(1)
clf
patch(X,Y,Z,C,'FaceColor','interp');
% patch(X,Y,Z,C,'FaceColor','interp','EdgeColor','none');
colorbar
colormap jet
caxis([0,Cmax])
axis equal
view(0,45)
xlim([-(max(r)+scale) max(r)+scale])
ylim([-(max(r)+scale) max(r)+scale])
zlim([-(max(z)+scale) max(z)+scale])
axis off

[ELizabeth Sharp]

Hi, Would this be adaptable for plotting the deformation on a cylinder for 200 degrees of rotation?

I am a bit confused with the rotate element section.

Thanks!

[Christopher Knuth]

Hi, Elizabeth,

yes it should be adaptable, as long as it is an axisymmetric structure modelled in a 2D plane. I defined the rotation for angles "th = linspace(0,2*pi,Nangle);", which you could change from 2pi to any angle. But if you do not rotate a full circle, you will have an open section and may need to find a workaround to close it, if that's what you want.

The "rotate element" section basically creates a mesh of the boundary of the structure by sweeping sectionwise over all angles defined in th using the sin/cosine distribution of the displacements as defined in finite element literature.

I hope this makes sense and helps.

Best, Christopher

[Romain Liechti]

Thank you so much! I had vibrometer measurements for a circular device, but only along a single radius. I was able to adapt the code to visualize the axisymmetric device's modal shape in 3D. It's fantastic!