Add a Maxima script deriving quadrature weights on element sides

src/util/fem/q1_3d_boundary.mac

@@ -0,0 +1,68 @@
+/* This script produces formulas for approximating surface integrals
+over faces of 3D Q1 elements in the case using a regular grid in the x
+and y directions. */
+kill(all);
+
+load(vect)$ /* now ~ is the vector cross product */
+
+/* coordinates of the nodes of the reference element */
+xis   : [-1,  1,  1, -1, -1,  1, 1, -1]$
+etas  : [-1, -1,  1,  1, -1, -1, 1,  1]$
+zetas : [-1, -1, -1, -1,  1,  1, 1,  1]$
+
+/* basis functions on the reference element */
+chi[i]:= (1/8) * (1 + xi*xis[i]) * (1 + eta*etas[i]) * (1 + zeta*zetas[i])$
+
+/* The map from the reference to a physical element.
+ *
+ * M(x) = x(xi, eta, zeta).
+ */
+M(v) := sum(v[j] * chi[j], j, 1, 8)$
+
+/* partial derivatives of M
+ *
+ * dM(x, xi) = diff(x(xi, eta, zeta), xi)
+*/
+dM(v, w) := sum(v[j] * diff(chi[j], w), j, 1, 8)$
+
+/* gradient of the map M with respect to xi, eta, zeta, i.e. on the
+ * reference element
+ *
+ * gradM(x) = [diff(x, xi), diff(x, eta), diff(x, zeta)]
+ */
+gradM(v) := [dM(v, xi), dM(v, eta), dM(v, zeta)];
+
+/* parameterizations of faces of the reference hexahedron */
+faces : [
+  [-1, s,  t], [1, s, t],
+  [s, -1,  t], [s, 1, t],
+  [s,  t, -1], [s, t, 1]
+]$
+
+/* equal spacing in the x and y directions */
+equal_spacing : [
+  x[3] = x[2], x[4] = x[1], x[5] = x[1], x[6] = x[2], x[7] = x[2], x[8] = x[1],
+  x[2] = x[1] + DX,
+  y[2] = y[1], y[4] = y[3], y[5] = y[1], y[6] = y[1], y[7] = y[3], y[8] = y[3],
+  y[3] = y[1] + DY
+]$
+depends(z, [s, t, xi, eta, zeta]);
+
+/* D(x, s, 1) = dx/ds on face 1 */
+D(v, w, n) := gradM(v) . diff(faces[n], w);
+
+/* partial derivative of the parameterization R(s, t) on the face n.
+dR(s, n) = dR/ds on face n */
+dR(v, n) := [D(x, v, n), D(y, v, n), diff(z, v)];
+
+/* tangential vectors on the face j */
+T[j] := ratsimp(subst(equal_spacing, [dR(s, j), dR(t, j)]))$
+/* normal vector to the face j */
+N[j] := express(T[j][1] ~ T[j][2])$
+
+/* partial derivative of z with respect to v on face n. For example,
+dz(t, 1) = dz/dt on face 1 */
+dz(v, n) := [diff(z, xi), diff(z, eta), diff(z, zeta)] . diff(faces[n], v);
+
+/* simplified normal vector to the face j */
+NN[j] := subst([diff(z, t) = dz(t, j), diff(z, s) = dz(s, j)], N[j]);