# How to plot contour on a curved surface for the below problem?

13 views (last 30 days)
Junaid on 11 Sep 2023
Edited: Junaid on 29 Sep 2023
Hello,
Can anyone here help me with the plotting contours on a curved surfaceThe code should be able to generate the plot as in ooutput image.
##### 3 CommentsShow 1 older commentHide 1 older comment
Nathan Hardenberg on 11 Sep 2023
Yeah, works as expected! But your "surface_x_value" function just does not give correct results. Since I don't know how you calculate those I can't help you further. I see that you do an optimization, but somewhere lies an error.
Two small things I found:
1. in "distance_from_pith" you use y^2. But you probably need y.^2
2. Your inital guess seems to be very weird. Your distance to pith should not be a good guess for your x-value. Maybe a better guess would be the surface of the cylinder.
As you can see below, if I draw the red circles on the cylinder it works. (except for the imaginary parts, since they are not on the cylinder)
... % Only changed the surface_x_value function (see below)
function x_new = surface_x_value(y, z, k, zc0, Ri0)
Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored
Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored
x_new = sqrt(Ri0^2-y.^2); % calculate points on cylinder (Ri0 = radius of cylinder)
end
Nathan Hardenberg on 14 Sep 2023
Here is all the code that is needed to function. Again only your code except the function at the end.
clear all
close all
z_root=0;
z_top=150;
dk_log=40; % extension of drawn knot, i.e. length in addition to Ri0
Ri0=80;
% Make a drawing to illustrate a part of a log (two circles in 3D space representing the circumference of the log
% and dashed lines representing the pith and a plane through the log.)
figure(1)
hold on
teta_mesh=[0:5:360]/180*2*pi;
plot3([-Ri0 Ri0 Ri0 -Ri0 -Ri0]',[0 0 0 0 0]',[z_root z_root z_top z_top z_root],':k')
plot3([Ri0*cos(teta_mesh)]',[Ri0*sin(teta_mesh)]',[teta_mesh*0+z_root]','-k')
plot3([Ri0*cos(teta_mesh)]',[Ri0*sin(teta_mesh)]',[teta_mesh*0+z_top]','-k')
plot3([0 0]',[0 0]',[z_root z_top],':r')
box on
axis equal
view(3)
% Define parameters for the geometry of a knot according to Foley's
% model/thesis.
zc0=40; % z-position of the knot apex at the pith
Ax100=41.5; %page 88
L100=28.6; %page 88
Dl=1.22e-2*L100+0.1252; %page 91
Cax=9.7e-3*Ax100+0.1725; %page 91
kappa=0.95; %page 93
c1=-1.458e-3; %page 91
c2=0.01*Ax100+0.1458; %page 92,
x_mesh=[0:4:(Ri0+dk_log)];
lx_mesh=-1e-7*x_mesh.^4+3e-5*x_mesh.^3-4e-3*x_mesh.^2+Dl*x_mesh; %page 91
bx_mesh=lx_mesh*kappa; %page 93
zc_mesh=c1*x_mesh.^2+Cax*x_mesh; %page 92
k=1/kappa; %page 94, eq. 5.6
% Add to the drawing lines representing the geometry of a knot
XYZp=zeros(181,3);
XYZp(:,1)=Ri0+dk_log;
zc=c1*(Ri0+dk_log)^2+Cax*(Ri0+dk_log);
plot3(x_mesh,x_mesh*0,zc_mesh+zc0,'--k')
plot3(x_mesh,x_mesh*0,zc_mesh+zc0-lx_mesh/2,'-k')
plot3(x_mesh,x_mesh*0,zc_mesh+zc0+lx_mesh/2,'-k')
knot_ovals_x_pos=[30:30:120];
for xov=knot_ovals_x_pos
XYZp=zeros(181,3);
XYZp(:,1)=xov;
pos=0;
zc=c1*xov^2+Cax*xov;
for v=0:2:360
pos=pos+1;
vr=v/180*pi;
ypos=dp*cos(vr)/k;
zpos=dp*sin(vr);
XYZp(pos,2:3)=[ypos zpos];
if pos>1
plot3(XYZp((pos-1):pos,1),XYZp((pos-1):pos,2),XYZp((pos-1):pos,3)+zc0+zc,':k');
end
end
end
% Add to the figure a number of contour curves (red color) in the 'bulged part' of a 'growth surface' with
% Ri0 = 80 mm, see pages 62-67 in Foley's thesis.
Bbump=2.0; %page 93
Abump=(0.0217*L100-0.2); %page 93
Aexp=0.0056*L100+1.96; %page 93
kappa_limit=0.99; %page 94
Rik=80;
%--- Add the missing code here ---
% Define the function A(x) using equation (5.3.a)
A = @(x) -1e-7*x.^4 + 3e-5*x.^3 - 4e-3*x.^2 + (1.22e-2*L100 + 0.1252).*x;
% Define the function zc(x) using equation (5.3.b)
zc = @(x) -1.458e-3*x.^2 + (9.7e-3*Ax100 + 0.1725).*x;
% Plotting the growth surface contours
p_vals = [15:5:75];
for zp = p_vals % page 65, fig 4.7
theta_vals = linspace(0, 2 * pi, 100);
% Calculate the distance from the pith to the contour curve
distance_to_pith = sqrt(Ri0^2 - zp^2);
% Adjust the radius based on the distance to the pith and aspect ratio k
r_vals = distance_to_pith * k; % k is the aspect ratio
y_vals = r_vals .* cos(theta_vals); % Adjusted for YZ plane
z_vals = 35 + zc0 + r_vals .* sin(theta_vals); % Adjusted for YZ plane
% Get the x values for the growth surface
x_new = arrayfun(@(y, z) surface_x_value(y, z, k, zc0, Ri0), y_vals, z_vals);
plot3(x_new, y_vals, z_vals, 'r'); % plotting the growth surface
end
Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored
Warning: Imaginary parts of complex X, Y, and/or Z arguments ignored
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
% Function to calculate the x-coordinate of the growth surface for a given y and z
%-----------------------------------
grid on
axis equal
function x_new = surface_x_value(y, z, k, zc0, Ri0)
x_new = sqrt(Ri0^2-y.^2);
end

J. Alex Lee on 16 Sep 2023
This is how I would do it with the bump.
% parameters
z_root=0;
z_top=150;
Ri0 = 80;
trunk.TRes = 360/5;
trunk.ZRes = 30;
dk_log=40; % extension of drawn knot, i.e. length in addition to Ri0
zc0=40; % z-position of the knot apex at the pith
Ax100=41.5; %page 88
L100=28.6; %page 88
Dl=1.22e-2*L100+0.1252; %page 91
Cax=9.7e-3*Ax100+0.1725; %page 91
c1=-1.458e-3; %page 91
c2=0.01*Ax100+0.1458; %page 92,
kappa = .7; % free choice how oval are the rings
% 1. the symmetrical gaussian bump
bumpH = 10; % free choice approx height of bump relative to Ri0
bumpVar = 200; % free choice how wide/gradual bump is
numRings = 15; % free choice number of rings to draw
ringRes = 200;
% knot geometry
knotCenterFcn = @(r)(polyval([c1,Cax,zc0],r));
knot.centerline.n = 50;
knot.centerline.r = linspace(0,Ri0 + dk_log,knot.centerline.n).';
knot.centerline.theta = 0 * ones(knot.centerline.n,1);
[knot.centerline.x,knot.centerline.y] = pol2cart(knot.centerline.theta,knot.centerline.r);
knot.centerline.z = knotCenterFcn(knot.centerline.r);
% create gridded surface coordinates (easy way)
[trunk.X,trunk.Y,trunk.Z] = cylinder(Ri0*ones(trunk.ZRes,1),trunk.TRes);
trunk.Z = trunk.Z*(z_top-z_root) + z_root;
% warning: this will break for knot.centerline.theta .ne. 0
knot.surf.Y = knot.surf.Y * kappa;
%--- Add the missing code here ---
% 1. model the bump as a gaussian
% 2. generate the contours using built in contourc
% 3. use Adam Danz's contour line extractor
% 4. wrap the contour line cooridnates on the trunk cylinder with trig
%
% steps 2 and 3 can be done manually, but why not use built-in
zPith = knotCenterFcn(Ri0 + bumpH);
% need to calculate the gaussian exponential factor A s.t. at the height of
% the bump for rPothAtTrunk = bumpH
% bumpH = A*exp(... + kappa*y^2/A)
A = bumpH/exp(-(rPithAtTrunk^2)/bumpVar);
% the double gaussian function
bumpFcn = @(z,y)(A*exp(-((z-zPith).^2+(y/kappa).^2)/bumpVar) + Ri0);
z = linspace(z_root,z_top,ringRes);
y = linspace(-Ri0,Ri0,ringRes).';
XFlat = bumpFcn(z,y);
% some logic to space the rings
skewedlevels = linspace(0,1,numRings).^.3;
levelLims = log10([1e-8,bumpH]);
loglevels = skewedlevels*diff(levelLims) + levelLims(1);
levels = Ri0 + 10.^loglevels;
% get the contour coordinates using built-ins plus FEX
[cm] = contourc(z,y,XFlat,levels);
[ct,~] = getContourLineCoordinates(cm);
grps = unique(ct.Group);
for i = numel(grps):-1:1
grp = ct(ct.Group==grps(i),:);
% contour assumed z-axis is height map
% so need to switch some coordinates up
x = grp.Level;
y = grp.Y;
% use trig to wrap around trunk
rings(i).x = x .* cos(y./x);
rings(i).y = x .* sin(y./x);
rings(i).z = grp.X;
end
% draw
figure(1);
tl = tiledlayout(2,2);
ax(1) = nexttile([1,2]);
ax(2) = nexttile(tl);
ax(3) = nexttile(tl);
mesh(ax(1),trunk.X,trunk.Y,trunk.Z,"EdgeAlpha",0.1,"FaceAlpha",0,"EdgeColor","k")
mesh(ax(1),knot.surf.X+knot.centerline.x,knot.surf.Y+knot.centerline.y,knot.surf.Z+knot.centerline.z,"EdgeAlpha",0.1,"FaceAlpha",0,"EdgeColor","k")
plot3(ax(1),knot.centerline.x,knot.centerline.y,knot.centerline.z,"--b")
for i = 1:numel(grps)
plot3(ax(1),rings(i).x,rings(i).y,rings(i).z,"-","Color",[0.8,0,0])
end
copyobj(allchild(ax(1)),ax(2))
copyobj(allchild(ax(1)),ax(3))
view(ax(1),[ 39.0438 17.8241])
view(ax(2),[90,0])
view(ax(3),[0,0])
xlabel(ax,"X-axis")
ylabel(ax,"Y-axis")
zlabel(ax,"Z-axis")