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

Question

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}]

Answer 1

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];

Answer 2

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}]

Answer 3

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]]

Answer 4

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.