Area/Volume of a 2D/3D object as it is filled up with water
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];
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 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]]