Measure of CylinderSurface assumes right-cylinder about meshes.jl HOT 10 CLOSED

I suspect that this existing symbolic function actually works correctly for right cylinders and also cylinders with parallel planes.

For the case of non-parallel planes, numerical integration could be a valid, if expensive, fallback method. There might also be a closed-form general solution to calculate this directly, if someone wants to work on deriving it.

from meshes.jl.

It seems very challenging indeed to get a closed-form solution to this. I think the best plan as you said is to absorb the numerical integration work after we finish the refactoring with support for CRS.

Feel free to close the PR if you feel that a simple solution doesn't exist at the moment. Thanks for the valuable investigation.

from meshes.jl.

I had an idea this morning for how to calculate this in the general (non-intersecting) case. After deriving it out ... it turns out that, so long as you can guarantee that the top and bottom disks don't intersect, the current measure calculations (e.g. $V=h \pi r^2$) actually remain accurate. Essentially, any perturbation in the angle of a plane causes one side of the cylinder to increase by some wedge-shaped volume while the other side of the cylinder decreases by the same. I'd read something like this before, but it was under the assumption of parallel planes both changing by the same angle. Turns out this holds so long as the disks don't intersect, which would then generate an imbalance.

tl:dr; these current measure functions should actually be accurate for any cylinder [surface] so long its top and bottom disks don't intersect.

from meshes.jl.

Got it. Thanks for clarifying it. Closing.

from meshes.jl.

@mikeingold your recent work on Cylinder/CylinderSurface provides a fix for this or we really need numerical integration in the oblique case?

from meshes.jl.

I had a few spare minutes this evening to think about this more. I think I was able to sketch up the math for how to compute this in the general case without integration so long as the planes don’t intersect.

from meshes.jl.

I'm going to have to re-work #806 slightly, but I think have a closed-form solution/algorithm that works in all but the intersecting-planes case (inspired by https://mathworld.wolfram.com/CylindricalSegment.html)

For V = measure(c::Cylinder)

1. If c is a right-cylinder, then use the current solution: $V = h ~ \pi r^2$
2. If c has intersecting planes ... throw a error indicating this case is currently unimplemented.
3. Otherwise
   1. Bisect cylinder c using a plane orthogonal to axis(c) that doesn't intersect c's top or bot planes within the cylinder radius.
   2. Treat the two partitioned sections as new Cylinders whose bot is the bisection plane and top is the outer plane on that side. These two new Cylinders are "cylindrical segments".
   3. Find each cylindrical segment's minimum and maximum side lengths, $h_1$ and $h_2$, using their: radius $r$, top center point $\bar{t}$, normalized central axis $\hat{a}$, and top plane's normal $\hat{n}$. Use that information to calculate each's volume, and return the sum of the two results.

$$[\bar{h_1}, \bar{h_2}] = \bar{t} \pm r ~ \hat{n} \times (\hat{a} \times \hat{n})$$

$$V = \frac{1}{2}\left( h_1 + h_2 \right) \pi r^2 $$

I would probably implement a helper function like _measure_cylindersegment(::Cylinder) to handle the partitioned components, where a cylinder segment is assumed to have an orthogonal bottom plane. The only real edge case to this solution that I'll have to catch is if a segment's top plane is also orthogonal (a right-cylinder), in which case that segments measure reduces to simply measure(c). Overall, this would simplify the third case above to just: partition into cyls = Cylinder[c_1, c_2], and then sum(_measurecylindersegment, cyls).

An analogous process could then be used for A = measure(c::CylinderSurface).

Edit: Update below. There's a significant edge case where this third technique fails.

from meshes.jl.

This is turning out to be a surprisingly-challenging problem. The third method described above (for cylinders with non-parallel planes) should work well in most cases, but it fails if an orthogonal bisection plane can't be determined. In effect: if you imagine a plane orthogonal to axis(c) transiting along the axis from the cylinder's bottom to top, there has to be a stopping point where the intersection of the plane and the cylinder is a circle. For many cylinders this won't be possible, especially with large radii and/or highly-slanted end planes.

I'm out of ideas for covering this non-parallel planes case at the moment. Throwing numerical integration at it, e.g. integral(p -> T(1), c) for this specific conditional case, would definitely work but seems a little heavy-handed. We'd also either need to add MeshIntegrals.jl as a dependency or sort out it's long-term path (transition under JuliaGeometry or directly absorbed into Meshes.jl).

from meshes.jl.

That is really nice. Feel free to submit the PR with the final solution. Looking forward to reviewing it.

from meshes.jl.

I think I could have made the statement above more clear. When I said that the "current" methods actually are accurate, I was referring to the ones already present in the master branch.

The only potential PR-able change we might consider might be to add an if/else condition that errors on hasintersectingplanes(c) since those results won't be accurate.

Other than that, and pending any dissent/feedback, I'd recommend closure of this Issue.

from meshes.jl.