Hello everyone.I'm new to the community, please help me to use matlab to calculate the volume of this 3D Shape

clc
u=linspace(0,1,30);
v=linspace(0,2*pi,30);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+cos(v)).*cos(4.*pi.*u);
y=(1-u).*(3+cos(v)).*sin(4.*pi.*u);
z=3.*u+(1-u).*sin(v);
surf(x,y,z);
xlabel('x')
ylabel('y')
zlabel('z')

[David Goodmanson]

Hi PDA,

Could you comment on whether this is a homework problem or if it has it come up for a different reason?

[Bruno Luong]

Easy, the volume is 0,since you plot a 2D manifold

Answer 1

[Bruno Luong]

u=linspace(0,1,31);
v=linspace(0,2*pi,37);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+cos(v)).*cos(4.*pi.*u);
y=(1-u).*(3+cos(v)).*sin(4.*pi.*u);
z=3.*u+(1-u).*sin(v);

V = [x(:) y(:) z(:)];
[m,n]=size(u);
K = reshape(1:m*n,[m,n]);

V1 = K(1:m-1,1:n-1);
V2 = K(2:m,  1:n-1);
V3 = K(2:m,  2:n);
V4 = K(1:m-1,2:n);
F0 = [V3(:) V2(:) V1(:);
      V1(:) V4(:) V3(:)];

V1 = K(2:m,  1);
V2 = K(1:m-1,1);
V3 = K(ones(m-1,1),1);
Fl = [V3(:) V2(:) V1(:)];

V1 = K(1:m-1,  n);
V2 = K(2:m,n);
V3 = K(ones(m-1,1),n);
Fr = [V3(:) V2(:) V1(:)];

F = [F0; Fl; Fr];
% Volume estimated 
VF = permute(reshape(V(F,:),[size(F) 3]),[3 1 2]);
Vol = 1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1))

Vol = 28.5915

[David Goodmanson]

Hi Bruno,

I believe the volume is 3*pi^2, so the answer above is short by a factor of .9656 which is not too bad. I will post an analytic calc depending on what the OP has to say about the homework issue.

[Bruno Luong]

Hi David I'm interested in how you derive the exact volume.

[David Goodmanson]

HI Bruno, you do not have Mathworks-mediated email and I'm kind of reluctant to post anything yet in case this is a straight homework problem.

[Matt J]

Something involving the divergence theorem perhaps?

[Bruno Luong]

Rather discrete Stoke theorem.

[Sam Chak]

Eureka! +1 to @Bruno Luong. I don't know such a "high school homework" problem to calculate a volume of a rigid spiral conical-like tunnel can be so difficult. I only learned how to calculate the Volume of a Solid of Revolution using the Disk Method in Calculus.

https://en.wikipedia.org/wiki/Solid_of_revolution

Gauss' divergence theorem as suggested by @Matt J can also be considered.

https://www.sangakoo.com/en/unit/volumes-calculation-using-gauss-theorem

To determine such a volume of irregular shape, I'd probably make 3D printing of the object from the STL file, and then use the ancient Archimedes' Water Displacement principle to estimate the volume.

[Bruno Luong]

I would think the discrete form I posted above must have a continuous equivalence.

Consider the 2D parametrization of the surface boundary

xyz : (0,1)x(0,2*pi) -> R^3
(u,v) +-> (x(u,v),y(u,v),z(u,v)) 

We should have this kind of formula

V = 1/3*integral over (0x1) x (0,2*pi) dot(cross(dxyz/du,dxyz/dv), xyz) dudv +
    1/3*integral over the cap of the disc (horn base) of similar cross/for. 

I'm to lazy to check it right now. May be later on this.

[Bruno Luong]

I did also my homework. Here is the continuous surface integral to find the volume. I do numerically rather than symbolic (a little messy since I don't see clearly a simplification), but it seems working just fine:

M=3;
N=4*pi;
P=3;

x=@(u,v) (1-u).*(M+cos(v)).*cos(N*u);
y=@(u,v) (1-u).*(M+cos(v)).*sin(N*u);
z=@(u,v) P.*u+(1-u).*sin(v);

dxdu=@(u,v) (M+cos(v)).*(-cos(N*u)-N*(1-u).*sin(N*u));
dydu=@(u,v) (M+cos(v)).*(-sin(N*u)+N*(1-u).*cos(N*u));
dzdu=@(u,v) P-sin(v);

dxdv=@(u,v) (1-u).*(-sin(v)).*cos(N*u);
dydv=@(u,v) (1-u).*(-sin(v)).*sin(N*u);
dzdv=@(u,v) (1-u).*cos(v);

crossgrad=@(u,v) cross([dxdu(u,v);dydu(u,v);dzdu(u,v)],...
                       [dxdv(u,v);dydv(u,v);dzdv(u,v)],1);
dxyz=@(u,v) [x(u,v);y(u,v);z(u,v)]-[x(0,0);y(0,0);z(0,0)];
dVhelper=@(u,v) dot(crossgrad(u,v),dxyz(u,v),1);
dV=@(u,v) reshape(dVhelper(u(:).',v(:).'),size(u));

Volume = integral2(dV,0,1,0,2*pi)/3

Volume = 29.6088

Answer 2

[Matt J]

The volume can be parametrized by introducing a third variable 0<=w<=1,

x=(1-u).*(3+(w*cos(v))).*cos(4.*pi.*u);

y=(1-u).*(3+(w*cos(v))).*sin(4.*pi.*u);

z=3.*u+(1-u).*(w*sin(v));

Sweeping over w produces concentric horns which partition the volume, as seen by the following,

close all
for w=0.1:0.1:1
plotIt(w);
hold on
end
hold off
view(50,15)

function plotIt(w)

u=linspace(0,1,30);
v=linspace(0,2*pi,30);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+(w*cos(v))).*cos(4.*pi.*u);
y=(1-u).*(3+(w*cos(v))).*sin(4.*pi.*u);
z=3.*u+(1-u).*(w*sin(v));
surf(x,y,z,'FaceAlpha',0.3, 'EdgeColor','none','FaceColor',rand(1,3));
xlabel('x')
ylabel('y')
zlabel('z')

end

To compute the volume, we can therefore integrate the Jacobian determinant of the mapping (x,y,z)-->(u,v,w), which leads to a volume of,

Volume = 0.3333*pi*(sin(2*pi) + 9*pi)

or . This is consistent with findings mentioned by @David Goodmanson, though of course I can't know if he used the same method.

[Bruno Luong]

This (u,v,w) parametrization of the volume is more understable for me. +1

[David Goodmanson]

I did use the same method as Matt, where a big key to appreciating the solution is working out the Jacobian, rather than, say, letting some syms engine do it. Since u goes from 0 to 1, it's a bit cleaner to let u --> 1-u and use the general expressions

x = u(N + rcosv) cos(Ku)    y = u(N + rcosv) sin(Ku)    z = M(1-u) + u rsin(v)    

where in this particular case K = 4pi, N=M=3. The third variable r has been inserted, necessary to define a 3d volume. I'll write out J, same as Matt probably did. It is(x,y,x across, d/du, d/dr, d/dv down)

(N+rcosv) cos(Ku) -u(N+rcosv)Ksin(Ku)  (N+rcosv)sin(Ku) +u(N+rcosv)Kcos(Ku)   -M+rsinv 
ucosv     cos(Ku)                       ucos(v) sin(Ku)                        usinv
-ursin(v) cos(Ku)                     -ursin(v) sin(Ku)                        urcosv

Multiply the first col by sin(Ku), the second col by cos(Ku), subtract them, get two zero elements, etc. The result is

J = K(N+rcosv) u^3 r

The v integral from 0 to 2pi integrates out the cosv term, and the result from the three integrations is

2pi 1/4 1/2 KN = (pi/4)KN.

which is 3pi^2 in this case.

oo From the Jacobian, the result does not depend on M. That's consistent because a translation in the z direction is in the plane of the disc generated by r and v and does not sweep out any area.

oo K is the sweep angle around the z axis and does not have to be an integral number of turns.

oo The expression N+rcosv is the distance of a given point from the z axis. The fact that the rcosv term integrates out means that the centroid of the disk is locally at x = 0, so the mean distance from the z axis is N. Then the v integration gives 2 pi N which is consistent with the theorem of Pappus. It seems that Pappus was a pretty smart guy.

oo How does one know to insert r here (same as w for Matt), instead of some other function? It doesn't matter. You could use r^n or any function of some parameter zeta that runs from a to b, such that h(zeta) is continuous, monotonic, differentiable and is 0 at a and 1 at b. In that case the Jacobian and the integration give

Int{a,b} h dh/dzeta  dzeta  = h^2(zeta) /2  | b,a = 1/2    same as if you use r   

Answer 3

[Matt J]

Just as a hint, it's pretty easy to show that the cross sections of the horn are circle with radius (1-u) and whose center is parametrized in u according to,

xc=3*(1-u)*cos(4.*pi.*u); %"center" of the horn
yc=3*(1-u)*sin(4.*pi.*u);
zc=3*u;

The distance from (xc,yc,yc) to the surface can then be calculated for each (u,v) according to,

x-xc=(1-u)*cos(v)*cos(4.*pi.*u);
y-yc=(1-u)*cos(v)*sin(4.*pi.*u);
z-zc=(1-u)*sin(v);

leading to,

which is a u-dependent constant independent of v. It is also stragithforward to show that the points corresponding to a fixed u but arbitrary v are all perpendicular to . They are therefore planar, equidistant from (xc,yc,yc), and therefore form a circle.

Therefore, the cross-sectional area is likewise a u-dependent quantity,

To get the volume, the task is then to integrate this area function along u with some kind of integration measure m(u).

The question though, is what is the functional form of m(u) that is appropriate here?

[Sam Chak]

+1 to @Matt J. Do you mean the function m(u) that describe the red curve along the surface of the spiral conical-like tunnel?

[Matt J]

@Sam Chak It is related to that, but m(u) is some scalar function of u, not a parametric 3D curve. It's something that measures the incremental thinkness of each of the circular cross-sections of the horn. I haven't worked out it's actual functional form.

[Bruno Luong]

I just wonder the validity of such cross section formulation,

For example if we sweep a constant square along another square, the volume intersecting of two ajadcent edges (interior side) must be removed to get the total volume correctly.

When we can integrate along the line to get the volume?

V = line integral A(l) dl ?

The above example shows that is not always the case if the line has 0 radius of curvature. Is it enough if the curve has second derivative everywhere or continuous to use the line integral to get the correct volume?

[Bruno Luong]

Also take an example of the torus volume, one can consider a circle sweep on a circle. The formula is correct is correct if the cross-section has its center aligned with the sweep line ("center" of the torus volume). But one also consider as the cross section sweeping on any other circle, and in general the length changes. Therefore it is not a question of V = cross-section area x line length in a more complex case of non-disc cross section and curved sweep line.

I would vbe very careful in applying such formula.

[Sam Chak]

@Matt J, Thanks for your explanation.

I'm trying to understand your proposition from the volume calculation of a Torus. Since the cross-sectional area of the horn can be determined, now we just need to trace the arc length from the base of the horn all the way up to the tip.

[Bruno Luong]

2*pi*R is the arclength of the center-line. But why the center line? I could also consider the same torus as sweeping of the disk on any other circle (just shifted the radius).

Also the fact if you consider tiny step, the volume of the infinitesimal thin section somewhat "intersecting" because the sweep line curves. It is not clear to me the volume integral can be derived directly because of that (for the straight line it's OK).

May be there is some theorem somewhere to backup this kind of calculation, but until someone shows it to me, I still have some doubt to apply such separation cross-section x line length formula.

[Matt J]

Since the cross-sectional area of the horn can be determined, now we just need to trace the arc length from the base of the horn all the way up to the tip.

@Sam Chak That works for a torus, but not for a general volume. As Bruno is hinting, the key is that the center-line is perpendicular to the circular cross-section in the case of the torus, but not for the case of the horn.

As a simpler example, consider the 3D convex hull of two circles:

The convex hull C of C1 and C2 is a sheared cylinder whose volume is π. However, if you divide the volume up into circular cuts parallel to C1 and integrate the area, π, of the cut along the center-line of C, you calculate an erroneous volume of .

What you really need to do is project the center-line onto the normals of the circular cross-sections (in this case the cross-section normals are all uniformly the z-axis). When the center-line is not straight, as in the case of the horn, you need to do this projection locally, for each infinitessimally thin cross-section.

[David Goodmanson]

If you describe the volume using a third parameter per Matt's and my previous comments (here is an excerpt from mine with parameter r)

Since u goes from 0 to 1, it's a bit cleaner to let u --> 1-u and use the general expressions

    x = u(N + rcosv) cos(Ku)    y = u(N + rcosv) sin(Ku)    z = M(1-u) + u rsin(v) 

and then use the Jacobian, there is no issue with any disks. But if you want to use the swept-disk model, note that during the rotation involving u, the z axis remains in the plane of the disk. That means that the M term above (which is 3 in the original problem) provides a shearing motion parallel to the disk and does not sweep out any area. So the motion involving M can be dropped, in which case the sweeping motion is always perpendicular to the disk. Consequently the arc length of the path that does not include the z component, works, and the result does not even depend on M. Dropping the shearing motion to begin with is easier than [1] calculating the arc length including z and [2] calculating and multiplying by the cosine of the angle between the path and the normal to the disk.

As far as 'why the center line', I already mentioned the thorem of Pappus but didn't state it explicitly. For an arbitrarily shaped 'disk', the volume equals the cross sectional area times the swept arc length r0*theta, where r0 is the distance from the axis of rotation (z in this case) to the centroid of the disk. For a circle the centroid is of course at the center, but in general it has to be calculated.

[Bruno Luong]

David, do you have a reference/link on the Papus theorem of sweep volume?

[David Goodmanson]

Hi Bruno,

It's in lots of books but you don't have to look any further than 'pappus centroid theorem wiki'

[Bruno Luong]

Thanks. Strangely I have never heard of Papus's theorem. Now I feel less ignorance.

[David Goodmanson]

Hi Bruno, 'ignorance' sounds kind of harsh, although it does seem unusual that you would not have run across the theorem, since it is in so many geometry and calculus books.

[Matt J]

I've never heard of Pappu's theorem either, though now that I have, it's a bit obscure to me why it applies here. Pappu's theorem talks about solids and surfaces of revolution. This volume has no rotational symmetry.

[Bruno Luong]

If you look at Wiki there is a Generalization, by Goodman and Goodman for general curve (and not only for a circle about the revolution axis), but with constant cross section. What we need is an even-more generalization of Pappus-Goddman-Godmann where the cross section could evolves along the curve (main assumptions are still perpendicular, centroid along the curve, non intersecting volume). Not sure if this generalization is done by someone, but intuitively it should stands.

[Matt J]

What we need is an even-more generalization of Pappus-Goddman-Godmann where the cross section could evolves along the curve

Seems pretty easy to get such a result for the case (like this one) where the region can be written in cylindrical coordinates. For the volume, we have

where for each fixed θ and z, the centroid with respect to r is while is the width of the region. In the special case (which holds here), where is independent of z, then, denoting the cross-sectional area as , the integration simplifies to

which is essentially the result that I was looking for in my original answer. So, this does allow for θ -varying cross sections, but admittedly requires that they have certain symmetry in r.

With this result though, it turns out to be very easy to calculate the volume of the horn. The horn has circular cross-sections of radius which are at a distance from the z-axis. We therefore have immediately that and leading to,