Implement TriaAccessor::measure() for wedges

diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc
index 1bdc089ebe3c99877f19dc48614ca0ace848596d..ec7da05b62eae6cc571852d9cbe79267135853c0 100644 (file)
--- a/source/grid/grid_tools_nontemplates.cc
+++ b/source/grid/grid_tools_nontemplates.cc
@@ -128,6 +128,127 @@ namespace GridTools
 
         return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d));
       }
+    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;")
+         */
+        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));
+      }
 
     AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell);
     // note that this is the