how to create a volume from the revolution of a variable area trapezoid

[Edit: fix typos in the text. Code and formulas are not changed.]

[Edit: add equation for, and calculation of, hs.]

Thank you for posting an interesting problem.

You say in your original post that you want to "create a volume by revolving a variable area trapezoid around the y axis". However, you say that Figure 1 (figure with an ellipse) in the original post shows the x-y plane. In that case, the +z axis points upward out of that plane of Figure 1, and goes through point S. Is that correct? (It would help if each 2D figure included labels for axes.) You say Fig. 2 in the original post shows the x-z plane, and you say that the axis of rotation passes vertically through point S in Figure 2. Therefore it seems the trapezoid revolves about the z axis. But that contradicts what you said earlier. Please clarify.

You write in the original post

It seems to me, based on your other explanations and comments, that the area given by the equaiton above is the area when

, because it is only when

that

. It seems to me that a more general formula for the area is

Am I understanding the problem correctly?

Your equation for p/R appears on my screen as a bunch of text that looks like it should be a Latex equation. So I will try to make it appear as a formatted equation:

And w, R, and θ are simply constants:

, R=0.0006, w=0.0007908. Then

and the equation above can be written

(Check my math there: did I do that right?) Then we can plot of

versus ϕ, and we can plot

in the x,y plane:

phi=0:pi/20:2*pi; w=7.908e-4; R=6e-4;

% Compute p(phi):

p=((w/2)*cos(phi)+(R^2*(sin(phi).^2+0.5*cos(phi).^2)-(w^2/2)*sin(phi).^2).^0.5)./(sin(phi).^2+0.5*cos(phi).^2);

% Plot results

figure;

subplot(211), plot(phi,p,'-r.')

grid on; xlabel('\phi'); ylabel('p'); title('p(\phi)');

subplot(212), plot(p.*cos(phi),p.*sin(phi),'-r.');

grid on; xlabel('x'); ylabel('y'); title('p(x,y)')

The lower plot above looks consistent with the plot of p(x,y) in Fig.1 of your original post.

Your Figure 3 from the original post is below:

I am not sure why use you use subscripts h_i and h_k, above. I think hi=hk when the angle

is small, so it would make more sense to label both h's in the figure as

. You provide the following equation for the total volume:

We can infer from the equation above that you believe that the volume of the 3-dimensional wedge in the figure above is

I do not think that is correct. IF the thick end of the wedge were perpendicular to the wedge long axis (and it's not perpendicular), then the volume would be

where p is a function of ϕ, and therefore dV is a function of ϕ. That's what I think, but you should check my work.

The thick end of the wedge is not perpendicular to the long axis of the wedge, so the preceding formula needs a correction factor which equals unity when the end is perpendicular, and the correction factor is less than 1, but more than 0, otherwise. The end of the wedge is perpendicular to the long axis only when

and when

.

You say

and

. Then

. It follows from the equation I provided above for

that

. Therefore

.

Robert Demyanovich on 17 Oct 2025

Edited: Robert Demyanovich on 17 Oct 2025

@William Rose: Thank you for taking the time to look at this problem. You are quite right that the vertical axis around which the variable area trapezoid rotates is the "z" axis and not the "y" axis (fixed in the original post). I have also fixed the general equation for the trapezoid area and the equation for

as recommended. Thank you for pointing out these inconsistencies.

Your eqn. 4 where

is set equal to 0.5 is correct. Your second plot is essentially the same as the first plot in my post. So far, so good.

With regards to my 3rd plot,

is the height of the trapezoid at ϕ = 0 degrees. The notation

is meant to show a height for a trapezoid that is a very short distance away. The area of trapezoid k is slightly different. It is less than the area of trapezoid i, which has the largest area of any trapezoid slice. The figure was meant to show that the differential thickness by which the trapezoid area is multiplied is

.

The equation that you show for dV I believe is correct and as you state is inferred from the integral equation in the original post. The thinking is that the volume is simply the integral of the area of the trapezoid (which is a function of ϕ) multiplied by the differential thickness, which I believe is

.

I do not understand what you mean by "IF the thick end of the wedge were perpendicular to the wedge long axis (and it's not perpendicular), then the volume would be". I don't know what this means. My understanding is that the volume is the entire rotation of the variable area trapezoid over 2pi. Then you also wrote "where p is a function of ϕ, and therefore dV is a function of ϕ". But the fact that dV is a function of ϕ is already indicated in your Eq. 6, so I guess I am missing the point.

By thick end of the wedge, do you mean the height? At a specific value of ϕ the height is given by:

so the height of the wedge varies as the trapezoid is rotated around the z axis. I'm not sure why a correction factor would be needed. I'm guessing I'm not understanding what you are saying with regards to the thick end of the wedge. My thinking is that for all trapezoid slices where ϕ is less than 90 degrees or greater than 270 degrees,

. When ϕ is equal to 90 or 270 degrees, then

. For

,

(see the second figure in my original post which shows the trapezoids at ϕ equal to 0 degrees (towards the right of point "S") and ϕ equal to 180 degrees (towards the left of point "S").

William Rose on 17 Oct 2025

Edited: William Rose on 17 Oct 2025

@Robert Demyanovich, thank you for checking my work.

[Edit: run the code so that it generates a figure. Change my conclusion: Robert D is right, no adjustment is needed.]

I now recognize that you are right and I was thinking incorrectly. As you said earlier, there is no need for an adjustment factor for the differential volume, dV.

The figure below shows a view of the object from above, i.e. a slice through the object in the x-y plane, through z=0. Three triangles are highlighted to show that the edge of the triangle that is opposite point S is generally not perpendicular to the long axis of the triangle. But no adjustment to the volume formula is needed, because the width of the triangle, perpendicular to the long axis, is p(t)*dt.

N=40; % number of angular steps in one rotation

w=7.908e-4; R=6e-4; % constants in equation for p(t) [units=length]

B=0.867; % constant in equation for h(t): h(t)=B*p(t) [dimensionless]

dt=2*pi/N; % angle step size [radians]

t=[0:dt:2*pi]'-dt/2; % angle of rotation (vector, length N+1) [radians]

% Compute p(t):

p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);

px=p.*cos(t); py=p.*sin(t);

figure

plot(px,py,'--k'); hold on

plot([0;px(1:2);0],[0;py(1:2);0],'-r')

plot([0;px(6:7);0],[0;py(6:7);0],'-g')

plot([0;px(11:12);0],[0;py(11:12);0],'-b')

xlabel('x'); ylabel('y')

text(0,-4e-5,'S')

grid on; axis equal

Robert Demyanovich on 22 Oct 2025

Edited: Robert Demyanovich on 22 Oct 2025

Well, in the image posted, what is it that is above the red trapezoid line at the top? The image is exactly in the x-z plane, there is no curvature due to the y-plane. So, which heights at other azimuthal angles could be larger than the height at d$\phi$ = 0, since by equation, those height are always less. That radial line is kind of like a backbone and everything drips off the upper red line, for example, in the x-y plane. I guess I just don't see what you are saying regarding the image that I posted (not the Seattle space needle example, which I don't think applies, because in the the image I posted, I can see all of it rather easily. That wouldn't be the case if I was at the foot of the space needle trying to get "perspective" on the entire space needle). It doesn't matter where I am located on the z plane between +7E-04 and -7E-04, I can see all of the image and the greenish blue is pretty much always above the upper red trapezoid line or below the lower red trapezoid line. Of course, if I have my face pressed up on the monitor, then maybe I can't see the entire image. But backing away even a little bit allows me to view the entire image with greenish blue above and below the red lines.

Matt J on 23 Oct 2025

Edited: Matt J on 23 Oct 2025

It looks like William Rose commented on this issue several days ago

You haven't provided a link to his comment, but you might be refering to this comment, describing that the top surface of your volume is not planar. An implication of this is that the central spine of your volume (at

) is in fact not the highest. For every fixed x, the points at the rim of your volume are higher than at the center.

That does indeed explain what you are seeing, but I still wanted to illustrate how, even when the top surface is planar, you can still see the concave envelope, depending on your display settings. In this first figure below, I plot the cutoff cylinder that I showed you how to generate using AxelRot. The top and bottom surfaces are planes perpendicular to the x-z plane, and therefore the central spine of the object really is maximally high in this case for any fixed x.

figure;

plotVolume();

view(-30,20)

Accordingly at view(0,0), the central trapezoid is snug against the blue envelope of the object, as long as the axes' Projection property is orthographic (the default).

figure;

plotVolume();

view(0,0)

However, if we switch to perspective projection, which is how peoples' eyes more naturally see the world, the central trapezoid will never be flush with the upper an lower envelopes of view(0,0), and you get the concave envelope. This is because the edge closer to the observer will always appear higher than the central spine.

figure;

plotVolume();

view(0,0)

set(gca,'Projection','perspective');

function plotVolume()

%Volume parameters (independent)

N=100; %discretize the perimeter into N points

diam=6; %tube diameter

ho=1;

hi=4;

xS=-2; %x-coordinate of S

%Volume parameters (derived)

dh=(hi-ho)/2;

angle=atand(dh/diam);

axMajor=hypot(diam,dh);

hS=ho+2*dh/diam*(xS+diam/2);

%Compute boundary vertices

t=linspace(0,360,N+1)'; t(end)=[];

V=0.5*[axMajor*cosd(t),diam*sind(t),0*t];

XYZ1=AxelRot( (V+[0,0,+ho/2])' ,+angle,[0,-1,0], [-diam/2,0,+ho/2])';

XYZ2=AxelRot( (V+[0,0,-ho/2])' ,-angle,[0,-1,0], [-diam/2,0,-ho/2])';

XYZ=[XYZ1;XYZ2];

%Display

trisurf(convhulln(XYZ), XYZ(:,1), XYZ(:,2), XYZ(:,3), ...

'FaceColor', 'cyan', 'EdgeColor', 'none','FaceAlpha',0.3);

axis equal; xlabel x; ylabel y; zlabel z

camlight

%Movie

i=1;

A=XYZ1(i,:);

B=XYZ2(i,:);

C=[xS,0,-hS/2];

D=[xS,0,+hS/2];

args=num2cell([A;B;C;D;A],1);

[x,y,z]=deal(args{:});

line(x,y,z,Color='r',LineWidth=2);

view(0,0)

camdolly(0,+30,0,'fixtarget','data');

end

Matt J on 24 Oct 2025

Edited: Matt J on 25 Oct 2025

Were you talking to me? I don't think I have a volume computation approach of my own anymore that is distinctive from one of @William Rose's. As far as I am aware there are only two volume computation methods that are finalists, and they are as below, from this comment.

wedge1verts=[hst;hsb;pt([1,2],:);pb([1,2],:)];  % array with 6 vertices of wedge 1
[~,wedge1vol]=convhull(wedge1verts);            % wedge volume from convhull
p1=mean(p(1:2));                                % mean p for wedge 1
dV1=(B*p1^3/3+hs*p1^2/6)*dt;                    % wedge volume from formula
fprintf('convhull volume=%.3e, formula volume=%.3e.\n',wedge1vol,dV1)

The one I would have the most confidence in is the computation of wedge1vol using convhull.

As for his analytical wedge volume formula,

dV1=(B*p1^3/3+hs*p1^2/6)*dt;

I sort of get it, but I have the impression that p1=mean(p(1:2)) is being used as a surrogate for the altitude of the triangular base of the wedge. That might be a legit approximation, but I don't immediately see what guarantees it.

Matt J on 24 Oct 2025

Edited: Matt J on 25 Oct 2025

Based on what I posted for your code, the difference seems very large (7.6E-10 from your code versus 1.29E-09 from the integral equation)

I don't know what you did to arrive at that, but below is my idea of an apples-to-apples comparison. You can see that using my original computation (convexVolume) is about 5% greater than the trueVolume (computed using wedges). That seems consistent with the plot, which visualizes of how much convhulln overestimates the volume of the object. The blue region is the convex hull (the object boundary as interpreted by convhulln). The parts of the red trapezoids that are hidden from view (they sort of go underwater into the blue area) are where the true boundaries of the object go strictly inside the hull.

N=60; % number of angular steps in one rotation

w=7.908e-4; R=6e-4; % constants in equation for p(t) [units=length]

B=0.867; % constant in equation for h(t): h(t)=B*p(t) [dimensionless]

dt=2*pi/N; % angle step size [radians]

t=[0:dt:2*pi]'; % angle of rotation (vector, length N+1) [radians]

% Compute p(t):

p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);

h=B*p; % [length]

hs=B*sqrt(R^2-w^2/2);

% Compute arrays of 3D points

pt=[p.*cos(t),p.*sin(t),h/2]; % x,y,z coords of points along top of shape

pb=[p.*cos(t),p.*sin(t),-h/2]; % x,y,z coords of points along bottom of shape

hst=[0,0,hs/2]; hsb=[0,0,-hs/2]; % x,y,z coords of top, bottom central points

% Plot the trapezoid in 3D at various rotation angles

figure; hold on

trueVolume=0;

for i=1:N

plot3([hst(1),pt(i,1),pb(i,1),hsb(1),hst(1)],...

[hst(2),pt(i,2),pb(i,2),hsb(2),hst(2)],...

[hst(3),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

k=[i,i+1];

if i==N, k=[N,1]; end

wedgeVerts=[hst;hsb;pt(k,:);pb(k,:)]; % array with 6 vertices of wedge 1

[~,dV]=convhull(wedgeVerts); % wedge volume from convhull

trueVolume=trueVolume+dV;

end

plot3(pt(:,1),pt(:,2),pt(:,3),'-g') % curve connecting top points

plot3(pb(:,1),pb(:,2),pb(:,3),'-b') % curve connecting bottom points

xlabel('x'); ylabel('y'); zlabel('z');

axis equal; grid on; view(-45,25)

%Display

V=[pt;pb;hst;hsb];

[K,convexVolume]=convhulln(V);

convexVolume,

convexVolume = 1.3392e-09

trueVolume

trueVolume = 1.2787e-09

trisurf(K, V(:,1), V(:,2), V(:,3), ...

'FaceColor', 'cyan', 'EdgeColor', 'none','FaceAlpha',1);

axis equal; xlabel x; ylabel y; zlabel z

camlight

view(-30,20)

hold off

Paul on 25 Oct 2025

The only contribution of this comment is to demonstrate how the volume integral can be determined and evaluated using the Symbolic Math Toolbox. In principal, one could also convert the integrand to an anonymous function for use with integral, but I didn't try that.

clearvars
syms r z phi z_u real
syms p(phi)

Define the volume integral. The factor of 2 out front follows from the lower bound of the innermost integral being 0 (as I see it the volume is symmetric around z = 0)

V = 2*int(int(r*int(1,z,0,z_u,'Hold',true),r,0,p(phi),'Hold',true),phi,0,2*sym(pi),'Hold',true)

V =

Inner integral

I = int(1,z,0,z_u)

I =

Define the height of the trapezoid, relative to z = 0, as a function of phi and the radial distance

syms h_s h(phi)
z_u = h_s/2 + (h(phi)/2 - h_s/2)/p(phi)*r;

Substitute

I = subs(I)

I =

Middle integral

I = int(r*I,r,0,p(phi))

I =

Define h(phi) and substitute, bring in the factor of 2

syms B

h(phi) = B*p(phi);

I = expand(2*subs(I))

I =

Hence the volume is

V = int(I,phi,0,2*sym(pi))

V =

Expressions for p(phi) and h(phi)

syms R theta w real

eqn = p(phi)/R == (w*sin(theta)^2*cos(phi)/R + sqrt(sin(phi)^2+sin(theta)^2*cos(phi)^2-w^2/R^2*sin(theta)^2*sin(phi)^2))/(sin(phi)^2+sin(theta)^2*cos(phi)^2)

eqn =

eqn = isolate(eqn,p)

eqn =

p(phi) = rhs(eqn);

h(phi) = B*p(phi);

Definition of h_s

h_s = h(sym(pi)/2);

Sub all of the above into the integral

V = subs(V)

V =

Define the numeric parameters of the problem

thetav = sym(pi)./[12;6;4;3;2];
Rv = 0.0006*ones(5,1);
Bv = [0.2601;0.5384;0.8670;1.2589;1.25];
wv = [0.0022692;0.0011838;0.0007908;0.0005308;0];

Sub in parameters

V = subs(V,{theta,R,w,B},{thetav,Rv,wv,Bv});

Evaluate

V = double(vpa(V));
T = table(thetav,Rv,Bv,wv,V,'VariableNames',{'theta','R','B','w','V'})
T = 5×5 table
    theta      R         B           w            V     
    _____    ______    ______    _________    __________

    pi/12    0.0006    0.2601    0.0022692    2.6555e-09
    pi/6     0.0006    0.5384    0.0011838    1.5383e-09
    pi/4     0.0006     0.867    0.0007908    1.2907e-09
    pi/3     0.0006    1.2589    0.0005308    1.2527e-09
    pi/2     0.0006      1.25            0    8.4823e-10

William Rose on 25 Oct 2025

Thank you @Paul.

Paul on 25 Oct 2025

I didn't see any solutions using integral3 (apologies if this has already been done in code attachments, I didn't look at those), which eliminates the need to explicitly define dV.

theta = pi/4;
R = 0.0006;
B = 0.8670;
w = 0.0007908;
p = @(phi) (w*sin(theta).^2.*cos(phi)/R + sqrt(sin(phi).^2+sin(theta).^2.*cos(phi).^2-w^2/R^2.*sin(theta).^2*sin(phi).^2))./(sin(phi).^2+sin(theta).^2.*cos(phi).^2)*R;
h = @(phi) B*p(phi);
h_s = h(pi/2);
z_u = @(phi,r) h_s/2 + (h(phi)/2 - h_s/2)./p(phi).*r;
I = @(phi,r,z) 2*r;
V = integral3(I,0,2*pi,0,p,0,z_u)
V = 1.2907e-09

Sign in to answer this question.

Follow Question

Answer 1

William Rose on 17 Oct 2025

Edited: William Rose on 20 Oct 2025

2 votes

rotateTrapezoid.m

@Robert Demyanovich

{Edit 10/20/2025: Fix error in code below by changing a 1 to a 3. It was

[hst(1),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

and it should be, and now is

[hst(3),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

}

I like the elegant solution of @Matt J.

You indicate in your comments that you are also interested in solving the problem by explicit rotation of the trapezoid.

Here is code that that generates a wire frame view of the shape, by rotating the trapezoid.

N=40; % number of angular steps in one rotation

w=7.908e-4; R=6e-4; % constants in equation for p(t) [units=length]

B=0.867; % constant in equation for h(t): h(t)=B*p(t) [dimensionless]

dt=2*pi/N; % angle step size [radians]

t=[0:dt:2*pi]'; % angle of rotation (vector, length N+1) [radians]

% Compute p(t):

p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);

h=B*p; % [length]

hs=B*sqrt(R^2-w^2/2);

% Compute arrays of 3D points

pt=[p.*cos(t),p.*sin(t),h/2]; % x,y,z coords of points along top of shape

pb=[p.*cos(t),p.*sin(t),-h/2]; % x,y,z coords of points along bottom of shape

hst=[0,0,hs/2]; hsb=[0,0,-hs/2]; % x,y,z coords of top, bottom central points

% Plot the trapezoid in 3D at various rotation angles

figure; hold on

for i=1:N

plot3([hst(1),pt(i,1),pb(i,1),hsb(1),hst(1)],...

[hst(2),pt(i,2),pb(i,2),hsb(2),hst(2)],...

[hst(3),pt(i,3),pb(i,3),hsb(3),hst(3)],'-r.')

end

plot3(pt(:,1),pt(:,2),pt(:,3),'-g') % curve connecting top points

plot3(pb(:,1),pb(:,2),pb(:,3),'-b') % curve connecting bottom points

xlabel('x'); ylabel('y'); zlabel('z');

hold off; axis equal; grid on; view(-45,25)

If you run the code locally, you will be able to rotate the 3D plot by clicking on the 3D rotate tool at top of the plot, then click and drag inside the plot.

Wrap a surface around the points, compute the volume, plot the surface:

pAll=[pt;pb;hst;hsb]; % all the 3D points

[k,v]=convhull(pAll); % wrap the points and compute volume

fprintf('Volume=%.3e\n',v)

Volume=1.323e-09

figure

trisurf(k,pAll(:,1),pAll(:,2),pAll(:,3),'FaceColor','c','EdgeColor','none')

axis equal; camlight

I am attaching a script which has the code above, with more comments and an additional plot.

12 Comments
Show 10 older comments Hide 10 older comments

William Rose on 19 Oct 2025

Edited: William Rose on 19 Oct 2025

[edit: correct typo]

I estimate volume two ways: with the convhull command

[k,v]=convhull(pAll);

and by the formula for dV which I suggested.

% Estimate total volume using formula for dV
dV=(B*p(1:N).^3/3+hs*p(1:N).^2/6)*dt;
Vtot=sum(dV);

Results, when I change N=number of steps in one full rotation, then re-run the script:

fprintf('convhull() volume=%.3e, sum(dV) volume=%.3e. N=%d.\n',v,Vtot,N)

>> rotateTrapezoid

convhull() volume=1.323e-09, sum(dV) volume=1.291e-09. N=40.

>> rotateTrapezoid

convhull() volume=1.345e-09, sum(dV) volume=1.291e-09. N=80.

>> rotateTrapezoid

convhull() volume=1.351e-09, sum(dV) volume=1.291e-09. N=180.

>> rotateTrapezoid

convhull() volume=1.352e-09, sum(dV) volume=1.291e-09. N=360.

>> rotateTrapezoid

convhull() volume=1.352e-09, sum(dV) volume=1.291e-09. N=720.

Discussion

My sum(dV) volume is a lot less than what you report, and pretty close to the convex hull volume.
volume=sum(dV) is constant (1.291e-9) as the number of slices (N) grows and each slice gets smaller.
The volume of the convex hull increases slightly from N=40 to N=180 slices, and after that it is stable to four significant figures.
There is a persistent difference, even at the highest N values, where the convex hull volume is ~5% more than sum(dV) volume. I don't know why.

Matt J on 19 Oct 2025

Edited: Matt J on 20 Oct 2025

Now that we know the solid is not convex and convhull is not applicable, you will have to use a tetrahedralization of the solid to compute the volume. Here is a modification of @William Rose's original loop for doing so,

N=360;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]';         % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
h=B*p;  % [length]
hs=B*sqrt(R^2-w^2/2);
% Compute arrays of 3D points
pt=[p.*cos(t),p.*sin(t),h/2];      % x,y,z coords of points along top of shape
pb=[p.*cos(t),p.*sin(t),-h/2];     % x,y,z coords of points along bottom of shape
hst=[0,0,hs/2]; hsb=[0,0,-hs/2];    % x,y,z coords of top, bottom central points
V=0;
wrappedList=[1:N,1];
for k=1:N
    
    i=wrappedList(k);  
    j=wrappedList(k+1);  
    
    x1=[hst(1),pt(i,1),pb(i,1),hsb(1)];
    y1=[hst(2),pt(i,2),pb(i,2),hsb(2)];
    z1=[hst(3),pt(i,3),pb(i,3),hsb(3)];
    x2=[pt(j,1),pb(j,1)];
    y2=[pt(j,2),pb(j,2)];
    z2=[pt(j,3),pb(j,3)];
    
    TR=delaunayTriangulation([x1,x2]',[y1,y2]',[z1,z2]');
    V=V+tetmeshVolume(TR);
end
volume=V
volume = 1.2903e-09
function V_total = tetmeshVolume(TR)
    % TR can be a delaunayTriangulation or struct with .ConnectivityList and .Points
    T = TR.ConnectivityList;
    X = TR.Points;
    % Extract vertex coordinates for all tetrahedra
    v1 = X(T(:,1), :);
    v2 = X(T(:,2), :);
    v3 = X(T(:,3), :);
    v4 = X(T(:,4), :);
    % Compute per-tetra volumes using vectorized cross and dot
    V = abs(dot(v1 - v4, cross(v2 - v4, v3 - v4, 2), 2)) / 6;
    % Total volume
    V_total = sum(V);
end

William Rose on 20 Oct 2025

@Matt J, @Robert Demyanovich,

Thank you @Matt J for showing how to tetrahedralize a dV, and find its volume as the sum of the volumes of the tetrahedra that make up each dV. The total volume you get by this method is quite a bit less than the volume I got by using a formula for each dV: 1.2144e-09 versus 1.291e-09. So I wonder if I made a mistake when I derived my formula for dV by integrating.

To get more insight, I compute the dVs and total volume three ways:

@Matt J's tetrahedralization way
Using the formula I derived by integrating, with the average value of p for each wedge.
Using convhull() to estimate the volume of each dV

N=360;                   % number of angular steps in one rotation
w=7.908e-4;  R=6e-4;    % constants in equation for p(t) [units=length]
B=0.867;                % constant in equation for h(t): h(t)=B*p(t)  [dimensionless]
dt=2*pi/N;              % angle step size [radians]
t=[0:dt:2*pi]';         % angle of rotation (vector, length N+1) [radians]
% Compute p(t):
p=((w/2)*cos(t)+(R^2*(sin(t).^2+0.5*cos(t).^2)-(w^2/2)*sin(t).^2).^0.5)./(sin(t).^2+0.5*cos(t).^2);
h=B*p;  % [length]
hs=B*sqrt(R^2-w^2/2);
% Compute arrays of 3D points
pt=[p.*cos(t),p.*sin(t),h/2];      % x,y,z coords of points along top of shape
pb=[p.*cos(t),p.*sin(t),-h/2];     % x,y,z coords of points along bottom of shape
hst=[0,0,hs/2]; hsb=[0,0,-hs/2];    % x,y,z coords of top, bottom central points
% 1. Volume by Matt J's tetrahedralization method
dV1=zeros(N,1);
wrappedList=[1:N,1];
for k=1:N
    
    i=wrappedList(k);  
    j=wrappedList(k+1);  
    
    x1=[hst(1),pt(i,1),pb(i,1),hsb(1)];
    y1=[hst(2),pt(i,2),pb(i,2),hsb(2)];
    z1=[hst(1),pt(i,3),pb(i,3),hsb(3)];
    x2=[pt(j,1),pb(j,1)];
    y2=[pt(j,2),pb(j,2)];
    z2=[pt(j,3),pb(j,3)];
    TR=delaunayTriangulation([x1,x2]',[y1,y2]',[z1,z2]');
    dV1(i)=tetmeshVolume(TR);
end
volume1=sum(dV1);
% 2. Volume using a formula for dV
pbar=(p(1:N)+p(2:end))/2;           % mean p for each wedge
dV2=(B*pbar.^3/3+hs*pbar.^2/6)*dt;
volume2=sum(dV2);
% 3. Volume by convhull() of each dV
dV3=zeros(N,1);
for i=1:N
    wedgeVerts=[hst;hsb;pt([i,i+1],:);pb([i,i+1],:)];
    [~,dV3(i)]=convhull(wedgeVerts);
end
volume3=sum(dV3);
fprintf('N=%d: volume 1=%.4e, volume 2=%.4e volume 3=%.4e\n',N,volume1,volume2,volume3)
N=360: volume 1=1.2401e-09, volume 2=1.2905e-09 volume 3=1.2903e-09
function V_total = tetmeshVolume(TR)
% TR can be a delaunayTriangulation or struct with .ConnectivityList and .Points
T = TR.ConnectivityList;
X = TR.Points;
% Extract vertex coordinates for all tetrahedra
v1 = X(T(:,1), :);
v2 = X(T(:,2), :);
v3 = X(T(:,3), :);
v4 = X(T(:,4), :);
% Compute per-tetra volumes using vectorized cross and dot
V = abs(dot(v1 - v4, cross(v2 - v4, v3 - v4, 2), 2)) / 6;
% Total volume
V_total = sum(V);
end

Here are the total volume results, for different values of N (multiply each volume by 1e-9):

N volume 1 volume 2 volume 3

40 1.2144 1.2747 1.2638

80 1.2338 1.2866 1.2839

360 1.2401 1.2905 1.2903

720 1.2403 1.2906 1.2906

Remember that vol.1=Matt J's tetra method; vol.2=formula; vol.3=convhull of each dV.

The tetrahedron method gives a smaller volume by about 4%. The other two methods converge to one another as N gets large. I don;t know why.

Does the difference between the methods vary as a function of rotation angle? To find out, plot dV2/dV1 and dV3/dV1, versus wedge number, when N=360. Therefore wedge number=degrees.

plot(1:N,dV2./dV1,'rx',1:N,dV3./dV1,'b+');

xlabel('Degrees'); ylabel('Ratio'); legend('dV2/dV1','dV3/dV1')

grid on; xlim([0 360]); xticks([0 90 180 270 360])

The plot shows that methods 2 and 3 give a 50% larger value than the tetrahedron method, at 180 degrees. At this angle, the wedges are very small. I'm not sure which method is right, and why the other(s) are wrong.

Sam Chak on 21 Oct 2025

Hi @William Rose

+1 👍 to you and @Matt J. Loved all three solutions as they are practical and well-explained. Initially, I was confused because I thought the OP wanted to estimate the air volume of a section of the circular ventilation duct with a 45° bend, rather than the elliptical duct.

William Rose on 21 Oct 2025

@Sam Chak, nice photo! THank you! Your comme nt means a lot, considering how many nice solutions and explanations you have posted.

Answer 2

Matt J on 17 Oct 2025

Edited: Matt J on 17 Oct 2025

3 votes

Along the lines of what @Mathieu NOE commented, I think it does make more sense to start with an elliptic cylinder and clip off the top and bottom with cutting planes. That will be much easier than approaching it as a solid of revolution:

%Volume parameters

N=1000;

ellipse=scale(nsidedpoly(N),[1.5,1]);

planeSlope=0.1;

planeHeight=0.2;

cylLength=0.5;

%Get the outside vertices

X=ellipse.Vertices(:,1);

Y=ellipse.Vertices(:,2);

Z=linspace(-cylLength,+cylLength,N);

[X,Z]=ndgrid(X,Z);

Y=repmat(Y,1,N);

keep=Z<=planeSlope*X+planeHeight & Z>=-planeSlope*X-planeHeight;

X=X(keep); Y=Y(keep); Z=Z(keep);

%Triangulate and compute volume

[K,volume]=convhulln([X,Y,Z]);

volume,

volume = 1.8849

%Display

trisurf(K, X, Y, Z, ...

'FaceColor', 'cyan', 'EdgeColor', 'none');

axis equal

camlight

view(-15,15)

●

5 Comments
Show 3 older comments Hide 3 older comments

Robert Demyanovich on 17 Oct 2025

Edited: Robert Demyanovich on 17 Oct 2025

@Matt J: There are two goals that I have. The main goal is to verify that my equations will yield the correct pictorial representation of the volume and correctly calculate the actual volume itself. Secondly, I would eventually like to have a "video" that shows a trapezoid slice revolving around the vertical line through "S" so that I can visually demonstrate my thinking of how the volume should be determined.

Visually, your chart appears to be correct. I see that you were able to determine the volume of the solid produced from your code. Can MatLab also measure the volume of the solid produced by rotating the trapezoid around the z axis (i.e. by not requiring the volume equation shown in my original post)? That would be a great help in verifying that the equations are, in fact, calculating the correct volume (or perhaps more accurately, that there is a consistency between the two). I would check this further by varying the angle θ, which then varies the constants that I listed in my original post.

Matt J on 17 Oct 2025

Edited: Matt J on 17 Oct 2025

Can MatLab also measure the volume of the solid produced by rotating the trapezoid around the z axis (i.e. by not requiring the volume equation shown in my original post)?

That is what my code has already done (in numerical approximation), assuming you agree that the solid I have graphed matches what you are looking for.

As for the rest of what you are trying to do, I still think you are barking up the wrong tree by analyzing this as a solid of revolution, including your analytical volume computation. The volume can be decomposed into 3 parts: a central elliptic cylinder section, whose volume is trivially

and two identical upper and lower sections with triangular x-z profiles (see below). The upper and lower sections should be fairly easy triple integrals in Cartesian coordinates. Note from symmetry that you really only need to integrate over one half (illustrated below) of one section

%keep=Z<=planeSlope*X+planeHeight & Z>=planeSlope*min(X)+planeHeight & Y>=0;

Matt J on 17 Oct 2025

There are many videos explaing how a volume can be determined by revolving an area around an axis;

But when the area is changing with angle? That seems non-standard.

Matt J on 17 Oct 2025

Edited: Matt J on 17 Oct 2025

https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection?s_tid=srchtitle

To make the movie, you could use this File Exchange tool,

which computes the intersection of two triangulated surfaces. To triangulate your solid you would give inputs from my code,

TR = triangulation(K,X,Y,Z);

while the second triangulated surface would just be the cut plane you want for a particular frame of your movie. You would represent the plane as a sufficiently big rectangle (two triangles) passing through S.

Answer 3

Matt J on 17 Oct 2025

Edited: Matt J on 19 Oct 2025

https://www.mathworks.com/matlabcentral/fileexchange/30864-3d-rotation-about-shifted-axis

3 votes

Here's a little bit more of a polished version of my original answer, which also generates a movie of the trapezoidal cross-sections. You can adjust the VideoWriter properties to your liking.

This code requires the download of AxelRot,

%Volume parameters (independent)

N=100; %discretize the perimeter into N points

diam=8; %tube diameter

ho=1;

hi=4;

xS=-3; %x-coordinate of S

%Volume parameters (derived)

dh=(hi-ho)/2;

angle=atand(dh/diam);

axMajor=hypot(diam,dh);

hS=ho+2*dh/diam*(xS+diam/2);

%Compute boundary vertices

t=linspace(0,360,N+1)'; t(end)=[];

V=0.5*[axMajor*cosd(t),diam*sind(t),0*t];

XYZ1=AxelRot( (V+[0,0,+ho/2])' ,+angle,[0,-1,0], [-diam/2,0,+ho/2])';

XYZ2=AxelRot( (V+[0,0,-ho/2])' ,-angle,[0,-1,0], [-diam/2,0,-ho/2])';

XYZ=[XYZ1;XYZ2];

%Triangulate and compute volume

[K,volume]=convhulln(XYZ);

volume,

volume = 122.1617

%Display

trisurf(K, XYZ(:,1), XYZ(:,2), XYZ(:,3), ...

'FaceColor', 'cyan', 'EdgeColor', 'none','FaceAlpha',0.5);

axis equal; xlabel x; ylabel y; zlabel z

camlight

view(-30,20)

%Movie

shg

v=VideoWriter('trapezoidMovie.mp4','MPEG-4');

v.Quality=95;

v.FrameRate=30;

v.open;

for i=1:N

A=XYZ1(i,:);

B=XYZ2(i,:);

C=[xS,0,-hS/2];

D=[xS,0,+hS/2];

args=num2cell([A;B;C;D;A],1);

[x,y,z]=deal(args{:});

h=line(x,y,z,Color='r',LineWidth=2);

drawnow;

v.writeVideo(getframe(gcf));

if i<N, delete(h); end

end

v.close

24 Comments
Show 22 older comments Hide 22 older comments

Robert Demyanovich on 19 Oct 2025

Edited: Robert Demyanovich on 19 Oct 2025

@William Rose: Matt J's answer is excellent for creating the video that I was looking for. However, the actual image generated is not quite correct (but it doesn't really matter for the video). The equations that I presented indicate that to the rear of the stagnation point (S), the height does not seem to indicate a linear change (if viewed in the x-z plane). There is a steeper dropoff near the edge of the ellipse at 180 degrees for example. This is not meant to be a criticism of Matt J's work on this problem, as he was working from the original blue figures that I posted (Figs 1 and 2). I really appreciate the help that both of you have given me.

I put pause(5) before the loop in your movie and then maiximized the figure. The resulting video has a size that is almost 3 times larger and appears to be a bit crisper even though the font sizes are smaller.

The only issue I still have is resolving why the equations seem to be calculating a larger volume than Matlab calculates.

William Rose on 19 Oct 2025

Edited: William Rose on 19 Oct 2025

rotateTrapezoid.m

[Edit: add a screenshot which I meant to include earlier.]

I see what you mean about the shape in the x-z plane. When I use the 3D rotate tool to rotate the 3D figure that my code makes, to view it looking down the y axis (az=0, el=0), I get this view:

This view makes clear something I had not realized, but which you obviously did realize: the "top" and "bottom" sides of the volume swept out by the trapezoid are not planes.

Non-planarity of the top and bottom also explains the persistent difference of about 5% between the volume of the convex hull and the volume obtained as the sum of the skinny slices. I should have realized this before. The figure below shows the trapezoids at 0, 45, and -45 degrees of rotation. The view is almost exactly down the y- axis, but not quite, because if it were exactly down the y-axis, the +45 and -45 trapezoids would exactly overlap. Notice how the top and bottom edges of the 0 degree trapezoid below the edges of the +,-45 degree trapezoids, in this view. That means the top surface dips down slightly, i.e. it is concave.

If you play around with the wireframe diagram, rotating with your mouse or touchpad in 3D, you can also appreciate the concave nature of the top and bottom surfaces. Here's a screenshot - but its easier to see the concavity when it's moving around, so try it on your machine, if you're skeptical.

The script I used to draw the overall wire frame and the convex hull and the figure with three trapezoids is attached: rotateTrapezoid.m. It does not save a video file.

Since the shape is concave in places, convhull() overestimates the volume. It also means that the volume obtained by taking planar slices through the top and bottom of the elliptical cylinder is not exactly the volume swept out by the rotating trapezoid.

William Rose on 20 Oct 2025

Edited: William Rose on 20 Oct 2025

[Edit: add .mpeg file made by the script]

This script makes a series of figures similar to the nice ones by @John D'Errico, but the surrounding volume is the convex hull around all the trapezoids, instead of being a cylinder with an elliptical cross section, sliced by two planes. I realize that this convex hull is not perfect, as discussed in an earlier comment. But this convex hull is exactly correct along the non-planar top and bottom edges, which you can see by rotating it to az=0, el=0. Screenshots from the .mpeg file during playback are shown, at the 8th and final frames.

William Rose on 20 Oct 2025

Edited: William Rose on 20 Oct 2025

@Robert Demyanovich, you wrote:

"Are you now of the mind that your original equation for dV is correct? Or are you just using it as a point of reference."

I think my equation for dV is correct,

,

to the extent that each dV has the same length, p, on each long sides. This is a pretty good approximation, when the wedges are skinny. In my code, I use the length of the trailing edge of each wedge to calculate dV, all the way around. This means I overestimate the volume of each dV, for phi=0 to 180, because the leading edge is shorter than the trailing edge. But from phi=180 to 360, I underestimate the volume of each wedge, because the leading edge is longer than the trailing edge. The errors in my dVs from phi=0 to 180 pretty much cancel out the errors from 180 to 360, and that is why vol=sum(dV) does not change (at 4 significant digits), when I vary N from 40 wedges (9 deg each) to 720 (0.5 deg each).

If you need each dV to be more accurate, use p=mean(leading & trailing edge lengths) in the formula for each dV, instead of using p=trailing edge length.

For one wedge, comprised of 6 vertices, I expect the volume from convhull() will be very close to the volume from the formula, if we use the mean value of p. Check this with the following commands, after running rotateTrapezoid, with N=80:

wedge1verts=[hst;hsb;pt([1,2],:);pb([1,2],:)];  % array with 6 vertices of wedge 1
[~,wedge1vol]=convhull(wedge1verts);            % wedge volume from convhull
p1=mean(p(1:2));                                % mean p for wedge 1
dV1=(B*p1^3/3+hs*p1^2/6)*dt;                    % wedge volume from formula
fprintf('convhull volume=%.3e, formula volume=%.3e.\n',wedge1vol,dV1)

convhull volume=1.052e-10, formula volume=1.053e-10.

The two volume estimates are very close, as expected.