Instantly share code, notes, and snippets.

# jwpeterson/pyramid13.maple Created Dec 8, 2013

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