Skip to content

Instantly share code, notes, and snippets.

@jwpeterson
Created December 8, 2013 04:50
Show Gist options
  • Save jwpeterson/7853458 to your computer and use it in GitHub Desktop.
Save jwpeterson/7853458 to your computer and use it in GitHub Desktop.
Compute and verify basis functions and derivatives for 13 node pyramid elements.
# "Serendipity" pyramid element functions, as defined in:
#
# G. Bedrosian, "Shape functions and integration formulas for
# three-dimensional finite element analysis", Int. J. Numerical
# Methods Engineering, vol 35, p. 95-108, 1992.
#
# This pyramid has 13 nodes and a "serendipity" quad on the bottom.
# Bedrosian uses the same size reference element as we want to use,
# but a different node numbering. The purpose of this script is to
# make sure his basis functions check out. Bedrosian also uses
# 1-based numbering so we will copy that here as well.
pyramid13 := proc()
uses LinearAlgebra;
local phi, x, phi_val, i, j, phi_sum, planes, phi_restrict, phi_xi,
phi_eta, phi_zeta, phi_xixi, phi_etaeta, phi_zetazeta,
phi_xieta, phi_xizeta, phi_etazeta, phi_sums, true_phi_sums,
bedrosian_to_libmesh_mapping;
phi := Vector(13);
# Node (1,1,0)
phi[1] := (1/4) * (xi + eta - 1) * ((1+xi)*(1+eta) - zeta + xi*eta*zeta/(1-zeta));
# Node (0,1,0)
phi[2] := (1/2) * (1 + xi - zeta) * (1 - xi - zeta) * (1 + eta - zeta) / (1-zeta);
# Node (-1,1,0)
phi[3] := (1/4) * (-xi + eta - 1) * ((1-xi)*(1+eta) - zeta - xi*eta*zeta/(1-zeta));
# Node (-1,0,0)
phi[4] := (1/2) * (1 + eta - zeta) * (1 - eta - zeta) * (1 - xi - zeta) / (1 - zeta);
# Node (-1,-1,0)
phi[5] := (1/4) * (-xi - eta - 1) * ((1-xi)*(1-eta) - zeta + xi*eta*zeta/(1-zeta));
# Node (0,-1,0)
phi[6] := (1/2) * (1 + xi - zeta) * (1 - xi - zeta) * (1 - eta - zeta) / (1 - zeta);
# Node (1,-1,0)
phi[7] := (1/4) * (xi - eta - 1) * ((1+xi)*(1-eta) - zeta - xi*eta*zeta/(1-zeta));
# Node (1,0,0)
phi[8] := (1/2) * (1 + eta - zeta) * (1 - eta - zeta) * (1 + xi - zeta) / (1-zeta);
# Node (1/2, 1/2, 1/2)
phi[9] := zeta * (1 + xi - zeta) * (1 + eta - zeta) / (1-zeta);
# Node (-1/2, 1/2, 1/2)
phi[10] := zeta * (1 - xi - zeta) * (1 + eta - zeta) / (1-zeta);
# Node (-1/2, 1/2, 1/2)
phi[11] := zeta * (1 - xi - zeta) * (1 - eta - zeta) / (1-zeta);
# Node (1/2, -1/2, 1/2)
phi[12] := zeta * (1 + xi - zeta) * (1 - eta - zeta) / (1-zeta);
# Node (0,0,1)
phi[13] := zeta*(2*zeta - 1);
# The nodal points.
x := Vector(13);
x[1] := Vector(3, [1, 1, 0]);
x[2] := Vector(3, [0, 1, 0]);
x[3] := Vector(3, [-1, 1, 0]);
x[4] := Vector(3, [-1, 0, 0]);
x[5] := Vector(3, [-1, -1, 0]);
x[6] := Vector(3, [0, -1, 0]);
x[7] := Vector(3, [1, -1, 0]);
x[8] := Vector(3, [1, 0, 0]);
x[9] := Vector(3, [1/2, 1/2, 1/2]);
x[10] := Vector(3, [-1/2, 1/2, 1/2]);
x[11] := Vector(3, [-1/2, -1/2, 1/2]);
x[12] := Vector(3, [1/2, -1/2, 1/2]);
x[13] := Vector(3, [0, 0, 1]);
# Check that the basis functions satisfy phi_i(x_j) = delta_ij property
for i from 1 to Dimension(phi) do
# Skip if not implemented
if (phi[i] <> 0) then
printf("Verifying interpolation property of phi[%d]...\n", i);
for j from 1 to Dimension(x) do
# Note that we can't do straight evaluation in zeta since Maple
# reports an error about dividing by zero, so we take limits...
phi_val := limit(subs(xi=x[j][1], eta=x[j][2], phi[i]), zeta=x[j][3]);
# printf("phi[%d] at node %d = %a\n", i, j, phi_val);
# Error checking
if (j = i) then
if (phi_val <> 1) then
printf("Error! phi[%d] has value %a (should be 1) at node %d\n", i, phi_val, j);
# error;
end if;
elif (phi_val <> 0) then
printf("Error! phi[%d] has value %a (should be 0) at node %d\n", i, phi_val, j);
# error;
end if;
end do;
# Print blank line between basis functions.
# printf("\n");
end if;
end do;
printf("\n");
printf("Checking interpolation properties of basis functions...\n");
true_phi_sums := Vector(10, [
1, # 1
xi, # 2
eta, # 3
zeta, # 4
xi^2, # 5
xi*eta, # 6
xi*zeta, # 7
eta^2, # 8
eta*zeta, # 9
zeta^2 # 10
]);
phi_sums := Vector(10);
for i from 1 to Dimension(phi) do
phi_sums[1] := phi_sums[1] + phi[i];
phi_sums[2] := phi_sums[2] + x[i][1]*phi[i];
phi_sums[3] := phi_sums[3] + x[i][2]*phi[i];
phi_sums[4] := phi_sums[4] + x[i][3]*phi[i];
phi_sums[5] := phi_sums[5] + x[i][1]*x[i][1]*phi[i];
phi_sums[6] := phi_sums[6] + x[i][1]*x[i][2]*phi[i];
phi_sums[7] := phi_sums[7] + x[i][1]*x[i][3]*phi[i];
phi_sums[8] := phi_sums[8] + x[i][2]*x[i][2]*phi[i];
phi_sums[9] := phi_sums[9] + x[i][2]*x[i][3]*phi[i];
phi_sums[10] := phi_sums[10] + x[i][3]*x[i][3]*phi[i];
end do;
# Substitute in eps=0
for i from 1 to Dimension(phi_sums) do
phi_sums[i] := simplify(phi_sums[i]);
end do;
# Print interpolation results
for i from 1 to Dimension(phi_sums) do
printf("phi_sums[%d] = %a\n", i, phi_sums[i]);
end do;
# Check interpolation results for errors
for i from 1 to Dimension(phi_sums) do
if (phi_sums[i] <> true_phi_sums[i]) then
printf("Error in phi_sums[%d] = %a\n", i);
error;
end if;
end do;
printf("success!\n");
printf("\n");
# Equations of the faces:
planes := Vector(5);
planes[1] := eta = 1 - zeta; # "back"; nodes 1,3,13; eta = 1 - zeta
planes[2] := xi = -1 + zeta; # "left"; nodes 3,5,13; xi = -1 + zeta
planes[3] := eta = -1 + zeta; # "front"; nodes 5,7,13; eta = -1 + zeta
planes[4] := xi = 1 - zeta; # "right"; nodes 7,1,13; xi = 1 - zeta
planes[5] := zeta = 0; # "bottom"; nodes 1--8; zeta = 0
# We can check the basis functions restricted to the faces are polynomials.
for i from 1 to Dimension(phi) do
printf("Restricting phi[%d]...\n", i);
for j from 1 to Dimension(planes) do
# Do the restriction
phi_restrict := simplify(subs(planes[j], phi[i]));
# printf("phi_restrict=%a\n", phi_restrict);
# Test the restriction for non-trivial denominator
if (type(denom(phi_restrict), integer)=false) then
printf("Error! phi_restrict is not a polynomial!\n");
error;
end if;
end do;
# printf("\n");
end do;
printf("\n");
# Compute first derivatives
phi_xi := Vector(Dimension(phi));
phi_eta := Vector(Dimension(phi));
phi_zeta := Vector(Dimension(phi));
for i from 1 to Dimension(phi) do
# Simplify actually does OK here...
phi_xi[i] := simplify(diff(phi[i], xi));
phi_eta[i] := simplify(diff(phi[i], eta));
# I debated about simplifying phi_zeta, in the end the terms you
# code up are much shorter so I went with it.
phi_zeta[i] := simplify(diff(phi[i], zeta));
end do:
# Compute second derivatives
phi_xixi := Vector(Dimension(phi));
phi_etaeta := Vector(Dimension(phi));
phi_zetazeta := Vector(Dimension(phi));
phi_xieta := Vector(Dimension(phi));
phi_xizeta := Vector(Dimension(phi));
phi_etazeta := Vector(Dimension(phi));
for i from 1 to Dimension(phi) do
phi_xixi[i] := diff(phi_xi[i], xi);
phi_etaeta[i] := diff(phi_eta[i], eta);
phi_zetazeta[i] := simplify(diff(phi_zeta[i], zeta));
phi_xieta[i] := simplify(diff(phi_xi[i], eta));
phi_xizeta[i] := simplify(diff(phi_xi[i], zeta));
phi_etazeta[i] := simplify(diff(phi_eta[i], zeta));
end do:
# Map from Bedrosian's numbering to the libmesh numbering.
# Use this when printing, so that we can get the results in
# an order that will be useful in libmesh.
bedrosian_to_libmesh_mapping := Vector(Dimension(phi),
[5, 7, 1, 3, 13, 6, 8, 2, 4, 11, 12, 9, 10]);
# Print the basis functions
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi[%d] = %a\n", i-1, phi[j]);
end do;
printf("\n");
# Print first derivatives
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_xi[%d] = %a \n", i-1, phi_xi[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_eta[%d] = %a \n", i-1, phi_eta[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_zeta[%d] = %a \n", i-1, phi_zeta[j]);
end do:
# Print second derivatives
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_xixi[%d] = %a \n", i-1, phi_xixi[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_etaeta[%d] = %a \n", i-1, phi_etaeta[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_zetazeta[%d] = %a \n", i-1, phi_zetazeta[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_xieta[%d] = %a \n", i-1, phi_xieta[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_xizeta[%d] = %a \n", i-1, phi_xizeta[j]);
end do:
printf("\n");
for i from 1 to Dimension(phi) do
j := bedrosian_to_libmesh_mapping[i];
printf("phi_etazeta[%d] = %a \n", i-1, phi_etazeta[j]);
end do:
return;
end proc;
# Results - basis functions
# phi[0] = 1/4*(-xi-eta-1)*((1-xi)*(1-eta)-zeta+xi*eta*zeta/(1-zeta))
# phi[1] = 1/4*(-eta+xi-1)*((1+xi)*(1-eta)-zeta-xi*eta*zeta/(1-zeta))
# phi[2] = 1/4*(xi+eta-1)*((1+xi)*(1+eta)-zeta+xi*eta*zeta/(1-zeta))
# phi[3] = 1/4*(eta-xi-1)*((1-xi)*(1+eta)-zeta-xi*eta*zeta/(1-zeta))
# phi[4] = zeta*(2*zeta-1)
# phi[5] = 1/2*(1+xi-zeta)*(1-xi-zeta)*(1-eta-zeta)/(1-zeta)
# phi[6] = 1/2*(1+eta-zeta)*(1-eta-zeta)*(1+xi-zeta)/(1-zeta)
# phi[7] = 1/2*(1+xi-zeta)*(1-xi-zeta)*(1+eta-zeta)/(1-zeta)
# phi[8] = 1/2*(1+eta-zeta)*(1-eta-zeta)*(1-xi-zeta)/(1-zeta)
# phi[9] = zeta*(1-xi-zeta)*(1-eta-zeta)/(1-zeta)
# phi[10] = zeta*(1+xi-zeta)*(1-eta-zeta)/(1-zeta)
# phi[11] = (1+eta-zeta)*(1+xi-zeta)*zeta/(1-zeta)
# phi[12] = zeta*(1-xi-zeta)*(1+eta-zeta)/(1-zeta)
# Results - first derivatives
# phi_xi[0] = 1/4*(-zeta-eta+2*eta*zeta-2*xi+2*zeta*xi+2*eta*xi+zeta^2+eta^2)/(-1+zeta)
# phi_xi[1] = -1/4*(-zeta-eta+2*eta*zeta+2*xi-2*zeta*xi-2*eta*xi+zeta^2+eta^2)/(-1+zeta)
# phi_xi[2] = -1/4*(-zeta+eta-2*eta*zeta+2*xi-2*zeta*xi+2*eta*xi+zeta^2+eta^2)/(-1+zeta)
# phi_xi[3] = 1/4*(-zeta+eta-2*eta*zeta-2*xi+2*zeta*xi-2*eta*xi+zeta^2+eta^2)/(-1+zeta)
# phi_xi[4] = 0
# phi_xi[5] = -(-1+eta+zeta)*xi/(-1+zeta)
# phi_xi[6] = 1/2*(-1+eta+zeta)*(1+eta-zeta)/(-1+zeta)
# phi_xi[7] = (1+eta-zeta)*xi/(-1+zeta)
# phi_xi[8] = -1/2*(-1+eta+zeta)*(1+eta-zeta)/(-1+zeta)
# phi_xi[9] = -(-1+eta+zeta)*zeta/(-1+zeta)
# phi_xi[10] = (-1+eta+zeta)*zeta/(-1+zeta)
# phi_xi[11] = -(1+eta-zeta)*zeta/(-1+zeta)
# phi_xi[12] = (1+eta-zeta)*zeta/(-1+zeta)
# phi_eta[0] = 1/4*(-zeta-2*eta+2*eta*zeta-xi+2*zeta*xi+2*eta*xi+zeta^2+xi^2)/(-1+zeta)
# phi_eta[1] = -1/4*(zeta+2*eta-2*eta*zeta-xi+2*zeta*xi+2*eta*xi-zeta^2-xi^2)/(-1+zeta)
# phi_eta[2] = -1/4*(-zeta+2*eta-2*eta*zeta+xi-2*zeta*xi+2*eta*xi+zeta^2+xi^2)/(-1+zeta)
# phi_eta[3] = 1/4*(zeta-2*eta+2*eta*zeta+xi-2*zeta*xi+2*eta*xi-zeta^2-xi^2)/(-1+zeta)
# phi_eta[4] = 0
# phi_eta[5] = -1/2*(-1+xi+zeta)*(1+xi-zeta)/(-1+zeta)
# phi_eta[6] = (1+xi-zeta)*eta/(-1+zeta)
# phi_eta[7] = 1/2*(-1+xi+zeta)*(1+xi-zeta)/(-1+zeta)
# phi_eta[8] = -(-1+xi+zeta)*eta/(-1+zeta)
# phi_eta[9] = -(-1+xi+zeta)*zeta/(-1+zeta)
# phi_eta[10] = (1+xi-zeta)*zeta/(-1+zeta)
# phi_eta[11] = -(1+xi-zeta)*zeta/(-1+zeta)
# phi_eta[12] = (-1+xi+zeta)*zeta/(-1+zeta)
# phi_zeta[0] = -1/4*(xi+eta+1)*(-1+2*zeta-zeta^2+eta*xi)/(-1+zeta)^2
# phi_zeta[1] = 1/4*(eta-xi+1)*(1-2*zeta+zeta^2+eta*xi)/(-1+zeta)^2
# phi_zeta[2] = 1/4*(xi+eta-1)*(-1+2*zeta-zeta^2+eta*xi)/(-1+zeta)^2
# phi_zeta[3] = -1/4*(eta-xi-1)*(1-2*zeta+zeta^2+eta*xi)/(-1+zeta)^2
# phi_zeta[4] = 4*zeta-1
# phi_zeta[5] = 1/2*(-2+eta+6*zeta+eta*xi^2+eta*zeta^2-6*zeta^2+2*zeta^3-2*eta*zeta)/(-1+zeta)^2
# phi_zeta[6] = -1/2*(2-6*zeta+xi+xi*zeta^2+eta^2*xi+6*zeta^2-2*zeta^3-2*zeta*xi)/(-1+zeta)^2
# phi_zeta[7] = -1/2*(2+eta-6*zeta+eta*xi^2+eta*zeta^2+6*zeta^2-2*zeta^3-2*eta*zeta)/(-1+zeta)^2
# phi_zeta[8] = 1/2*(-2+6*zeta+xi+xi*zeta^2+eta^2*xi-6*zeta^2+2*zeta^3-2*zeta*xi)/(-1+zeta)^2
# phi_zeta[9] = (1-eta-4*zeta-xi-xi*zeta^2-eta*zeta^2+eta*xi+5*zeta^2-2*zeta^3+2*eta*zeta+2*zeta*xi)/(-1+zeta)^2
# phi_zeta[10] = -(-1+eta+4*zeta-xi-xi*zeta^2+eta*zeta^2+eta*xi-5*zeta^2+2*zeta^3-2*eta*zeta+2*zeta*xi)/(-1+zeta)^2
# phi_zeta[11] = (1+eta-4*zeta+xi+xi*zeta^2+eta*zeta^2+eta*xi+5*zeta^2-2*zeta^3-2*eta*zeta-2*zeta*xi)/(-1+zeta)^2
# phi_zeta[12] = -(-1-eta+4*zeta+xi+xi*zeta^2-eta*zeta^2+eta*xi-5*zeta^2+2*zeta^3+2*eta*zeta-2*zeta*xi)/(-1+zeta)^2
# Results - second derivatives
# phi_xixi[0] = 1/4*(-2+2*zeta+2*eta)/(-1+zeta)
# phi_xixi[1] = -1/4*(2-2*eta-2*zeta)/(-1+zeta)
# phi_xixi[2] = -1/4*(2+2*eta-2*zeta)/(-1+zeta)
# phi_xixi[3] = 1/4*(-2+2*zeta-2*eta)/(-1+zeta)
# phi_xixi[4] = 0
# phi_xixi[5] = -(-1+eta+zeta)/(-1+zeta)
# phi_xixi[6] = 0
# phi_xixi[7] = (1+eta-zeta)/(-1+zeta)
# phi_xixi[8] = 0
# phi_xixi[9] = 0
# phi_xixi[10] = 0
# phi_xixi[11] = 0
# phi_xixi[12] = 0
# phi_etaeta[0] = 1/4*(-2+2*zeta+2*xi)/(-1+zeta)
# phi_etaeta[1] = -1/4*(2-2*zeta+2*xi)/(-1+zeta)
# phi_etaeta[2] = -1/4*(2-2*zeta+2*xi)/(-1+zeta)
# phi_etaeta[3] = 1/4*(-2+2*zeta+2*xi)/(-1+zeta)
# phi_etaeta[4] = 0
# phi_etaeta[5] = 0
# phi_etaeta[6] = (1+xi-zeta)/(-1+zeta)
# phi_etaeta[7] = 0
# phi_etaeta[8] = -(-1+xi+zeta)/(-1+zeta)
# phi_etaeta[9] = 0
# phi_etaeta[10] = 0
# phi_etaeta[11] = 0
# phi_etaeta[12] = 0
# phi_zetazeta[0] = 1/2*(xi+eta+1)*eta*xi/(-1+zeta)^3
# phi_zetazeta[1] = -1/2*(eta-xi+1)*eta*xi/(-1+zeta)^3
# phi_zetazeta[2] = -1/2*(xi+eta-1)*eta*xi/(-1+zeta)^3
# phi_zetazeta[3] = 1/2*(eta-xi-1)*eta*xi/(-1+zeta)^3
# phi_zetazeta[4] = 4
# phi_zetazeta[5] = -(1-3*zeta+3*zeta^2-zeta^3+eta*xi^2)/(-1+zeta)^3
# phi_zetazeta[6] = (-1+3*zeta-3*zeta^2+zeta^3+eta^2*xi)/(-1+zeta)^3
# phi_zetazeta[7] = (-1+3*zeta-3*zeta^2+zeta^3+eta*xi^2)/(-1+zeta)^3
# phi_zetazeta[8] = -(1-3*zeta+3*zeta^2-zeta^3+eta^2*xi)/(-1+zeta)^3
# phi_zetazeta[9] = -2*(-1+3*zeta-3*zeta^2+zeta^3+eta*xi)/(-1+zeta)^3
# phi_zetazeta[10] = 2*(1-3*zeta+3*zeta^2-zeta^3+eta*xi)/(-1+zeta)^3
# phi_zetazeta[11] = -2*(-1+3*zeta-3*zeta^2+zeta^3+eta*xi)/(-1+zeta)^3
# phi_zetazeta[12] = 2*(1-3*zeta+3*zeta^2-zeta^3+eta*xi)/(-1+zeta)^3
# phi_xieta[0] = 1/4*(-1+2*zeta+2*xi+2*eta)/(-1+zeta)
# phi_xieta[1] = -1/4*(-1+2*zeta-2*xi+2*eta)/(-1+zeta)
# phi_xieta[2] = -1/4*(1-2*zeta+2*xi+2*eta)/(-1+zeta)
# phi_xieta[3] = 1/4*(1-2*zeta-2*xi+2*eta)/(-1+zeta)
# phi_xieta[4] = 0
# phi_xieta[5] = -xi/(-1+zeta)
# phi_xieta[6] = eta/(-1+zeta)
# phi_xieta[7] = xi/(-1+zeta)
# phi_xieta[8] = -eta/(-1+zeta)
# phi_xieta[9] = -zeta/(-1+zeta)
# phi_xieta[10] = zeta/(-1+zeta)
# phi_xieta[11] = -zeta/(-1+zeta)
# phi_xieta[12] = zeta/(-1+zeta)
# phi_xizeta[0] = -1/4*(-1+2*zeta-zeta^2+eta+2*eta*xi+eta^2)/(-1+zeta)^2
# phi_xizeta[1] = 1/4*(-1+2*zeta-zeta^2+eta-2*eta*xi+eta^2)/(-1+zeta)^2
# phi_xizeta[2] = 1/4*(-1+2*zeta-zeta^2-eta+2*eta*xi+eta^2)/(-1+zeta)^2
# phi_xizeta[3] = -1/4*(-1+2*zeta-zeta^2-eta-2*eta*xi+eta^2)/(-1+zeta)^2
# phi_xizeta[4] = 0
# phi_xizeta[5] = eta*xi/(-1+zeta)^2
# phi_xizeta[6] = -1/2*(1+zeta^2+eta^2-2*zeta)/(-1+zeta)^2
# phi_xizeta[7] = -eta*xi/(-1+zeta)^2
# phi_xizeta[8] = 1/2*(1+zeta^2+eta^2-2*zeta)/(-1+zeta)^2
# phi_xizeta[9] = (-1-zeta^2+eta+2*zeta)/(-1+zeta)^2
# phi_xizeta[10] = -(-1-zeta^2+eta+2*zeta)/(-1+zeta)^2
# phi_xizeta[11] = (1+zeta^2+eta-2*zeta)/(-1+zeta)^2
# phi_xizeta[12] = -(1+zeta^2+eta-2*zeta)/(-1+zeta)^2
# phi_etazeta[0] = -1/4*(-1+2*zeta-zeta^2+xi+2*eta*xi+xi^2)/(-1+zeta)^2
# phi_etazeta[1] = 1/4*(1-2*zeta+zeta^2+xi+2*eta*xi-xi^2)/(-1+zeta)^2
# phi_etazeta[2] = 1/4*(-1+2*zeta-zeta^2-xi+2*eta*xi+xi^2)/(-1+zeta)^2
# phi_etazeta[3] = -1/4*(1-2*zeta+zeta^2-xi+2*eta*xi-xi^2)/(-1+zeta)^2
# phi_etazeta[4] = 0
# phi_etazeta[5] = 1/2*(1+xi^2+zeta^2-2*zeta)/(-1+zeta)^2
# phi_etazeta[6] = -eta*xi/(-1+zeta)^2
# phi_etazeta[7] = -1/2*(1+xi^2+zeta^2-2*zeta)/(-1+zeta)^2
# phi_etazeta[8] = eta*xi/(-1+zeta)^2
# phi_etazeta[9] = (-1-zeta^2+xi+2*zeta)/(-1+zeta)^2
# phi_etazeta[10] = -(1+zeta^2+xi-2*zeta)/(-1+zeta)^2
# phi_etazeta[11] = (1+zeta^2+xi-2*zeta)/(-1+zeta)^2
# phi_etazeta[12] = -(-1-zeta^2+xi+2*zeta)/(-1+zeta)^2
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment