# Area/Volume of a 2D/3D object as it is filled up with water

How could I get the volume of water on a region as it is filled up with water from below? Assume gravity points in some appropriate downward axis like -y (or -z in 3D) so the water fills upward from the base of the object.

In addition, I not only consider closed objects but also open ones, so the water should stop filling when it's about to spill out of the opening.

For some objects this is quite easy as demonstrated below, but I'd like a way to calculate these volumes for a wider class of 2D and even 3D objects:

(* Mathematica messes up the padding on the rhs *)
GraphicsRow[
{Graphics[{
LightBlue,
DiskSegment[{0, 0}, 1, {-\[Pi]/2 - 0.6, -\[Pi]/2 + 0.6}],
Black, AbsoluteThickness[2],
Circle[{0, 0}, 1, {-4 \[Pi]/3, \[Pi]/6}],
AbsoluteThickness[1], Arrow[{{0, 0}, {0, 0.5}}]
}],
Graphics[{
LightBlue, DiskSegment[{0, 0}, 1, {-\[Pi] - \[Pi]/6, \[Pi]/6}],
Black, AbsoluteThickness[2],
Circle[{0, 0}, 1, {-4 \[Pi]/3, \[Pi]/6}],
AbsoluteThickness[1], Arrow[{{0, 0}, {0, 0.5}}]
}]
}
]


I thought of a way to do this for simple closed 3D objects using RegionIntersection with a large cuboid that grows in height. The volume can be calculated using RegionMeasure. In the case of objects with holes / openings however, the interior to intersect is missing and the appropriate point to stop increasing z needs to be figured out.

Manipulate[
Block[{c = Cylinder[{{0, 0, 0}, {1, 3, 5}}, 1/2],
r = Cuboid[{-10, -10, -10}, {10, 10, z}]},
Show[RegionIntersection[DiscretizeGraphics@c, DiscretizeGraphics@r],
Graphics3D[{Opacity[.1],
Cylinder[{{0, 0, 0}, {1, 3, 5}}, 1/2]
}]]], {z, 0, 7}]


Here's a potential bowl-like 3D mesh I'm interested in filling - but note I want a general solution that works on non-convex objects too and arbitrary .obj meshes I can load from disk.

SeedRandom[1234];
(*Generate a random polyhedron with an opening near the top *)
r = RegionUnion[
If[Mean[#[[1]]][[3]] < 0.85, #, Nothing] & /@
MeshPrimitives[RandomPolyhedron[100], 2]];
Graphics3D[{
Red, Arrow[{{0, 0, 0}, {0, 0, 1}}],
Green, r}]


And here's an example random polygon with an opening for the 2D case:

SeedRandom[1234];
(*Generate a random polygon with an opening near the top *)
r = RegionUnion[
If[Max[#[[1]][[All, 2]]] < 0.85, #, Nothing] & /@
MeshPrimitives[RandomPolygon[20], 1]];
Graphics[{Red, Arrow[{{0, 0}, {0, 1}}], Green, r}]


• What version of Mathematica are you using? Your second code block (the Manipulate) generates a bunch of messages but still works. Your final code block does not work at all for me (except the Arrow does display). Commented May 28, 2020 at 15:30
• @bRost03 v12.1 Windows10. Restarting from a fresh kernel there are no errors on my side, although the Manipulate is a bit slow. Commented May 28, 2020 at 15:41
• Interesting, I'm on 12.0 Linux - I do see the documentation says some Geometric functions have been updated from 12.0 to 12.1. It's high time I updated anyway. I'll take another look once 12.1 is up and running. Commented May 28, 2020 at 15:47
• Which filling method are you seeking? Examples for your 2D case: imgur.com/GkanV7o or imgur.com/qYhM4Fj? My current thinking is the first case will be much more difficult than the second. Commented May 28, 2020 at 16:50
• @bRost03 good point! The second method is perfectly fine. If an algorithm exists for the first case I'd be interested too. Commented May 28, 2020 at 17:01

Here is a fairly general solution that should work in any dimensions and will not allow you to overflow.

SeedRandom[1234];
P = RandomPolyhedron[100];
{ps, fs} = {P[[1]], P[[2]]};(*points, faces*)

pos = Flatten@Position[ps, a_List /; a[[3]] > 0.95];(*points with z>0.95*)
fs2 = DeleteCases[fs, a_ /; ContainsAny[a, pos]]; (* remove faces containing points with z>0.95 *)

shell = RegionBoundary[P]; (* get the 2D boundary of the 3D region *)
openShell = Polygon[ps, fs2]; (* make a 2D region from remaining faces *)
missingFs = Complement[shell[[2]], openShell[[2]]]; (* retain discarded faces to identify "spill point" *)

zmin = Min[Map[ps[[#]] &, missingFs, {2}][[All, All, 3]]]; (* "spill point" is smallest z-value in discarded faces *)

cube[z_] = Cuboid[{0, 0, 0}, {1, 1, z}];
volume[z_] = Volume@RegionIntersection[P, cube[z]]; (* volume of intersection between cube and 3D region *)

Manipulate[Graphics3D[{Red, Arrow[{{0, 0, 0}, {0, 0, 1}}], PointSize[Large],
Map[Point[openShell[[1, #]]] &, openShell[[2]], {2}], Green, openShell, Blue,
cube[z],Text[Style[V == volume[z], Black, Bold, 18], {1/2, 1/2, 1.05}]}], {z, 0, zmin}]


If you are just given openShell as a starting point, you may find zmin as follows

Ps=openShell[[1]];
Fs=openShell[[2]];
edges = Sort /@ Partition[Flatten[Subsets[#, {2}] & /@ Fs], 2];
boundaryEdgePos = Position[Tally[edges][[All, -1]], 1];(* boundary edges are part of only 1 polygon *)
boundaryEdges = Extract[Tally[edges][[All, 1]], boundaryEdgePos];
zmin = Min@(Ps[[Union@Flatten@boundaryEdges]])[[All, 3]];
(* get the minimum z position of all the points on the boundary *)


If you are further just given a mesh region, just replace Ps and Fs above with:

openShellMesh = MeshRegion[openShell];
Ps = First /@ MeshCells[openShellMesh, 2];
Fs = MeshCoordinates[openShellMesh];

• Really nice. How would you suggest getting the spill point when we don't have 'discarded faces' e.g with an .obj mesh? Commented May 28, 2020 at 21:05
• Just find the points on the boundary and take the one with the minimum z value. See my edit Commented May 28, 2020 at 22:24
• In fact this is better, because it's possible to have a discarded point not on the boundary - and further possible that this point has a z value below the true spill point. E.g. if the shape were a cube plus a point in the center (so the 'top' face was an upside down pyramid), and the center point were discarded. Commented May 28, 2020 at 22:30
• I'm getting some errors when I replace openShell with an obj mesh: Part specification ...openShell [[2]]is longer than object. If I remove the [[2]] then Flatten throws an error. I need a way to split the mesh into points+faces Commented May 28, 2020 at 22:35
• Can you share the specific code you use to generate your mesh object? Can I just take it to be MeshRegion[openShell]? Commented May 28, 2020 at 23:10

Here is a stab at the 2D problem:

pts = {{0.0828049, 0.790215}, {0.245349, 0.759896},
{0.0862234, 0.377913}, {0.40815, 0.678676},
{0.401549, 0.632741}, {0.543757, 0.479332},
{0.471262, 0.309999}, {0.856038, 0.00781796},
{0.824395, 0.105538}, {0.781802, 0.216368},
{0.583854, 0.263973}, {0.651802, 0.323889},
{0.984993, 0.217045}, {0.91956, 0.423835},
{0.876608, 0.521964}, {0.98729, 0.587943},
{0.696159, 0.751866}};

ClearAll[area]
area[h_?(0 <= # <= 1 &)] :=
{
Show[
{Region[#], Graphics[{Black, Opacity[0.1], Polygon[pts]}]},
PlotRange -> All
],
Area[#]
}&@
RegionIntersection[Polygon[pts], Rectangle[{0, 0}, {1, h}]]

area[0.6]


Animate[area[h], {h, 0, 1}]


• This is a concise and well presented solution to the closed polygon problem. It's similar to my 3D example with the cylinder and cuboid, but do you know how I could fill up a polygon with some part of the top removed? Like if I gave you a BSplineCurve smile shape for instance could you fill it up to the lowest point on the spline that prevents a spill? Commented May 28, 2020 at 21:11

This is adapted from @bRost03 's answer - I take no credit. I've noticed that Mathematica cannot handle the RegionIntersection very well for certain meshes even though I've checked that SolidRegionQ@makesolid[RepairMesh[openShell]] returns True and the capped region has RegionEmbeddingDimension 3.

pts={{{-17.9936,-40.3473,-40.3135},{6.87295,-58.345,-11.2017},{-17.9936,-54.0963,-18.074}},{{100.,85.0578,52.5849},{30.9017,35.0465,133.48},{6.87295,-58.345,-11.2017}},{{-80.9017,54.1491,102.581},{-80.9017,115.966,2.58923},{-17.9936,-40.3473,-40.3135}},{{-80.9017,115.966,2.58923},{30.9017,135.069,-28.3098},{6.87295,-36.0986,-47.1858}},{{30.9017,35.0465,133.48},{-80.9017,54.1491,102.581},{-17.9936,-54.0963,-18.074}},{{30.9017,135.069,-28.3098},{100.,85.0578,52.5849},{22.2413,-47.2218,-29.1938}},{{22.2413,-47.2218,-29.1938},{100.,85.0578,52.5849},{6.87295,-58.345,-11.2017}},{{-17.9936,-54.0963,-18.074},{-80.9017,54.1491,102.581},{-17.9936,-40.3473,-40.3135}},{{-17.9936,-40.3473,-40.3135},{-80.9017,115.966,2.58923},{6.87295,-36.0986,-47.1858}},{{6.87295,-58.345,-11.2017},{30.9017,35.0465,133.48},{-17.9936,-54.0963,-18.074}},{{6.87295,-36.0986,-47.1858},{30.9017,135.069,-28.3098},{22.2413,-47.2218,-29.1938}},{{22.2413,-47.2218,-29.1938},{6.87295,-58.345,-11.2017},{-17.9936,-40.3473,-40.3135}},{{6.87295,-36.0986,-47.1858},{22.2413,-47.2218,-29.1938},{-17.9936,-40.3473,-40.3135}}};
openShell = RegionUnion[Polygon /@ prim];
(* makesolid tries to construct a mesh such that SolidRegionQ returns true *)
makesolid[mesh_] := BoundaryMeshRegion[MeshCoordinates[mesh], MeshCells[mesh, 2]]
(* get the polygons of the object *)
polygons = MeshPrimitives[openShell, 2][[All, 1]];
(* get the edges (all pairs of polygon coordinates each sorted *)
edgesOfPolygons = Flatten[(Sort /@ Subsets[#, {2}]) & /@ polygons, 1];
(* count the edges, select edges that appear once i.e on boundary *)
exposedEdge = Select[Tally[edgesOfPolygons], #[[2]] == 1 &][[All, 1]];
(* the spilling z is the smallest z coordinate in the exposed edges *)
spillz = Min[Flatten[exposedEdge, 1][[All, 3]]];
(* get the bounding box of the object *)
bbox = BoundingRegion[openShell];
minz = bbox[[1, 3]] - 10.;
maxz = bbox[[2, 3]];
(* glue the hole in the mesh shut to create a solid region *)
capped = makesolid[RepairMesh[openShell]];
(* intersect the region with a cuboid from minz-spillz - Mathematica
has problems here if you go all the way to the spill point so
subtract a little from spillz *)
tiny = 10^5 \$MachineEpsilon;
Show[RegionIntersection[capped, ReplacePart[bbox, {2, 3} -> spillz - tiny]],
Graphics3D[{Arrow[{{0, 0, 0}, {0, 0, 200}}], Opacity[.1], openShell}]]


However, even if Mathematica cannot produce a RegionIntersection to take the Volume, it's still possible to use RandomPoint to get a Monte-Carlo estimate of the volume provided capped is SolidRegionQ:

(* Monte-Carlo volume *)
With[{n = 50000, c = ReplacePart[bbox, {2, 3} -> spillz]},
Total[If[RegionMember[capped, #], 1, 0] & /@ RandomPoint[c, n]]/n*
Volume@c]

(* True volume - may not work with some meshes + might need to sub tiny from spillz*)
Volume@RegionIntersection[capped, ReplacePart[bbox, {2, 3} -> spillz]]


There is a built-in option for Plot or Graphics that does the trick very well.

Examples:

Plot3D[Sin[x + y^2], {x, -2, 2}, {y, -2, 2},
RegionFunction -> (#1^2 + #2^2 < 4 &), Filling -> Bottom,
FillingStyle -> Directive[Opacity[0.4], Red]]


The surface is swimming atop the fluid. This nicely opaque.

Plot3D[{Cos[x], Cos[x] + 1}, {x, 0, 2 Pi}, {y, 0, 2 Pi},
Filling -> {1 -> {Bottom, Blue}, 2 -> {Top, Red}}]


Filled boundaries of a wavy channel.

Table[Plot3D[Sin[x + y^2], {x, -2, 2}, {y, -2, 2}, Filling -> Bottom,
FillingStyle -> fs], {fs, {Opacity[0.5], Orange}}]


Different opacities.

The magic is this option:

Plot[{Sin[x], Cos[x]}, {x, 0, 2 Pi}, Filling -> {1 -> {{2}, Yellow}}]


This works with curves, surfaces as boundaries.

This works with Regions too.

For a more realistic appearance of bootles or flask or others:

Plot3D[Sin[x], {x, 0, 2 Pi}, {y, 0, \[Pi]/2},
PlotTheme -> "ThickSurface"]


Hypograph and Epigraph are available:

Plot3D[x^2 + y^2, {x, -3, 3}, {y, -2, 2}, Filling -> {1 -> Top}, PlotRange -> 2]

Plot3D[(x^2/9 + y^2/9) - 9, {x, -3, 3}, {y, -2, 2},
RegionFunction -> (#1^2 + #2^2 < 4 &), Filling -> {1 -> Top},
FillingStyle -> Directive[Opacity[0.85], Red],
PlotRange -> {2, 2, 2}]


Plot3D[(x^2/9 + y^2/9) - 9, {x, -3, 3}, {y, -2, 2},
RegionFunction -> (#1^2 + #2^2 < 4 &), Filling -> {1 -> Bottom},
FillingStyle -> Directive[Opacity[0.85], Red],
PlotRange -> {2, 2, 2}]


Mind the use of the PlotTheme with Thickness requires the separations of the 3D object via Show. The built-ins are really fast and can be used in combination with Manipulate.

The obj format are input as surfaces and there is a difference with them in the procedures to the ones shown here. Again mind not all of the example obj have closed surfaces.

• Is it possible to use this method to extract the area/volume as the question asks? Commented May 28, 2020 at 19:19
• @user2432923: But this doesn't answer the OP's question at all! Commented May 28, 2020 at 23:17