Compute and verify basis functions and derivatives for 13 node pyramid elements.

pyramid13.maple

# "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