Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Implement pyramid volume and fix wedge volume. #13766

Merged
merged 1 commit into from
May 22, 2022
Merged

Implement pyramid volume and fix wedge volume. #13766

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
254 changes: 137 additions & 117 deletions source/grid/grid_tools_nontemplates.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,126 +130,146 @@ namespace GridTools

return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d));
}
else if (vertex_indices.size() == 5) // pyramid
{
// This remarkably simple formula comes from Equation 4 of
// "Calculation of the volume of a general hexahedron for flow
// predictions", Davies and Salmond, AIAA Journal vol. 23 no. 6.
const auto &x0 = all_vertices[vertex_indices[0]];
const auto &x1 = all_vertices[vertex_indices[1]];
const auto &x2 = all_vertices[vertex_indices[2]];
const auto &x3 = all_vertices[vertex_indices[3]];
const auto &x4 = all_vertices[vertex_indices[4]];

const auto v01 = x1 - x0;
const auto v02 = x2 - x0;
const auto v03 = x3 - x0;
const auto v04 = x4 - x0;
const auto v21 = x2 - x1;

// doing high - low consistently puts us off by -1 from the original
// paper in the first term
return -v04 * cross_product_3d(v21, v03) / 6.0 +
v03 * cross_product_3d(v01, v02) / 12.0;
}
else if (vertex_indices.size() == 6) // wedge
{
/*
The following python/sympy script was used:

#!/usr/bin/env python
# coding: utf-8
import sympy as sp
from sympy.simplify.cse_main import cse
xs = list(sp.symbols(" ".join(["x{}".format(i) for i in range(6)])))
ys = list(sp.symbols(" ".join(["y{}".format(i) for i in range(6)])))
zs = list(sp.symbols(" ".join(["z{}".format(i) for i in range(6)])))
xi, eta, zeta = sp.symbols("xi eta zeta")
tphi = [(1 - xi - eta)*(1 - zeta),
(xi)*(1 - zeta),
(eta)*(1 - zeta),
(1 - xi - zeta)*(zeta),
(xi)*(zeta),
(eta)*(zeta)]
x_real = sum(xs[i]*tphi[i] for i in range(len(xs)))
y_real = sum(ys[i]*tphi[i] for i in range(len(xs)))
z_real = sum(zs[i]*tphi[i] for i in range(len(xs)))
J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]]
for var in [x_real, y_real, z_real]])
detJ = J.det()
detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi)
for x in xs:
detJ2 = detJ2.collect(x)
for y in ys:
detJ2 = detJ2.collect(y)
for z in zs:
detJ2 = detJ2.collect(z)
measure = sp.integrate(sp.integrate(sp.integrate(detJ2, (xi, 0, 1)),
(eta, 0, 1)), (zeta, 0, 1))
measure2 = measure
for vs in [xs, ys, zs]:
for v in vs:
measure2 = measure2.collect(v)
pairs, expression = cse(measure2)
for pair in pairs:
print("const double " + sp.ccode(pair[0]) + " = "
+ sp.ccode(pair[1]) + ";")
print("const double result = " + sp.ccode(expression[0]) + ";")
print("return result;")
/* Script used to generate volume code:

#!/usr/bin/env python
# coding: utf-8
import sympy as sp
from sympy.simplify.cse_main import cse
n_vertices = 6
xs = list(sp.symbols(" ".join(["x{}".format(i)
for i in range(n_vertices)])))
ys = list(sp.symbols(" ".join(["y{}".format(i)
for i in range(n_vertices)])))
zs = list(sp.symbols(" ".join(["z{}".format(i)
for i in range(n_vertices)])))
xi, eta, zeta = sp.symbols("xi eta zeta")
tphi = [(1 - xi - eta)*(1 - zeta),
(xi)*(1 - zeta),
(eta)*(1 - zeta),
(1 - xi - eta)*(zeta),
(xi)*(zeta),
(eta)*(zeta)]
x_real = sum(xs[i]*tphi[i] for i in range(n_vertices))
y_real = sum(ys[i]*tphi[i] for i in range(n_vertices))
z_real = sum(zs[i]*tphi[i] for i in range(n_vertices))
J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]]
for var in [x_real, y_real, z_real]])
detJ = J.det()
detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi)
for x in xs:
detJ2 = detJ2.collect(x)
for y in ys:
detJ2 = detJ2.collect(y)
for z in zs:
detJ2 = detJ2.collect(z)
measure = sp.integrate(sp.integrate(
sp.integrate(detJ2, (eta, 0, 1 - xi)),
(xi, 0, 1)), (zeta, 0, 1))
measure2 = measure
for vs in [xs, ys, zs]:
for v in vs:
measure2 = measure2.collect(v)

pairs, expression = cse(measure2)
for vertex_no in range(n_vertices):
for (coordinate, index) in [('x', 0), ('y', 1), ('z', 2)]:
print(
"const double {}{} = all_vertices[vertex_indices[{}]][{}];"
.format(coordinate, vertex_no, vertex_no, index))

for pair in pairs:
print("const double " + sp.ccode(pair[0]) + " = "
+ sp.ccode(pair[1]) + ";")
print("const double result = " + sp.ccode(expression[0]) + ";")
print("return result;")
*/
const double x0 = all_vertices[vertex_indices[0]](0);
const double y0 = all_vertices[vertex_indices[0]](1);
const double z0 = all_vertices[vertex_indices[0]](2);
const double x1 = all_vertices[vertex_indices[1]](0);
const double y1 = all_vertices[vertex_indices[1]](1);
const double z1 = all_vertices[vertex_indices[1]](2);
const double x2 = all_vertices[vertex_indices[2]](0);
const double y2 = all_vertices[vertex_indices[2]](1);
const double z2 = all_vertices[vertex_indices[2]](2);
const double x3 = all_vertices[vertex_indices[3]](0);
const double y3 = all_vertices[vertex_indices[3]](1);
const double z3 = all_vertices[vertex_indices[3]](2);
const double x4 = all_vertices[vertex_indices[4]](0);
const double y4 = all_vertices[vertex_indices[4]](1);
const double z4 = all_vertices[vertex_indices[4]](2);
const double x5 = all_vertices[vertex_indices[5]](0);
const double y5 = all_vertices[vertex_indices[5]](1);
const double z5 = all_vertices[vertex_indices[5]](2);

const double x6 = (1.0 / 6.0) * z5;
const double x7 = (1.0 / 12.0) * z1;
const double x8 = -x7;
const double x9 = (1.0 / 12.0) * z2;
const double x10 = -x9;
const double x11 = (1.0 / 4.0) * z5;
const double x12 = -x11;
const double x13 = (1.0 / 12.0) * z0;
const double x14 = x12 + x13;
const double x15 = (1.0 / 4.0) * z2;
const double x16 = (1.0 / 6.0) * z4;
const double x17 = (1.0 / 4.0) * z1;
const double x18 = (1.0 / 6.0) * z0;
const double x19 = x17 - x18;
const double x20 = -x6;
const double x21 = (1.0 / 4.0) * z0;
const double x22 = -x21;
const double x23 = -x17;
const double x24 = -x15;
const double x25 = (1.0 / 6.0) * z3;
const double x26 = x24 - x25;
const double x27 = x18 + x23;
const double x28 = (1.0 / 3.0) * z2;
const double x29 = (1.0 / 12.0) * z5;
const double x30 = (1.0 / 12.0) * z3;
const double x31 = -x30;
const double x32 = (1.0 / 4.0) * z4;
const double x33 = x31 + x32;
const double x34 = (1.0 / 3.0) * z1;
const double x35 = (1.0 / 12.0) * z4;
const double x36 = -x16;
const double x37 = x15 + x25;
const double x38 = -x13;
const double x39 = x11 + x38;
const double x40 = -x32;
const double x41 = x30 + x40;
const double x42 = (1.0 / 3.0) * z0;
const double x43 = (1.0 / 4.0) * z3;
const double x44 = x32 - x43;
const double x45 = x40 + x43;
return x0 * (y1 * (-x28 + x29 + x33) + y2 * (x12 + x31 + x34 - x35) +
y3 * (x20 + x7 + x9) + y4 * (x23 + x6 + x9) +
y5 * (x36 + x37 + x8)) +
x1 * (y0 * (x28 - x29 + x41) + y2 * (x11 + x33 - x42) +
y3 * (x39 + x9) + y4 * (x12 + x21 + x24) +
y5 * (x13 + x24 + x44)) +
x2 * (y0 * (x11 + x30 - x34 + x35) + y1 * (x12 + x41 + x42) +
y3 * (x39 + x8) + y4 * (x12 + x17 + x38) +
y5 * (x17 + x22 + x44)) +
x3 * (-x6 * y4 + y0 * (x10 + x6 + x8) + y1 * (x10 + x14) +
y2 * (x14 + x7) + y5 * (x15 + x16 + x19)) +
x4 * (x6 * y3 + y0 * (x10 + x17 + x20) + y1 * (x11 + x15 + x22) +
y2 * (x11 + x13 + x23) + y5 * (x26 + x27)) +
x5 * (y0 * (x16 + x26 + x7) + y1 * (x15 + x38 + x45) +
y2 * (x21 + x23 + x45) + y3 * (x24 + x27 + x36) +
y4 * (x19 + x37));
const double x0 = all_vertices[vertex_indices[0]][0];
const double y0 = all_vertices[vertex_indices[0]][1];
const double z0 = all_vertices[vertex_indices[0]][2];
const double x1 = all_vertices[vertex_indices[1]][0];
const double y1 = all_vertices[vertex_indices[1]][1];
const double z1 = all_vertices[vertex_indices[1]][2];
const double x2 = all_vertices[vertex_indices[2]][0];
const double y2 = all_vertices[vertex_indices[2]][1];
const double z2 = all_vertices[vertex_indices[2]][2];
const double x3 = all_vertices[vertex_indices[3]][0];
const double y3 = all_vertices[vertex_indices[3]][1];
const double z3 = all_vertices[vertex_indices[3]][2];
const double x4 = all_vertices[vertex_indices[4]][0];
const double y4 = all_vertices[vertex_indices[4]][1];
const double z4 = all_vertices[vertex_indices[4]][2];
const double x5 = all_vertices[vertex_indices[5]][0];
const double y5 = all_vertices[vertex_indices[5]][1];
const double z5 = all_vertices[vertex_indices[5]][2];
const double x6 = (1.0 / 12.0) * z1;
const double x7 = -x6;
const double x8 = (1.0 / 12.0) * z3;
const double x9 = x7 + x8;
const double x10 = (1.0 / 12.0) * z2;
const double x11 = -x8;
const double x12 = x10 + x11;
const double x13 = (1.0 / 6.0) * z2;
const double x14 = (1.0 / 12.0) * z4;
const double x15 = (1.0 / 6.0) * z1;
const double x16 = (1.0 / 12.0) * z5;
const double x17 = -x16;
const double x18 = x16 + x7;
const double x19 = -x14;
const double x20 = x10 + x19;
const double x21 = (1.0 / 12.0) * z0;
const double x22 = x19 + x21;
const double x23 = -x10;
const double x24 = x14 + x23;
const double x25 = (1.0 / 6.0) * z0;
const double x26 = x17 + x21;
const double x27 = x23 + x8;
const double x28 = -x21;
const double x29 = x16 + x28;
const double x30 = x17 + x6;
const double x31 = x14 + x28;
const double x32 = x11 + x6;
const double x33 = (1.0 / 6.0) * z5;
const double x34 = (1.0 / 6.0) * z4;
const double x35 = (1.0 / 6.0) * z3;
const double result =
x0 * (x12 * y5 + x9 * y4 + y1 * (-x13 + x14 + x8) +
y2 * (x11 + x15 + x17) + y3 * (x18 + x20)) +
x1 * (x22 * y3 + x24 * y5 + y0 * (x11 + x13 + x19) +
y2 * (x14 + x16 - x25) + y4 * (x26 + x27)) +
x2 * (x29 * y3 + x30 * y4 + y0 * (-x15 + x16 + x8) +
y1 * (x17 + x19 + x25) + y5 * (x31 + x32)) +
x3 * (x26 * y2 + x31 * y1 + y0 * (x24 + x30) + y4 * (x28 + x33 + x7) +
y5 * (x10 + x21 - x34)) +
x4 * (x18 * y2 + x32 * y0 + y1 * (x12 + x29) + y3 * (x21 - x33 + x6) +
y5 * (x23 + x35 + x7)) +
x5 * (x20 * y1 + x27 * y0 + y2 * (x22 + x9) + y3 * (x23 + x28 + x34) +
y4 * (x10 - x35 + x6));
return result;
}

AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell);
Expand Down