11
$\begingroup$

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

filling up a circle

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

filling up a cylinder

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

polyhedron with top cut out

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

random polygon with opening

$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – bRost03 May 28 at 15:30
  • $\begingroup$ @bRost03 v12.1 Windows10. Restarting from a fresh kernel there are no errors on my side, although the Manipulate is a bit slow. $\endgroup$ – flinty May 28 at 15:41
  • $\begingroup$ 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. $\endgroup$ – bRost03 May 28 at 15:47
  • 2
    $\begingroup$ 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. $\endgroup$ – bRost03 May 28 at 16:50
  • $\begingroup$ @bRost03 good point! The second method is perfectly fine. If an algorithm exists for the first case I'd be interested too. $\endgroup$ – flinty May 28 at 17:01
7
$\begingroup$

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

Animation of code working

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];
| improve this answer | |
$\endgroup$
  • $\begingroup$ Really nice. How would you suggest getting the spill point when we don't have 'discarded faces' e.g with an .obj mesh? $\endgroup$ – flinty May 28 at 21:05
  • $\begingroup$ Just find the points on the boundary and take the one with the minimum z value. See my edit $\endgroup$ – bRost03 May 28 at 22:24
  • $\begingroup$ 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. $\endgroup$ – bRost03 May 28 at 22:30
  • $\begingroup$ 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 $\endgroup$ – flinty May 28 at 22:35
  • $\begingroup$ Can you share the specific code you use to generate your mesh object? Can I just take it to be MeshRegion[openShell]? $\endgroup$ – bRost03 May 28 at 23:10
6
$\begingroup$

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]

example of output

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

animated gif of filling

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – flinty May 28 at 21:11
2
$\begingroup$

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

mesh with water filling

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]]
| improve this answer | |
$\endgroup$
-3
$\begingroup$

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

Graphics

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

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

Graphics

Graphics Different opacities.

The magic is this option:

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

Grphics

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

Graphics

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

Graphics

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

Graphics

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Is it possible to use this method to extract the area/volume as the question asks? $\endgroup$ – bRost03 May 28 at 19:19
  • 4
    $\begingroup$ @user2432923: But this doesn't answer the OP's question at all! $\endgroup$ – David G. Stork May 28 at 23:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.