New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement pyramid volume and fix wedge volume. #13766
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I assume you have verified the correctness of the answers in the output file?
return -1.0 * v04 * cross_product_3d(v21, v03) / 6.0 + | ||
v03 * cross_product_3d(v01, v02) / 12.0; | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a curious formula whose origin I don't understand. Is it only true for pyramids for which the four points of the base are in a single plane? In that case, the cross-product would have something to do with the area of the base, and the dot product with v04 would be like multiplying the base area by something relating to the height of the pyramid. But I don't understand the origin of the second term. It only contains points of the base. Is it a correct for the nonplanarity of the base? (The second term is zero for a planar base.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is derived in the paper I mention in the comment - I don't have a good intuition for it either. It works for non-planar bases, which is what I picked for the 'general pyramid' case:
dealii/tests/simplex/cell_measure_01.cc
Lines 208 to 213 in 6ec9ecb
std::vector<Point<spacedim>> vertices; | |
vertices.emplace_back(-2, -1, -1.5); | |
vertices.emplace_back(1.2, -0.5, -1); | |
vertices.emplace_back(-0.3, 1.2, -0.5); | |
vertices.emplace_back(1, 1, 0.0); | |
vertices.emplace_back(0.5, 0.5, 1.0); |
the first four points are definitely nonplanar.
DEAL:: | ||
DEAL:basic tetrahedron::dim=3 spacedim=3: | ||
DEAL:basic tetrahedron::measure: 0.166667 | ||
DEAL:basic tetrahedron::measure (via quadrature): 0.166667 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I didn't verify these results by hand - I am assuming that since we get the same values with the measure formula and quadrature that this measure formula is correct.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well, at least, for the arbitrary objects below - for the unit shapes its clear that they give the correct results.
6ec9ecb
to
cee71ac
Compare
You're going to have to fix the indentation. It wants to completely garble your python script in the comment. If I understand you right, you've checked that the volume is correct at least for the case of the reference pyramid. OK to merge if so. |
cee71ac
to
9e316ae
Compare
I was hoping it wouldn't do that. The present version is now in the intersection of valid python and things clang-format likes.
Yup! |
9e316ae
to
deb9174
Compare
Implement pyramid volume and fix wedge volume.
I found a couple of mistakes in the code I wrote to compute wedge volumes - by coincidence it was still correct for the reference wedge.