Verification of basis functions for the Pyramid14 element.

pyramid14.maple

# Quadratic shape functions, as defined in R. Graglia, "Higher order
# bases on pyramidal elements", IEEE Trans Antennas and Propagation,
# vol 47, no 5, May 1999.
pyramid14 := proc()

uses LinearAlgebra;

local p, phi, x, phi_val, i, j, k, phi_sum, n, xi_vec, true_phi_sums, phi_sums,
      phi_restrict, eqns, phi_xi, phi_eta, phi_zeta, phi_xixi,
      phi_etaeta, phi_zetazeta, phi_xieta, phi_xizeta, phi_etazeta,
      plane_eqns, Jacobian, x_physical;

# Normalized "xi_i" coordinates defined by Graglia.  They are actually
# the planes which define the different faces of the pyramid, hence
# the name "p".
p := Vector(5,
  [ (1/2)*(1 - eta - zeta), # "back",  1
    (1/2)*(1 + xi  - zeta), # "left",  2
    (1/2)*(1 + eta - zeta), # "front", 3
    (1/2)*(1 - xi  - zeta), # "right", 4
    zeta  # "bottom",                  5
  ]);

phi := Vector(14);

# Basis functions for base corner (nodes 0, 1, 2, 3)
j := 4;
for i from 1 to 4 do
  # (i,j) = {(4,1), (1,2), (2,3), (3,4)}
  phi[i] := p[i]*p[j]*((2*p[i]/(1 - p[5] + eps) - 1)*(2*p[j]/(1 - p[5] + eps) - 1) - p[5]/(1 - p[5] + eps));

  # Increment j, but let it wrap around if necessary.
  j := j + 1;

  if (j=5) then
      j := 1;
  end if;
end do;

# Apex basis function (node 4)
phi[4 + 1] := p[5]*(2*p[5] - 1);

# Base edge midpoint basis functions (nodes 5, 6, 7, 8)
# 4, 1, 2
# 1, 2, 3
# 2, 3, 4
# 3, 4, 1
j := 4;
k := 2;
for i from 1 to 4 do
  phi[i + 5] := 4 * p[j] * p[i] * p[k] / (1 - p[5] + eps) * (2*p[i]/(1 - p[5] + eps) - 1);

  # Increment j and k, but let them wrap around if necessary.
  j := j + 1;
  k := k + 1;

  if (j=5) then
      j := 1;
  end if;

  if (k=5) then
      k := 1;
  end if;
end do;

# Lateral edge midpoint basis functions (nodes 9, 10, 11, 12)
# 1, 4
# 2, 1
# 3, 2
# 4, 3
j := 4;
for i from 1 to 4 do
  phi[i + 9] := 4 * p[i] * p[j] * p[5] / (1 - p[5] + eps);

  j := j+1;
  if (j=5) then
      j := 1;
  end if;
end do;

# Base middle node (node 13)
phi[13 + 1] := 16*p[1]*p[2]*p[3]*p[4]/(1 - p[5] + eps)^2;

# Print the basis functions
for i from 1 to Dimension(phi) do
  # Simplification seems to help these...
  phi[i] := simplify(phi[i]);

  # We will print these later.
  # printf("phi[%d] = %a\n", i-1, phi[i]);
end do;

printf("\n");

# The nodal points - we number these so they nest with the PYRAMID13 numbering.
x := Vector(14);

# Base vertices
x[0 + 1]  := Vector(3, [-1, -1,   0]);
x[1 + 1]  := Vector(3, [1,  -1,   0]);
x[2 + 1]  := Vector(3, [1,   1,   0]);
x[3 + 1]  := Vector(3, [-1,  1,   0]);

# apex
x[4 + 1] := Vector(3, [0,    0,    1]);

# mid-edge nodes on base
x[5 + 1]  := Vector(3, [0,  -1,   0]);
x[6 + 1]  := Vector(3, [1,   0,   0]);
x[7 + 1]  := Vector(3, [0,   1,   0]);
x[8 + 1]  := Vector(3, [-1,  0,   0]);

# mid-height vertices
x[9 + 1] := Vector(3, [-1/2, -1/2, 1/2]);
x[10 + 1] := Vector(3, [1/2,  -1/2, 1/2]);
x[11 + 1]  := Vector(3, [1/2,   1/2, 1/2]);
x[12 + 1]  := Vector(3, [-1/2,  1/2, 1/2]);

# Base center node
x[13 + 1] := Vector(3, [0,    0,    0]);


# 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
    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 := subs(eps=0, 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-1, j-1, 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-1, phi_val, j-1);
              # error;
          end if;
      elif (phi_val <> 0) then
         printf("Error! phi[%d] has value %a (should be 0) at node %d\n", i-1, phi_val, j-1);
         # error;
      end if;
    end do;

  # Print blank line between basis functions.
  # printf("\n");
  end if;
end do;


printf("Checking interpolation properties of basis functions...");

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(subs(eps=0, 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");

# We can check whether the basis functions restricted to the faces are
# polynomials.  Note that algsubs, which does a bit more than subs'
# text replacement, is needed since the equations of our planes are
# not in the best form for Maple to use...
printf("Checking that restricting basis functions to faces yields polynomials...");
for i from 1 to Dimension(phi) do
  # printf("Restricting phi[%d]...\n", i-1);
  for j from 1 to Dimension(p) do
    # Do the restriction
    phi_restrict := simplify(subs(eps=0, algsubs(p[j]=0, 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");
       printf("phi_restrict=%a\n", phi_restrict);
       error;
    end if;
  end do;
  # printf("\n");
end do;

printf("success!\n");
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
  phi_xi[i]  := diff(phi[i], xi);
  phi_eta[i]  := diff(phi[i], eta);
  phi_zeta[i]  := diff(phi[i], zeta);
end do:

# Show that the Jacobian of an arbitrary pyramid element with
# _parallelogram base_ is a constant, that is, the element map in this
# case is (at most) affine.
# Note: the Jacobian matrix for the reference element itself should be
# the identity matrix!

# The points of a "physical" element.  We start with the points of the
# reference element, move some of them, and see what happens to the
# Jacobian...
x_physical := x;

## Case 1: Stretch the base to [-2,2]x[-2,2].
## Result: Jacobian=Matrix(3, 3, [[2,0,0],[0,2,0],[0,0,1]])
## det(Jacobian)=4
# x_physical[0 + 1]  := Vector(3, [-2, -2,   0]);
# x_physical[1 + 1]  := Vector(3, [2,  -2,   0]);
# x_physical[2 + 1]  := Vector(3, [2,   2,   0]);
# x_physical[3 + 1]  := Vector(3, [-2,  2,   0]);

## Case 2: Nodes 2 and 3 move to the right by 1 so the base is a parallelogram.
## Jacobian=Matrix(3, 3, [[1,1/2,-1/2],[0,1,0],[0,0,1]])
## det(Jacobian)=1
# x_physical[2 + 1]  := Vector(3, [2,   1,   0]);
# x_physical[3 + 1]  := Vector(3, [0,   1,   0]);

## Case 3: Node 2 _only_ moves, so the base is still planar, but no longer a parallelogram.
## Result: (The Jacobian is no longer constant!)
## Jacobian=Matrix(3, 3, [
## [-1/4*(-5*zeta+5+eta)/(-1+zeta), -1/4*(1+xi-zeta)/(-1+zeta)   , 1/4*(-zeta^2+2*zeta-1+eta*xi)/(-1+zeta)^2],
## [-1/4*(1+eta-zeta)/(-1+zeta),    -1/4*(-5*zeta+xi+5)/(-1+zeta), 1/4*(-zeta^2+2*zeta-1+eta*xi)/(-1+zeta)^2],
## [0,0,1]])
## det(Jacobian)=-1/4*(-6*zeta+xi+eta+6)/(-1+zeta)
# x_physical[2 + 1]  := Vector(3, [2,   2,   0]);

## Case 4: the pyramid apex moves off-center.
## Result: The Jacobian is still constant!
## Jacobian=Matrix(3, 3, [[1,0,1/2],[0,1,1/2],[0,0,1/2]])
## det(Jacobian)=1/2
# x_physical[4 + 1]  := Vector(3, [1/2,   1/2,   1/2]);

## Case 5: Node 2 moves down slightly, so the base is no longer planar
## Result: (The Jacobian is no longer constant!)
## Jacobian=Matrix(3, 3, [
## [1,0,0],
## [0,1,0],
## [1/8*(1+eta-zeta)/(-1+zeta), 1/8*(1+xi-zeta)/(-1+zeta), -1/8*(-9*zeta^2+18*zeta-9+eta*xi)/(-1+zeta)^2]])
## det(Jacobian)=-1/8*(-9*zeta^2+18*zeta-9+eta*xi)/(-1+zeta)^2
## Apex Jacobian=Matrix(3, 3, [[1,0,0],[0,1,0],[-1/8,-1/8,9/8]])
# x_physical[2 + 1]  := Vector(3, [1,   1,   -1/2]);

# Update all the mid-edge nodes so the pyramid still has straight sides.
# The base center node gets moved to the average of the four base vertices.
x_physical[5 + 1] := (1/2)*(x_physical[0 + 1] + x_physical[1 + 1]);
x_physical[6 + 1] := (1/2)*(x_physical[1 + 1] + x_physical[2 + 1]);
x_physical[7 + 1] := (1/2)*(x_physical[2 + 1] + x_physical[3 + 1]);
x_physical[8 + 1] := (1/2)*(x_physical[3 + 1] + x_physical[0 + 1]);

x_physical[9 + 1]  := (1/2)*(x_physical[0 + 1] + x_physical[4 + 1]);
x_physical[10 + 1] := (1/2)*(x_physical[1 + 1] + x_physical[4 + 1]);
x_physical[11 + 1] := (1/2)*(x_physical[2 + 1] + x_physical[4 + 1]);
x_physical[12 + 1] := (1/2)*(x_physical[3 + 1] + x_physical[4 + 1]);

x_physical[13 + 1] := (1/4)*(x_physical[0 + 1] + x_physical[1 + 1] + x_physical[2 + 1] + x_physical[3 + 1]);

# Print locations of midedge/midface nodes
# for i from 6 to 14 do
# printf("x[%d] = %a\n", i-1, x_physical[i]);
# end do;

Jacobian := Matrix(3,3);
for i from 1 to Dimension(phi) do
  # dx/d(xi)
  Jacobian[1,1] := Jacobian[1,1] + x_physical[i][1]*phi_xi[i];

  # dx/d(eta)
  Jacobian[1,2] := Jacobian[1,2] + x_physical[i][1]*phi_eta[i];

  # dx/d(zeta)
  Jacobian[1,3] := Jacobian[1,3] + x_physical[i][1]*phi_zeta[i];

  # dy/d(xi)
  Jacobian[2,1] := Jacobian[2,1] + x_physical[i][2]*phi_xi[i];

  # dy/d(eta)
  Jacobian[2,2] := Jacobian[2,2] + x_physical[i][2]*phi_eta[i];

  # dy/d(zeta)
  Jacobian[2,3] := Jacobian[2,3] + x_physical[i][2]*phi_zeta[i];

  # dz/d(xi)
  Jacobian[3,1] := Jacobian[3,1] + x_physical[i][3]*phi_xi[i];

  # dz/d(eta)
  Jacobian[3,2] := Jacobian[3,2] + x_physical[i][3]*phi_eta[i];

  # dz/d(zeta)
  Jacobian[3,3] := Jacobian[3,3] + x_physical[i][3]*phi_zeta[i];
end do;

# Substitute in eps=0 and simplify the Jacobian
Jacobian := simplify(subs(eps=0, Jacobian));

# Print the result
printf("Jacobian=%a\n", Jacobian);
printf("det(Jacobian)=%a\n", Determinant(Jacobian));

# Compute Jacobian at apex using limits to show there is no singularity
for i from 1 to 3 do:
  for j from 1 to 3 do:
    Jacobian[i,j] := limit(subs(xi=0,eta=0,Jacobian[i,j]), zeta=1);
  end do;
end do;

# Print Jacobian matrix evaluated at pyramid apex
printf("Apex Jacobian=%a\n", Jacobian);
printf("\n");


# 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]  := diff(phi_zeta[i], zeta);

  phi_xieta[i]  := diff(phi_xi[i], eta);
  phi_xizeta[i]  := diff(phi_xi[i], zeta);
  phi_etazeta[i]  := diff(phi_eta[i], zeta);
end do:

# Define a set of substitutions we can subsequently use in all
# the basis functions and derivatives below.
plane_eqns := {-1+eta+zeta=-2*p1, 1+xi-zeta=2*p2, -1+xi+zeta=-2*p4, 1+eta-zeta=2*p3};

# Do some simplifications of the basis functions before printing.
for i from 1 to Dimension(phi) do
  phi[i] := subs(plane_eqns, phi[i]);

  # Subs in eps=0 everywhere.  There was no really convenient way to
  # just substitute eps=0 in the numerator... When implementing, just
  # keep in mind that (-1+zeta+eps) should appear in the denominators
  # rather than (-1+zeta).
  phi[i] := subs(eps=0, phi[i]);

  printf("phi[%d] = %a\n", i-1, phi[i]);
end do;

# Do some simplifications of first derivatives before printing first derivatives.
for i from 1 to Dimension(phi) do
  # Do some specific text substitutions
  phi_xi[i] := subs(plane_eqns, phi_xi[i]);
  phi_eta[i] := subs(plane_eqns, phi_eta[i]);
  phi_zeta[i] := subs(plane_eqns, phi_zeta[i]);

  # Do plain simplification
  phi_xi[i] := simplify(phi_xi[i]);
  phi_eta[i] := simplify(phi_eta[i]);
  # phi_zeta[i] := simplify(phi_zeta[i]); # simplifying phi_zeta makes it a bit worse...

  # One other thing I noticed after plain simplification.
  phi_xi[i] := subs(-p4+p2=xi, phi_xi[i]);
  phi_eta[i] := subs(-p3+p1=-eta, phi_eta[i]);
end do:


# Do some simplifications of second derivatives before printing.
for i from 1 to Dimension(phi) do
  # Do some specific text substitutions
  phi_xixi[i]     := subs(plane_eqns, phi_xixi[i]);
  phi_etaeta[i]   := subs(plane_eqns, phi_etaeta[i]);
  phi_zetazeta[i] := subs(plane_eqns, phi_zetazeta[i]);
  phi_xieta[i]    := subs(plane_eqns, phi_xieta[i]);
  phi_xizeta[i]   := subs(plane_eqns, phi_xizeta[i]);
  phi_etazeta[i]  := subs(plane_eqns, phi_etazeta[i]);
end do:

# Subs in eps=0 everywhere.  There was no really convenient way to
# just substitute eps=0 in the numerator... When implementing, just
# keep in mind that (-1+zeta+eps) should appear in the denominators
# rather than (-1+zeta).
for i from 1 to Dimension(phi) do
  phi_xi[i]  := subs(eps=0, phi_xi[i]);
  phi_eta[i]  := subs(eps=0, phi_eta[i]);
  phi_zeta[i]  := subs(eps=0, phi_zeta[i]);
end do:

# Print first derivatives
for i from 1 to Dimension(phi) do
  printf("phi_xi[%d] = %a \n", i-1, phi_xi[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_eta[%d] = %a \n", i-1, phi_eta[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_zeta[%d] = %a \n", i-1, phi_zeta[i]);
end do:


# Subs in eps=0 everywhere.  There was no really convenient way to
# just substitute eps=0 in the numerator... When implementing, just
# keep in mind that (-1+zeta+eps) should appear in the denominators
# rather than (-1+zeta).
for i from 1 to Dimension(phi) do
  phi_xixi[i]  := subs(eps=0, phi_xixi[i]);
  phi_etaeta[i]  := subs(eps=0, phi_etaeta[i]);
  phi_zetazeta[i]  := subs(eps=0, phi_zetazeta[i]);

  phi_xieta[i]  := subs(eps=0, phi_xieta[i]);
  phi_xizeta[i]  := subs(eps=0, phi_xizeta[i]);
  phi_etazeta[i]  := subs(eps=0, phi_etazeta[i]);
end do:


printf("\n");

# Print second derivatives
for i from 1 to Dimension(phi) do
  printf("phi_xixi[%d] = %a \n", i-1, phi_xixi[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_etaeta[%d] = %a \n", i-1, phi_etaeta[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_zetazeta[%d] = %a \n", i-1, phi_zetazeta[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_xieta[%d] = %a \n", i-1, phi_xieta[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_xizeta[%d] = %a \n", i-1, phi_xizeta[i]);
end do:

printf("\n");

for i from 1 to Dimension(phi) do
  printf("phi_etazeta[%d] = %a \n", i-1, phi_etazeta[i]);
end do:

printf("\n");


return;

end proc;

# Results - basis functions

# phi[0] = p4*p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^2
# phi[1] = -p1*p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^2
# phi[2] = p2*p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^2
# phi[3] = -p3*p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2
# phi[4] = zeta*(2*zeta-1)
# phi[5] = -4*p2*p1*p4*eta/(-1+zeta)^2
# phi[6] = 4*p1*p2*p3*xi/(-1+zeta)^2
# phi[7] = 4*p2*p3*p4*eta/(-1+zeta)^2
# phi[8] = -4*p3*p4*p1*xi/(-1+zeta)^2
# phi[9] = -4*p1*p4*zeta/(-1+zeta)
# phi[10] = -4*p2*p1*zeta/(-1+zeta)
# phi[11] = -4*p3*p2*zeta/(-1+zeta)
# phi[12] = -4*p4*p3*zeta/(-1+zeta)
# phi[13] = 16*p1*p2*p3*p4/(-1+zeta)^2




# Results - first derivatives

# phi_xi[0] = 1/2*p1*(-xi*eta+zeta-zeta^2+2*p4*eta)/(-1+zeta)^2
# phi_xi[1] = -1/2*p1*(xi*eta+zeta-zeta^2+2*p2*eta)/(-1+zeta)^2
# phi_xi[2] = 1/2*p3*(xi*eta-zeta+zeta^2+2*p2*eta)/(-1+zeta)^2
# phi_xi[3] = -1/2*p3*(-xi*eta-zeta+zeta^2+2*p4*eta)/(-1+zeta)^2
# phi_xi[4] = 0
# phi_xi[5] = 2*p1*eta*xi/(-1+zeta)^2
# phi_xi[6] = 2*p1*p3*(xi+2*p2)/(-1+zeta)^2
# phi_xi[7] = -2*p3*eta*xi/(-1+zeta)^2
# phi_xi[8] = -2*p1*p3*(-xi+2*p4)/(-1+zeta)^2
# phi_xi[9] = 2*p1*zeta/(-1+zeta)
# phi_xi[10] = -2*p1*zeta/(-1+zeta)
# phi_xi[11] = -2*p3*zeta/(-1+zeta)
# phi_xi[12] = 2*p3*zeta/(-1+zeta)
# phi_xi[13] = -8*p1*p3*xi/(-1+zeta)^2
#
# phi_eta[0] = -1/2*p4*(xi*eta-zeta+zeta^2-2*p1*xi)/(-1+zeta)^2
# phi_eta[1] = 1/2*p2*(xi*eta+zeta-zeta^2-2*p1*xi)/(-1+zeta)^2
# phi_eta[2] = 1/2*p2*(xi*eta-zeta+zeta^2+2*p3*xi)/(-1+zeta)^2
# phi_eta[3] = -1/2*p4*(xi*eta+zeta-zeta^2+2*p3*xi)/(-1+zeta)^2
# phi_eta[4] = 0
# phi_eta[5] = 2*p2*p4*(eta-2*p1)/(-1+zeta)^2
# phi_eta[6] = -2*p2*xi*eta/(-1+zeta)^2
# phi_eta[7] = 2*p2*p4*(eta+2*p3)/(-1+zeta)^2
# phi_eta[8] = 2*p4*xi*eta/(-1+zeta)^2
# phi_eta[9] = 2*p4*zeta/(-1+zeta)
# phi_eta[10] = 2*p2*zeta/(-1+zeta)
# phi_eta[11] = -2*p2*zeta/(-1+zeta)
# phi_eta[12] = -2*p4*zeta/(-1+zeta)
# phi_eta[13] = -8*p2*p4*eta/(-1+zeta)^2
#
# phi_zeta[0] = -1/2*p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-1/2*p4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2+p4*p1*(2*zeta-1)/(-1+zeta)^2-2*p4*p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^3
# phi_zeta[1] = 1/2*p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+1/2*p1*(xi*eta+zeta-zeta^2)/(-1+zeta)^2-p1*p2*(1-2*zeta)/(-1+zeta)^2+2*p1*p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^3
# phi_zeta[2] = -1/2*p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-1/2*p2*(xi*eta-zeta+zeta^2)/(-1+zeta)^2+p2*p3*(2*zeta-1)/(-1+zeta)^2-2*p2*p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^3
# phi_zeta[3] = 1/2*p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+1/2*p3*(xi*eta+zeta-zeta^2)/(-1+zeta)^2-p3*p4*(1-2*zeta)/(-1+zeta)^2+2*p3*p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^3
# phi_zeta[4] = 4*zeta-1
# phi_zeta[5] = 2*p4*p1*eta/(-1+zeta)^2+2*p2*p4*eta/(-1+zeta)^2+2*p1*p2*eta/(-1+zeta)^2+8*p2*p1*p4*eta/(-1+zeta)^3
# phi_zeta[6] = -2*p2*p3*xi/(-1+zeta)^2-2*p1*p3*xi/(-1+zeta)^2-2*p1*p2*xi/(-1+zeta)^2-8*p1*p2*p3*xi/(-1+zeta)^3
# phi_zeta[7] = -2*p3*p4*eta/(-1+zeta)^2-2*p2*p4*eta/(-1+zeta)^2-2*p2*p3*eta/(-1+zeta)^2-8*p2*p3*p4*eta/(-1+zeta)^3
# phi_zeta[8] = 2*p4*p1*xi/(-1+zeta)^2+2*p1*p3*xi/(-1+zeta)^2+2*p3*p4*xi/(-1+zeta)^2+8*p3*p4*p1*xi/(-1+zeta)^3
# phi_zeta[9] = 2*p4*zeta/(-1+zeta)+2*p1*zeta/(-1+zeta)-4*p1*p4/(-1+zeta)+4*p1*p4*zeta/(-1+zeta)^2
# phi_zeta[10] = 2*p1*zeta/(-1+zeta)+2*p2*zeta/(-1+zeta)-4*p2*p1/(-1+zeta)+4*p2*p1*zeta/(-1+zeta)^2
# phi_zeta[11] = 2*p2*zeta/(-1+zeta)+2*p3*zeta/(-1+zeta)-4*p3*p2/(-1+zeta)+4*p3*p2*zeta/(-1+zeta)^2
# phi_zeta[12] = 2*p3*zeta/(-1+zeta)+2*p4*zeta/(-1+zeta)-4*p4*p3/(-1+zeta)+4*p4*p3*zeta/(-1+zeta)^2
# phi_zeta[13] = -8*p2*p3*p4/(-1+zeta)^2-8*p3*p4*p1/(-1+zeta)^2-8*p2*p1*p4/(-1+zeta)^2-8*p1*p2*p3/(-1+zeta)^2-32*p1*p2*p3*p4/(-1+zeta)^3




# Results - second derivatives

# phi_xixi[0] = -p1*eta/(-1+zeta)^2
# phi_xixi[1] = -p1*eta/(-1+zeta)^2
# phi_xixi[2] = p3*eta/(-1+zeta)^2
# phi_xixi[3] = p3*eta/(-1+zeta)^2
# phi_xixi[4] = 0
# phi_xixi[5] = 2*p1*eta/(-1+zeta)^2
# phi_xixi[6] = 4*p1*p3/(-1+zeta)^2
# phi_xixi[7] = -2*p3*eta/(-1+zeta)^2
# phi_xixi[8] = 4*p1*p3/(-1+zeta)^2
# phi_xixi[9] = 0
# phi_xixi[10] = 0
# phi_xixi[11] = 0
# phi_xixi[12] = 0
# phi_xixi[13] = -8*p1*p3/(-1+zeta)^2

# phi_etaeta[0] = -p4*xi/(-1+zeta)^2
# phi_etaeta[1] = p2*xi/(-1+zeta)^2
# phi_etaeta[2] = p2*xi/(-1+zeta)^2
# phi_etaeta[3] = -p4*xi/(-1+zeta)^2
# phi_etaeta[4] = 0
# phi_etaeta[5] = 4*p2*p4/(-1+zeta)^2
# phi_etaeta[6] = -2*p2*xi/(-1+zeta)^2
# phi_etaeta[7] = 4*p2*p4/(-1+zeta)^2
# phi_etaeta[8] = 2*p4*xi/(-1+zeta)^2
# phi_etaeta[9] = 0
# phi_etaeta[10] = 0
# phi_etaeta[11] = 0
# phi_etaeta[12] = 0
# phi_etaeta[13] = -8*p2*p4/(-1+zeta)^2

# phi_zetazeta[0] = 1/2*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-p1*(2*zeta-1)/(-1+zeta)^2+2*p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-p4*(2*zeta-1)/(-1+zeta)^2+2*p4*(xi*eta-zeta+zeta^2)/(-1+zeta)^3+2*p4*p1/(-1+zeta)^2-4*p4*p1*(2*zeta-1)/(-1+zeta)^3+6*p4*p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^4
# phi_zetazeta[1] = -1/2*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+p2*(1-2*zeta)/(-1+zeta)^2-2*p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+p1*(1-2*zeta)/(-1+zeta)^2-2*p1*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+2*p1*p2/(-1+zeta)^2+4*p1*p2*(1-2*zeta)/(-1+zeta)^3-6*p1*p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^4
# phi_zetazeta[2] = 1/2*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-p3*(2*zeta-1)/(-1+zeta)^2+2*p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-p2*(2*zeta-1)/(-1+zeta)^2+2*p2*(xi*eta-zeta+zeta^2)/(-1+zeta)^3+2*p2*p3/(-1+zeta)^2-4*p2*p3*(2*zeta-1)/(-1+zeta)^3+6*p2*p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^4
# phi_zetazeta[3] = -1/2*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+p4*(1-2*zeta)/(-1+zeta)^2-2*p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+p3*(1-2*zeta)/(-1+zeta)^2-2*p3*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+2*p3*p4/(-1+zeta)^2+4*p3*p4*(1-2*zeta)/(-1+zeta)^3-6*p3*p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^4
# phi_zetazeta[4] = 4
# phi_zetazeta[5] = -2*p1*eta/(-1+zeta)^2-2*p4*eta/(-1+zeta)^2-8*p4*p1*eta/(-1+zeta)^3-2*p2*eta/(-1+zeta)^2-8*p2*p4*eta/(-1+zeta)^3-8*p1*p2*eta/(-1+zeta)^3-24*p2*p1*p4*eta/(-1+zeta)^4
# phi_zetazeta[6] = 2*p3*xi/(-1+zeta)^2+2*p2*xi/(-1+zeta)^2+8*p2*p3*xi/(-1+zeta)^3+2*p1*xi/(-1+zeta)^2+8*p1*p3*xi/(-1+zeta)^3+8*p1*p2*xi/(-1+zeta)^3+24*p1*p2*p3*xi/(-1+zeta)^4
# phi_zetazeta[7] = 2*p4*eta/(-1+zeta)^2+2*p3*eta/(-1+zeta)^2+8*p3*p4*eta/(-1+zeta)^3+2*p2*eta/(-1+zeta)^2+8*p2*p4*eta/(-1+zeta)^3+8*p2*p3*eta/(-1+zeta)^3+24*p2*p3*p4*eta/(-1+zeta)^4
# phi_zetazeta[8] = -2*p1*xi/(-1+zeta)^2-2*p4*xi/(-1+zeta)^2-8*p4*p1*xi/(-1+zeta)^3-2*p3*xi/(-1+zeta)^2-8*p1*p3*xi/(-1+zeta)^3-8*p3*p4*xi/(-1+zeta)^3-24*p3*p4*p1*xi/(-1+zeta)^4
# phi_zetazeta[9] = -2*zeta/(-1+zeta)+4*p4/(-1+zeta)-4*p4*zeta/(-1+zeta)^2+4*p1/(-1+zeta)-4*p1*zeta/(-1+zeta)^2+8*p4*p1/(-1+zeta)^2-8*p1*p4*zeta/(-1+zeta)^3
# phi_zetazeta[10] = -2*zeta/(-1+zeta)+4*p1/(-1+zeta)-4*p1*zeta/(-1+zeta)^2+4*p2/(-1+zeta)-4*p2*zeta/(-1+zeta)^2+8*p1*p2/(-1+zeta)^2-8*p2*p1*zeta/(-1+zeta)^3
# phi_zetazeta[11] = -2*zeta/(-1+zeta)+4*p2/(-1+zeta)-4*p2*zeta/(-1+zeta)^2+4*p3/(-1+zeta)-4*p3*zeta/(-1+zeta)^2+8*p2*p3/(-1+zeta)^2-8*p3*p2*zeta/(-1+zeta)^3
# phi_zetazeta[12] = -2*zeta/(-1+zeta)+4*p3/(-1+zeta)-4*p3*zeta/(-1+zeta)^2+4*p4/(-1+zeta)-4*p4*zeta/(-1+zeta)^2+8*p3*p4/(-1+zeta)^2-8*p4*p3*zeta/(-1+zeta)^3
# phi_zetazeta[13] = 8*p3*p4/(-1+zeta)^2+8*p2*p4/(-1+zeta)^2+8*p2*p3/(-1+zeta)^2+32*p2*p3*p4/(-1+zeta)^3+8*p4*p1/(-1+zeta)^2+8*p1*p3/(-1+zeta)^2+32*p3*p4*p1/(-1+zeta)^3+8*p1*p2/(-1+zeta)^2+32*p2*p1*p4/(-1+zeta)^3+32*p1*p2*p3/(-1+zeta)^3+96*p1*p2*p3*p4/(-1+zeta)^4

# phi_xieta[0] = 1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-1/2*p1*xi/(-1+zeta)^2-1/2*p4*eta/(-1+zeta)^2+p4*p1/(-1+zeta)^2
# phi_xieta[1] = 1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2-1/2*p1*xi/(-1+zeta)^2+1/2*p2*eta/(-1+zeta)^2-p1*p2/(-1+zeta)^2
# phi_xieta[2] = 1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2+1/2*p3*xi/(-1+zeta)^2+1/2*p2*eta/(-1+zeta)^2+p2*p3/(-1+zeta)^2
# phi_xieta[3] = 1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+1/2*p3*xi/(-1+zeta)^2-1/2*p4*eta/(-1+zeta)^2-p3*p4/(-1+zeta)^2
# phi_xieta[4] = 0
# phi_xieta[5] = p4*eta/(-1+zeta)^2-2*p4*p1/(-1+zeta)^2-p2*eta/(-1+zeta)^2+2*p1*p2/(-1+zeta)^2
# phi_xieta[6] = -p3*xi/(-1+zeta)^2+p1*xi/(-1+zeta)^2-2*p2*p3/(-1+zeta)^2+2*p1*p2/(-1+zeta)^2
# phi_xieta[7] = p4*eta/(-1+zeta)^2+2*p3*p4/(-1+zeta)^2-p2*eta/(-1+zeta)^2-2*p2*p3/(-1+zeta)^2
# phi_xieta[8] = -p3*xi/(-1+zeta)^2+p1*xi/(-1+zeta)^2-2*p4*p1/(-1+zeta)^2+2*p3*p4/(-1+zeta)^2
# 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_xieta[13] = 4*p4*p1/(-1+zeta)^2-4*p3*p4/(-1+zeta)^2+4*p2*p3/(-1+zeta)^2-4*p1*p2/(-1+zeta)^2

# phi_xizeta[0] = 1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-1/2*p1*(2*zeta-1)/(-1+zeta)^2+p1*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-1/2*p1*eta/(-1+zeta)^2-1/2*p4*eta/(-1+zeta)^2-2*p4*p1*eta/(-1+zeta)^3
# phi_xizeta[1] = 1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2-1/2*p1*(1-2*zeta)/(-1+zeta)^2+p1*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+1/2*p2*eta/(-1+zeta)^2+1/2*p1*eta/(-1+zeta)^2+2*p1*p2*eta/(-1+zeta)^3
# phi_xizeta[2] = -1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2+1/2*p3*(2*zeta-1)/(-1+zeta)^2-p3*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-1/2*p3*eta/(-1+zeta)^2-1/2*p2*eta/(-1+zeta)^2-2*p2*p3*eta/(-1+zeta)^3
# phi_xizeta[3] = -1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+1/2*p3*(1-2*zeta)/(-1+zeta)^2-p3*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+1/2*p4*eta/(-1+zeta)^2+1/2*p3*eta/(-1+zeta)^2+2*p3*p4*eta/(-1+zeta)^3
# phi_xizeta[4] = 0
# phi_xizeta[5] = p4*eta/(-1+zeta)^2+4*p4*p1*eta/(-1+zeta)^3-p2*eta/(-1+zeta)^2-4*p1*p2*eta/(-1+zeta)^3
# phi_xizeta[6] = -p3*xi/(-1+zeta)^2-p1*xi/(-1+zeta)^2-4*p1*p3*xi/(-1+zeta)^3-2*p2*p3/(-1+zeta)^2-2*p1*p3/(-1+zeta)^2-2*p1*p2/(-1+zeta)^2-8*p1*p2*p3/(-1+zeta)^3
# phi_xizeta[7] = -p4*eta/(-1+zeta)^2-4*p3*p4*eta/(-1+zeta)^3+p2*eta/(-1+zeta)^2+4*p2*p3*eta/(-1+zeta)^3
# phi_xizeta[8] = -p3*xi/(-1+zeta)^2-p1*xi/(-1+zeta)^2-4*p1*p3*xi/(-1+zeta)^3+2*p4*p1/(-1+zeta)^2+2*p1*p3/(-1+zeta)^2+2*p3*p4/(-1+zeta)^2+8*p3*p4*p1/(-1+zeta)^3
# phi_xizeta[9] = -zeta/(-1+zeta)+2*p1/(-1+zeta)-2*p1*zeta/(-1+zeta)^2
# phi_xizeta[10] = zeta/(-1+zeta)-2*p1/(-1+zeta)+2*p1*zeta/(-1+zeta)^2
# phi_xizeta[11] = zeta/(-1+zeta)-2*p3/(-1+zeta)+2*p3*zeta/(-1+zeta)^2
# phi_xizeta[12] = -zeta/(-1+zeta)+2*p3/(-1+zeta)-2*p3*zeta/(-1+zeta)^2
# phi_xizeta[13] = -4*p4*p1/(-1+zeta)^2-4*p3*p4/(-1+zeta)^2-16*p3*p4*p1/(-1+zeta)^3+4*p2*p3/(-1+zeta)^2+4*p1*p2/(-1+zeta)^2+16*p1*p2*p3/(-1+zeta)^3

# phi_etazeta[0] = 1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2-1/2*p4*(2*zeta-1)/(-1+zeta)^2+p4*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-1/2*p1*xi/(-1+zeta)^2-1/2*p4*xi/(-1+zeta)^2-2*p4*p1*xi/(-1+zeta)^3
# phi_etazeta[1] = -1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2+1/2*p2*(1-2*zeta)/(-1+zeta)^2-p2*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+1/2*p2*xi/(-1+zeta)^2+1/2*p1*xi/(-1+zeta)^2+2*p1*p2*xi/(-1+zeta)^3
# phi_etazeta[2] = -1/4*(xi*eta-zeta+zeta^2)/(-1+zeta)^2+1/2*p2*(2*zeta-1)/(-1+zeta)^2-p2*(xi*eta-zeta+zeta^2)/(-1+zeta)^3-1/2*p3*xi/(-1+zeta)^2-1/2*p2*xi/(-1+zeta)^2-2*p2*p3*xi/(-1+zeta)^3
# phi_etazeta[3] = 1/4*(xi*eta+zeta-zeta^2)/(-1+zeta)^2-1/2*p4*(1-2*zeta)/(-1+zeta)^2+p4*(xi*eta+zeta-zeta^2)/(-1+zeta)^3+1/2*p4*xi/(-1+zeta)^2+1/2*p3*xi/(-1+zeta)^2+2*p3*p4*xi/(-1+zeta)^3
# phi_etazeta[4] = 0
# phi_etazeta[5] = -p4*eta/(-1+zeta)^2-p2*eta/(-1+zeta)^2-4*p2*p4*eta/(-1+zeta)^3+2*p4*p1/(-1+zeta)^2+2*p2*p4/(-1+zeta)^2+2*p1*p2/(-1+zeta)^2+8*p2*p1*p4/(-1+zeta)^3
# phi_etazeta[6] = p3*xi/(-1+zeta)^2+4*p2*p3*xi/(-1+zeta)^3-p1*xi/(-1+zeta)^2-4*p1*p2*xi/(-1+zeta)^3
# phi_etazeta[7] = -p4*eta/(-1+zeta)^2-p2*eta/(-1+zeta)^2-4*p2*p4*eta/(-1+zeta)^3-2*p3*p4/(-1+zeta)^2-2*p2*p4/(-1+zeta)^2-2*p2*p3/(-1+zeta)^2-8*p2*p3*p4/(-1+zeta)^3
# phi_etazeta[8] = p1*xi/(-1+zeta)^2+4*p4*p1*xi/(-1+zeta)^3-p3*xi/(-1+zeta)^2-4*p3*p4*xi/(-1+zeta)^3
# phi_etazeta[9] = -zeta/(-1+zeta)+2*p4/(-1+zeta)-2*p4*zeta/(-1+zeta)^2
# phi_etazeta[10] = -zeta/(-1+zeta)+2*p2/(-1+zeta)-2*p2*zeta/(-1+zeta)^2
# phi_etazeta[11] = zeta/(-1+zeta)-2*p2/(-1+zeta)+2*p2*zeta/(-1+zeta)^2
# phi_etazeta[12] = zeta/(-1+zeta)-2*p4/(-1+zeta)+2*p4*zeta/(-1+zeta)^2
# phi_etazeta[13] = 4*p3*p4/(-1+zeta)^2+4*p2*p3/(-1+zeta)^2+16*p2*p3*p4/(-1+zeta)^3-4*p4*p1/(-1+zeta)^2-4*p1*p2/(-1+zeta)^2-16*p2*p1*p4/(-1+zeta)^3


# Results with nonzero epsilon
#
# phi[0] = p1*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2
# phi[1] = p2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2
# phi[2] = p3*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2
# phi[3] = p4*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2
# phi[4] = zeta*(2*zeta-1)
# phi[5] = -4*p4*p1*p2*(eta+eps)/(1-zeta+eps)^2
# phi[6] = -4*p1*p2*p3*(-xi+eps)/(1-zeta+eps)^2
# phi[7] = -4*p2*p3*p4*(-eta+eps)/(1-zeta+eps)^2
# phi[8] = -4*p3*p4*p1*(xi+eps)/(1-zeta+eps)^2
# phi[9] = 4*p1*p4*zeta/(1-zeta+eps)
# phi[10] = 4*p2*p1*zeta/(1-zeta+eps)
# phi[11] = 4*p3*p2*zeta/(1-zeta+eps)
# phi[12] = 4*p4*p3*zeta/(1-zeta+eps)
# phi[13] = 16*p1*p2*p3*p4/(1-zeta+eps)^2

# phi_xi[0] = -1/2*p1*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps-2*p4*eta-2*p4*eps)/(1-zeta+eps)^2
# phi_xi[1] = 1/2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps-2*p2*eta-2*p2*eps)/(1-zeta+eps)^2
# phi_xi[2] = 1/2*p3*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps+2*p2*eta-2*p2*eps)/(1-zeta+eps)^2
# phi_xi[3] = -1/2*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps+2*p4*eta-2*p4*eps)/(1-zeta+eps)^2
# phi_xi[4] = 0
# phi_xi[5] = 2*p1*(eta+eps)*xi/(1-zeta+eps)^2
# phi_xi[6] = -2*p1*p3*(-xi+eps-2*p2)/(1-zeta+eps)^2
# phi_xi[7] = 2*p3*(-eta+eps)*xi/(1-zeta+eps)^2
# phi_xi[8] = 2*p1*p3*(xi+eps-2*p4)/(1-zeta+eps)^2
# phi_xi[9] = -2*p1*zeta/(1-zeta+eps)
# phi_xi[10] = 2*p1*zeta/(1-zeta+eps)
# phi_xi[11] = 2*p3*zeta/(1-zeta+eps)
# phi_xi[12] = -2*p3*zeta/(1-zeta+eps)
# phi_xi[13] = -8*p1*p3*xi/(1-zeta+eps)^2

# phi_eta[0] = -1/2*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps-2*p1*xi-2*p1*eps)/(1-zeta+eps)^2
# phi_eta[1] = -1/2*p2*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps+2*p1*xi-2*p1*eps)/(1-zeta+eps)^2
# phi_eta[2] = 1/2*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps+2*p3*xi-2*p3*eps)/(1-zeta+eps)^2
# phi_eta[3] = 1/2*p4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps-2*p3*xi-2*p3*eps)/(1-zeta+eps)^2
# phi_eta[4] = 0
# phi_eta[5] = 2*p4*p2*(eta+eps-2*p1)/(1-zeta+eps)^2
# phi_eta[6] = 2*p2*(-xi+eps)*eta/(1-zeta+eps)^2
# phi_eta[7] = -2*p4*p2*(-eta+eps-2*p3)/(1-zeta+eps)^2
# phi_eta[8] = 2*p4*(xi+eps)*eta/(1-zeta+eps)^2
# phi_eta[9] = -2*p4*zeta/(1-zeta+eps)
# phi_eta[10] = -2*p2*zeta/(1-zeta+eps)
# phi_eta[11] = 2*p2*zeta/(1-zeta+eps)
# phi_eta[12] = 2*p4*zeta/(1-zeta+eps)
# phi_eta[13] = -8*p4*p2*eta/(1-zeta+eps)^2

# phi_zeta[0] = -1/2*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p1*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+p1*p4*(-1+2*zeta-eps)/(1-zeta+eps)^2+2*p1*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3
# phi_zeta[1] = -1/2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p2*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+p2*p1*(-1+2*zeta-eps)/(1-zeta+eps)^2+2*p2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3
# phi_zeta[2] = -1/2*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p3*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+p3*p2*(-1+2*zeta-eps)/(1-zeta+eps)^2+2*p3*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3
# phi_zeta[3] = -1/2*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+p4*p3*(-1+2*zeta-eps)/(1-zeta+eps)^2+2*p4*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3
# phi_zeta[4] = 4*zeta-1
# phi_zeta[5] = 2*p1*p2*(eta+eps)/(1-zeta+eps)^2+2*p4*p2*(eta+eps)/(1-zeta+eps)^2+2*p1*p4*(eta+eps)/(1-zeta+eps)^2-8*p4*p1*p2*(eta+eps)/(1-zeta+eps)^3
# phi_zeta[6] = 2*p2*p3*(-xi+eps)/(1-zeta+eps)^2+2*p1*p3*(-xi+eps)/(1-zeta+eps)^2+2*p2*p1*(-xi+eps)/(1-zeta+eps)^2-8*p1*p2*p3*(-xi+eps)/(1-zeta+eps)^3
# phi_zeta[7] = 2*p4*p3*(-eta+eps)/(1-zeta+eps)^2+2*p2*p4*(-eta+eps)/(1-zeta+eps)^2+2*p2*p3*(-eta+eps)/(1-zeta+eps)^2-8*p2*p3*p4*(-eta+eps)/(1-zeta+eps)^3
# phi_zeta[8] = 2*p1*p4*(xi+eps)/(1-zeta+eps)^2+2*p3*p1*(xi+eps)/(1-zeta+eps)^2+2*p3*p4*(xi+eps)/(1-zeta+eps)^2-8*p3*p4*p1*(xi+eps)/(1-zeta+eps)^3
# phi_zeta[9] = -2*p4*zeta/(1-zeta+eps)-2*p1*zeta/(1-zeta+eps)+4*p1*p4/(1-zeta+eps)+4*p1*p4*zeta/(1-zeta+eps)^2
# phi_zeta[10] = -2*p1*zeta/(1-zeta+eps)-2*p2*zeta/(1-zeta+eps)+4*p2*p1/(1-zeta+eps)+4*p2*p1*zeta/(1-zeta+eps)^2
# phi_zeta[11] = -2*p2*zeta/(1-zeta+eps)-2*p3*zeta/(1-zeta+eps)+4*p3*p2/(1-zeta+eps)+4*p3*p2*zeta/(1-zeta+eps)^2
# phi_zeta[12] = -2*p3*zeta/(1-zeta+eps)-2*p4*zeta/(1-zeta+eps)+4*p4*p3/(1-zeta+eps)+4*p4*p3*zeta/(1-zeta+eps)^2
# phi_zeta[13] = -8*p2*p3*p4/(1-zeta+eps)^2-8*p3*p4*p1/(1-zeta+eps)^2-8*p4*p1*p2/(1-zeta+eps)^2-8*p1*p2*p3/(1-zeta+eps)^2+32*p1*p2*p3*p4/(1-zeta+eps)^3

# phi_xixi[0] = -p1*(eta+eps)/(1-zeta+eps)^2
# phi_xixi[1] = p1*(-eta-eps)/(1-zeta+eps)^2
# phi_xixi[2] = p3*(eta-eps)/(1-zeta+eps)^2
# phi_xixi[3] = -p3*(-eta+eps)/(1-zeta+eps)^2
# phi_xixi[4] = 0
# phi_xixi[5] = 2*p1*(eta+eps)/(1-zeta+eps)^2
# phi_xixi[6] = 4*p1*p3/(1-zeta+eps)^2
# phi_xixi[7] = 2*p3*(-eta+eps)/(1-zeta+eps)^2
# phi_xixi[8] = 4*p1*p3/(1-zeta+eps)^2
# phi_xixi[9] = 0
# phi_xixi[10] = 0
# phi_xixi[11] = 0
# phi_xixi[12] = 0
# phi_xixi[13] = -8*p1*p3/(1-zeta+eps)^2

# phi_etaeta[0] = -p4*(xi+eps)/(1-zeta+eps)^2
# phi_etaeta[1] = -p2*(-xi+eps)/(1-zeta+eps)^2
# phi_etaeta[2] = p2*(xi-eps)/(1-zeta+eps)^2
# phi_etaeta[3] = p4*(-xi-eps)/(1-zeta+eps)^2
# phi_etaeta[4] = 0
# phi_etaeta[5] = 4*p4*p2/(1-zeta+eps)^2
# phi_etaeta[6] = 2*p2*(-xi+eps)/(1-zeta+eps)^2
# phi_etaeta[7] = 4*p4*p2/(1-zeta+eps)^2
# phi_etaeta[8] = 2*p4*(xi+eps)/(1-zeta+eps)^2
# phi_etaeta[9] = 0
# phi_etaeta[10] = 0
# phi_etaeta[11] = 0
# phi_etaeta[12] = 0
# phi_etaeta[13] = -8*p4*p2/(1-zeta+eps)^2

# phi_zetazeta[0] = 1/2*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-p4*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-p1*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p1*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3+2*p1*p4/(1-zeta+eps)^2+4*p1*p4*(-1+2*zeta-eps)/(1-zeta+eps)^3+6*p1*p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^4
# phi_zetazeta[1] = 1/2*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-p1*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-p2*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p2*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3+2*p2*p1/(1-zeta+eps)^2+4*p2*p1*(-1+2*zeta-eps)/(1-zeta+eps)^3+6*p2*p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^4
# phi_zetazeta[2] = 1/2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-p2*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-p3*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p3*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3+2*p3*p2/(1-zeta+eps)^2+4*p3*p2*(-1+2*zeta-eps)/(1-zeta+eps)^3+6*p3*p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^4
# phi_zetazeta[3] = 1/2*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-p3*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-p4*(-1+2*zeta-eps)/(1-zeta+eps)^2-2*p4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3+2*p4*p3/(1-zeta+eps)^2+4*p4*p3*(-1+2*zeta-eps)/(1-zeta+eps)^3+6*p4*p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^4
# phi_zetazeta[4] = 4
# phi_zetazeta[5] = -2*p2*(eta+eps)/(1-zeta+eps)^2-2*p1*(eta+eps)/(1-zeta+eps)^2+8*p1*p2*(eta+eps)/(1-zeta+eps)^3-2*p4*(eta+eps)/(1-zeta+eps)^2+8*p4*p2*(eta+eps)/(1-zeta+eps)^3+8*p1*p4*(eta+eps)/(1-zeta+eps)^3-24*p4*p1*p2*(eta+eps)/(1-zeta+eps)^4
# phi_zetazeta[6] = -2*p3*(-xi+eps)/(1-zeta+eps)^2-2*p2*(-xi+eps)/(1-zeta+eps)^2+8*p2*p3*(-xi+eps)/(1-zeta+eps)^3-2*p1*(-xi+eps)/(1-zeta+eps)^2+8*p1*p3*(-xi+eps)/(1-zeta+eps)^3+8*p2*p1*(-xi+eps)/(1-zeta+eps)^3-24*p1*p2*p3*(-xi+eps)/(1-zeta+eps)^4
# phi_zetazeta[7] = -2*p3*(-eta+eps)/(1-zeta+eps)^2-2*p4*(-eta+eps)/(1-zeta+eps)^2+8*p4*p3*(-eta+eps)/(1-zeta+eps)^3-2*p2*(-eta+eps)/(1-zeta+eps)^2+8*p2*p4*(-eta+eps)/(1-zeta+eps)^3+8*p2*p3*(-eta+eps)/(1-zeta+eps)^3-24*p2*p3*p4*(-eta+eps)/(1-zeta+eps)^4
# phi_zetazeta[8] = -2*p4*(xi+eps)/(1-zeta+eps)^2-2*p1*(xi+eps)/(1-zeta+eps)^2+8*p1*p4*(xi+eps)/(1-zeta+eps)^3-2*p3*(xi+eps)/(1-zeta+eps)^2+8*p3*p1*(xi+eps)/(1-zeta+eps)^3+8*p3*p4*(xi+eps)/(1-zeta+eps)^3-24*p3*p4*p1*(xi+eps)/(1-zeta+eps)^4
# phi_zetazeta[9] = 2*zeta/(1-zeta+eps)-4*p4/(1-zeta+eps)-4*p4*zeta/(1-zeta+eps)^2-4*p1/(1-zeta+eps)-4*p1*zeta/(1-zeta+eps)^2+8*p1*p4/(1-zeta+eps)^2+8*p1*p4*zeta/(1-zeta+eps)^3
# phi_zetazeta[10] = 2*zeta/(1-zeta+eps)-4*p1/(1-zeta+eps)-4*p1*zeta/(1-zeta+eps)^2-4*p2/(1-zeta+eps)-4*p2*zeta/(1-zeta+eps)^2+8*p2*p1/(1-zeta+eps)^2+8*p2*p1*zeta/(1-zeta+eps)^3
# phi_zetazeta[11] = 2*zeta/(1-zeta+eps)-4*p2/(1-zeta+eps)-4*p2*zeta/(1-zeta+eps)^2-4*p3/(1-zeta+eps)-4*p3*zeta/(1-zeta+eps)^2+8*p3*p2/(1-zeta+eps)^2+8*p3*p2*zeta/(1-zeta+eps)^3
# phi_zetazeta[12] = 2*zeta/(1-zeta+eps)-4*p3/(1-zeta+eps)-4*p3*zeta/(1-zeta+eps)^2-4*p4/(1-zeta+eps)-4*p4*zeta/(1-zeta+eps)^2+8*p4*p3/(1-zeta+eps)^2+8*p4*p3*zeta/(1-zeta+eps)^3
# phi_zetazeta[13] = 8*p4*p3/(1-zeta+eps)^2+8*p4*p2/(1-zeta+eps)^2+8*p3*p2/(1-zeta+eps)^2-32*p2*p3*p4/(1-zeta+eps)^3+8*p1*p4/(1-zeta+eps)^2+8*p1*p3/(1-zeta+eps)^2-32*p3*p4*p1/(1-zeta+eps)^3+8*p2*p1/(1-zeta+eps)^2-32*p4*p1*p2/(1-zeta+eps)^3-32*p1*p2*p3/(1-zeta+eps)^3+96*p1*p2*p3*p4/(1-zeta+eps)^4

# phi_xieta[0] = 1/4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p1*(xi+eps)/(1-zeta+eps)^2-1/2*p4*(eta+eps)/(1-zeta+eps)^2+p1*p4/(1-zeta+eps)^2
# phi_xieta[1] = -1/4*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p1*(-xi+eps)/(1-zeta+eps)^2-1/2*p2*(-eta-eps)/(1-zeta+eps)^2-p2*p1/(1-zeta+eps)^2
# phi_xieta[2] = 1/4*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p3*(xi-eps)/(1-zeta+eps)^2+1/2*p2*(eta-eps)/(1-zeta+eps)^2+p3*p2/(1-zeta+eps)^2
# phi_xieta[3] = -1/4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p3*(-xi-eps)/(1-zeta+eps)^2+1/2*p4*(-eta+eps)/(1-zeta+eps)^2-p4*p3/(1-zeta+eps)^2
# phi_xieta[4] = 0
# phi_xieta[5] = -p2*(eta+eps)/(1-zeta+eps)^2+2*p2*p1/(1-zeta+eps)^2+p4*(eta+eps)/(1-zeta+eps)^2-2*p1*p4/(1-zeta+eps)^2
# phi_xieta[6] = p3*(-xi+eps)/(1-zeta+eps)^2-p1*(-xi+eps)/(1-zeta+eps)^2-2*p3*p2/(1-zeta+eps)^2+2*p2*p1/(1-zeta+eps)^2
# phi_xieta[7] = -p4*(-eta+eps)/(1-zeta+eps)^2+2*p4*p3/(1-zeta+eps)^2+p2*(-eta+eps)/(1-zeta+eps)^2-2*p3*p2/(1-zeta+eps)^2
# phi_xieta[8] = p1*(xi+eps)/(1-zeta+eps)^2-p3*(xi+eps)/(1-zeta+eps)^2-2*p1*p4/(1-zeta+eps)^2+2*p4*p3/(1-zeta+eps)^2
# phi_xieta[9] = zeta/(1-zeta+eps)
# phi_xieta[10] = -zeta/(1-zeta+eps)
# phi_xieta[11] = zeta/(1-zeta+eps)
# phi_xieta[12] = -zeta/(1-zeta+eps)
# phi_xieta[13] = 4*p1*p4/(1-zeta+eps)^2-4*p4*p3/(1-zeta+eps)^2+4*p3*p2/(1-zeta+eps)^2-4*p2*p1/(1-zeta+eps)^2

# phi_xizeta[0] = 1/4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p1*(-1+2*zeta-eps)/(1-zeta+eps)^2-p1*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p4*(eta+eps)/(1-zeta+eps)^2-1/2*p1*(eta+eps)/(1-zeta+eps)^2+2*p1*p4*(eta+eps)/(1-zeta+eps)^3
# phi_xizeta[1] = -1/4*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p1*(-1+2*zeta-eps)/(1-zeta+eps)^2+p1*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p1*(-eta-eps)/(1-zeta+eps)^2-1/2*p2*(-eta-eps)/(1-zeta+eps)^2+2*p2*p1*(-eta-eps)/(1-zeta+eps)^3
# phi_xizeta[2] = -1/4*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p3*(-1+2*zeta-eps)/(1-zeta+eps)^2+p3*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p2*(eta-eps)/(1-zeta+eps)^2-1/2*p3*(eta-eps)/(1-zeta+eps)^2+2*p3*p2*(eta-eps)/(1-zeta+eps)^3
# phi_xizeta[3] = 1/4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p3*(-1+2*zeta-eps)/(1-zeta+eps)^2-p3*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p3*(-eta+eps)/(1-zeta+eps)^2-1/2*p4*(-eta+eps)/(1-zeta+eps)^2+2*p4*p3*(-eta+eps)/(1-zeta+eps)^3
# phi_xizeta[4] = 0
# phi_xizeta[5] = -p2*(eta+eps)/(1-zeta+eps)^2+4*p1*p2*(eta+eps)/(1-zeta+eps)^3+p4*(eta+eps)/(1-zeta+eps)^2-4*p1*p4*(eta+eps)/(1-zeta+eps)^3
# phi_xizeta[6] = p3*(-xi+eps)/(1-zeta+eps)^2+p1*(-xi+eps)/(1-zeta+eps)^2-4*p1*p3*(-xi+eps)/(1-zeta+eps)^3-2*p3*p2/(1-zeta+eps)^2-2*p1*p3/(1-zeta+eps)^2-2*p2*p1/(1-zeta+eps)^2+8*p1*p2*p3/(1-zeta+eps)^3
# phi_xizeta[7] = p4*(-eta+eps)/(1-zeta+eps)^2-4*p4*p3*(-eta+eps)/(1-zeta+eps)^3-p2*(-eta+eps)/(1-zeta+eps)^2+4*p2*p3*(-eta+eps)/(1-zeta+eps)^3
# phi_xizeta[8] = -p1*(xi+eps)/(1-zeta+eps)^2-p3*(xi+eps)/(1-zeta+eps)^2+4*p3*p1*(xi+eps)/(1-zeta+eps)^3+2*p1*p4/(1-zeta+eps)^2+2*p1*p3/(1-zeta+eps)^2+2*p4*p3/(1-zeta+eps)^2-8*p3*p4*p1/(1-zeta+eps)^3
# phi_xizeta[9] = zeta/(1-zeta+eps)-2*p1/(1-zeta+eps)-2*p1*zeta/(1-zeta+eps)^2
# phi_xizeta[10] = -zeta/(1-zeta+eps)+2*p1/(1-zeta+eps)+2*p1*zeta/(1-zeta+eps)^2
# phi_xizeta[11] = -zeta/(1-zeta+eps)+2*p3/(1-zeta+eps)+2*p3*zeta/(1-zeta+eps)^2
# phi_xizeta[12] = zeta/(1-zeta+eps)-2*p3/(1-zeta+eps)-2*p3*zeta/(1-zeta+eps)^2
# phi_xizeta[13] = -4*p1*p4/(1-zeta+eps)^2-4*p4*p3/(1-zeta+eps)^2+16*p3*p4*p1/(1-zeta+eps)^3+4*p3*p2/(1-zeta+eps)^2+4*p2*p1/(1-zeta+eps)^2-16*p1*p2*p3/(1-zeta+eps)^3

# phi_etazeta[0] = 1/4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p4*(-1+2*zeta-eps)/(1-zeta+eps)^2-p4*(eta*xi+eta*eps+eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p4*(xi+eps)/(1-zeta+eps)^2-1/2*p1*(xi+eps)/(1-zeta+eps)^2+2*p1*p4*(xi+eps)/(1-zeta+eps)^3
# phi_etazeta[1] = 1/4*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2-1/2*p2*(-1+2*zeta-eps)/(1-zeta+eps)^2-p2*(-eta*xi-eps*xi+eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p1*(-xi+eps)/(1-zeta+eps)^2-1/2*p2*(-xi+eps)/(1-zeta+eps)^2+2*p2*p1*(-xi+eps)/(1-zeta+eps)^3
# phi_etazeta[2] = -1/4*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p2*(-1+2*zeta-eps)/(1-zeta+eps)^2+p2*(eta*xi-eta*eps-eps*xi+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p2*(xi-eps)/(1-zeta+eps)^2-1/2*p3*(xi-eps)/(1-zeta+eps)^2+2*p3*p2*(xi-eps)/(1-zeta+eps)^3
# phi_etazeta[3] = -1/4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^2+1/2*p4*(-1+2*zeta-eps)/(1-zeta+eps)^2+p4*(-eta*xi+eps*xi-eta*eps+eps^2-zeta+zeta^2-zeta*eps)/(1-zeta+eps)^3-1/2*p3*(-xi-eps)/(1-zeta+eps)^2-1/2*p4*(-xi-eps)/(1-zeta+eps)^2+2*p4*p3*(-xi-eps)/(1-zeta+eps)^3
# phi_etazeta[4] = 0
# phi_etazeta[5] = -p2*(eta+eps)/(1-zeta+eps)^2-p4*(eta+eps)/(1-zeta+eps)^2+4*p4*p2*(eta+eps)/(1-zeta+eps)^3+2*p2*p1/(1-zeta+eps)^2+2*p4*p2/(1-zeta+eps)^2+2*p1*p4/(1-zeta+eps)^2-8*p4*p1*p2/(1-zeta+eps)^3
# phi_etazeta[6] = -p3*(-xi+eps)/(1-zeta+eps)^2+4*p2*p3*(-xi+eps)/(1-zeta+eps)^3+p1*(-xi+eps)/(1-zeta+eps)^2-4*p2*p1*(-xi+eps)/(1-zeta+eps)^3
# phi_etazeta[7] = p4*(-eta+eps)/(1-zeta+eps)^2+p2*(-eta+eps)/(1-zeta+eps)^2-4*p2*p4*(-eta+eps)/(1-zeta+eps)^3-2*p4*p3/(1-zeta+eps)^2-2*p4*p2/(1-zeta+eps)^2-2*p3*p2/(1-zeta+eps)^2+8*p2*p3*p4/(1-zeta+eps)^3
# phi_etazeta[8] = p1*(xi+eps)/(1-zeta+eps)^2-4*p1*p4*(xi+eps)/(1-zeta+eps)^3-p3*(xi+eps)/(1-zeta+eps)^2+4*p3*p4*(xi+eps)/(1-zeta+eps)^3
# phi_etazeta[9] = zeta/(1-zeta+eps)-2*p4/(1-zeta+eps)-2*p4*zeta/(1-zeta+eps)^2
# phi_etazeta[10] = zeta/(1-zeta+eps)-2*p2/(1-zeta+eps)-2*p2*zeta/(1-zeta+eps)^2
# phi_etazeta[11] = -zeta/(1-zeta+eps)+2*p2/(1-zeta+eps)+2*p2*zeta/(1-zeta+eps)^2
# phi_etazeta[12] = -zeta/(1-zeta+eps)+2*p4/(1-zeta+eps)+2*p4*zeta/(1-zeta+eps)^2
# phi_etazeta[13] = 4*p4*p3/(1-zeta+eps)^2+4*p3*p2/(1-zeta+eps)^2-16*p2*p3*p4/(1-zeta+eps)^3-4*p2*p1/(1-zeta+eps)^2-4*p1*p4/(1-zeta+eps)^2+16*p4*p1*p2/(1-zeta+eps)^3