Skip to content

Commit

Permalink
Implement pyramid volume and fix wedge volume.
Browse files Browse the repository at this point in the history
  • Loading branch information
@drwells
drwells committed May 21, 2022
1 parent 821ac20 commit cee71ac
Show file tree
Hide file tree
Showing 3 changed files with 287 additions and 153 deletions.
248 changes: 131 additions & 117 deletions source/grid/grid_tools_nontemplates.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,126 +130,140 @@ 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

0 comments on commit cee71ac

Please sign in to comment.