Given basis functions for the pyramid 5, compute first and second derivatives and check various qualities that each should have.

pyramid5.maple

# Shape functions, first, and second derivatives of PYRAMID5 basis functions using Maple

pyramid5 := proc()

uses LinearAlgebra;

local
 phi, i, t, phi_xi, phi_eta, phi_zeta, phi_xixi, phi_etaeta, phi_zetazeta,
 phi_xieta, phi_xizeta, phi_etazeta, x, n, xi_vec, planes, phi_sum, 
 phi_xi_sum, phi_eta_sum, phi_zeta_sum;

# Stop reading when any type of error is encountered
interface(errorbreak=2):

phi := Vector(5);

phi[0+1] := (1/4)*(zeta + xi - 1)*(zeta + eta - 1)/((1 - zeta) + eps);
phi[1+1] := (1/4)*(zeta - xi - 1)*(zeta + eta - 1)/((1 - zeta) + eps);
phi[2+1] := (1/4)*(zeta - xi - 1)*(zeta - eta - 1)/((1 - zeta) + eps);
phi[3+1] := (1/4)*(zeta + xi - 1)*(zeta - eta - 1)/((1 - zeta) + eps);
phi[4+1] := zeta;

# Maple does not seem to be able to compute limits when there are multiple
# indeterminates present...
# t := limit(phi[1], zeta=1);
# printf("t=%a\n", t);

# t := limit(phi[1], {xi=0, eta=0, zeta=1});
# printf("t=%a\n", t);

# But it can do it if you substitute out xi and eta first
for i from 1 to 4 do
  t := limit(subs(xi=0, eta=0, phi[i]), zeta=1);
  printf("limit(phi[%d]) as zeta->1 = %a\n", i-1, t);
end do;

printf("\n");

# Simple evaluation does result in division by zero, though,
# which triggers a Maple error...
# for i from 1 to 4 do
#   t := subs(xi=0, eta=0, zeta=1, phi[i]);
#   printf("Evaluating phi[%d] at (0,0,1): %a\n", i-1, t);
# end do;

# Check the "partition of unity" property of these basis functions
printf("Checking partition of unity property of basis functions...");
phi_sum := 0;
for i from 1 to 5 do
  phi_sum := phi_sum + phi[i];
end do;

phi_sum := simplify(phi_sum);
phi_sum := subs(eps=0, phi_sum);

if (phi_sum <> 1) then
  printf("phi_sum = %a\n", phi_sum);
  error;
end if;

printf("success!\n");



# Compute first derivatives
phi_xi := Vector(5);
phi_eta := Vector(5);
phi_zeta := Vector(5);

for i from 1 to 5 do
  phi_xi[i] := diff(phi[i], xi);
  phi_eta[i] := diff(phi[i], eta);
  phi_zeta[i] := diff(phi[i], zeta);
end do;

# The phi_zeta derivatives are pretty long... the standard Maple
# simplify command works pretty well on them though.
for i from 1 to 5 do
  phi_zeta[i] := simplify(phi_zeta[i]);

  # Let eps=0 in the numerator of each phi_zeta term
  phi_zeta[i] := subs(eps=0, numer(phi_zeta[i])) / denom(phi_zeta[i]);
end do;

printf("\n");
printf("Confirming equality of some derivatives...\n");

# Confirm that phi_zeta[0] == phi_zeta[2]
t := simplify(phi_zeta[0+1] - phi_zeta[2+1]);
printf("phi_zeta[0] - phi_zeta[2] = %a\n", t);

# Confirm that phi_zeta[1] == phi_zeta[3]
t := simplify(phi_zeta[1+1] - phi_zeta[3+1]);
printf("phi_zeta[1] - phi_zeta[3] = %a\n", t);

printf("\n");

# Check limits of derivatives as xi->0, eta->0, zeta->1
# (They all remain bounded.)
for i from 1 to 5 do
  t := limit(subs(xi=0, eta=0, phi_xi[i]), zeta=1);
  printf("limit(phi_xi[%d]) as zeta->1 = %a\n", i-1, t);
end do;

printf("\n");

for i from 1 to 5 do
  t := limit(subs(xi=0, eta=0, phi_eta[i]), zeta=1);
  printf("limit(phi_eta[%d]) as zeta->1 = %a\n", i-1, t);
end do;

printf("\n");

for i from 1 to 5 do
  t := limit(subs(xi=0, eta=0, phi_zeta[i]), zeta=1);
  printf("limit(phi_zeta[%d]) as zeta->1 = %a\n", i-1, t);
end do;

printf("\n");


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

printf("\n");

# Print results
for i from 1 to 5 do
  printf("phi_eta[%d] = %a \n", i-1, phi_eta[i]);
end do;

printf("\n");

# Print results
for i from 1 to 5 do
  printf("phi_zeta[%d] = %a \n", i-1, phi_zeta[i]);
end do;

# Test that first derivaties sum to zero.
printf("\n");
printf("Checking that first derivatives sum to zero...");

phi_xi_sum := 0;
phi_eta_sum := 0;
phi_zeta_sum := 0;

for i from 1 to 5 do
  phi_xi_sum := phi_xi_sum + phi_xi[i];
  phi_eta_sum := phi_eta_sum + phi_eta[i];
  phi_zeta_sum := phi_zeta_sum + phi_zeta[i];
end do;

phi_xi_sum := subs(eps=0, simplify(phi_xi_sum));
phi_eta_sum := subs(eps=0, simplify(phi_eta_sum));
phi_zeta_sum := subs(eps=0, simplify(phi_zeta_sum));

if (phi_xi_sum <> 0) then
  printf("phi_xi_sum = %a\n", phi_xi_sum);
  error;
end if;

if (phi_eta_sum <> 0) then
  printf("phi_eta_sum = %a\n", phi_eta_sum);
  error;
end if;

if (phi_zeta_sum <> 0) then
  printf("phi_zeta_sum = %a\n", phi_zeta_sum);
  error;
end if;

printf("success!\n");


# Compute second derivatives
phi_xixi := Vector(5);
phi_etaeta := Vector(5);
phi_zetazeta := Vector(5);

phi_xieta := Vector(5);
phi_xizeta := Vector(5);
phi_etazeta := Vector(5);

for i from 1 to 5 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:

# Do some simplifications
for i from 1 to 5 do
  phi_zetazeta[i] := simplify(phi_zetazeta[i]);
  phi_xizeta[i] := simplify(phi_xizeta[i]);
  phi_etazeta[i] := simplify(phi_etazeta[i]);

  # Let eps=0 in the numerator of each term
  phi_zetazeta[i] := subs(eps=0, numer(phi_zetazeta[i])) / denom(phi_zetazeta[i]);
  phi_xizeta[i] := subs(eps=0, numer(phi_xizeta[i])) / denom(phi_xizeta[i]);
  phi_etazeta[i] := subs(eps=0, numer(phi_etazeta[i])) / denom(phi_etazeta[i]);
end do:


# Print second derivatives

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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

# Compute equations of planes defining the triangular sides of the pyramid
x := Vector(5);

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]);
x[4 + 1]  := Vector(3, [ 0,   0,  1]);

# None of this did anything...
#geom3d[plane](p0, [x[0+1], x[1+1], x[4+1]]);
#detail(p0);
#Equation(p0, [xi,eta,zeta]);
#printf("p0 = %a\n", p0);

# Normalize uses the infinity-norm by default... actually we don't care
# what the norm is as long as its consistent, so pick something that gives
# you 'nice' numbers in the normal vectors?
n := Vector(5);
n[1] := Normalize(CrossProduct(x[1+1] - x[0+1], x[4+1] - x[0+1]));
n[2] := Normalize(CrossProduct(x[2+1] - x[1+1], x[4+1] - x[1+1]));
n[3] := Normalize(CrossProduct(x[3+1] - x[2+1], x[4+1] - x[2+1]));
n[4] := Normalize(CrossProduct(x[0+1] - x[3+1], x[4+1] - x[3+1]));
n[5] := Vector(3, [0, 0, -1]);

printf("\n");

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

printf("\n");

# Compute equation of the plane for each side.
xi_vec := Vector(3, [xi, eta, zeta]);
planes := Vector(Dimension(n));
for i from 1 to 4 do
  planes[i] := DotProduct(n[i], (xi_vec - x[i]));
end do;

# The bottom plane uses the point (0,0,0) which is on the bottom of
# the reference pyramid, for convenience.
planes[5] := DotProduct(n[i], (xi_vec - Vector(3, [0,0,0])));

# Print equations of planes
for i from 1 to Dimension(planes) do
  printf("plane[%d] = %a\n", i-1, planes[i]);
end do;


return;

end proc;

# Results

# phi_xi[0] = 1/4*(zeta+eta-1)/(1-zeta)
# phi_xi[1] = -1/4*(zeta+eta-1)/(1-zeta)
# phi_xi[2] = -1/4*(zeta-eta-1)/(1-zeta)
# phi_xi[3] = 1/4*(zeta-eta-1)/(1-zeta)
# phi_xi[4] = 0
#
# phi_eta[0] = 1/4*(zeta+xi-1)/(1-zeta)
# phi_eta[1] = 1/4*(zeta-xi-1)/(1-zeta)
# phi_eta[2] = -1/4*(zeta-xi-1)/(1-zeta)
# phi_eta[3] = -1/4*(zeta+xi-1)/(1-zeta)
# phi_eta[4] = 0
#
# With eps identically = 0:
# phi_zeta[0] = 1/4*(-zeta^2+2*zeta-1+xi*eta)/(zeta-1)^2
# phi_zeta[1] = -1/4*(zeta^2-2*zeta+1+xi*eta)/(zeta-1)^2
# phi_zeta[2] = 1/4*(-zeta^2+2*zeta-1+xi*eta)/(zeta-1)^2
# phi_zeta[3] = -1/4*(zeta^2-2*zeta+1+xi*eta)/(zeta-1)^2
# phi_zeta[4] = 1
#
# With general eps:
# phi_zeta[0] = 1/4*(2*zeta-zeta^2+2*zeta*eps+eta*eps-1-2*eps+xi*eps+xi*eta)/(1-zeta+eps)^2
# phi_zeta[1] = 1/4*(2*zeta-zeta^2+2*zeta*eps+eta*eps-1-2*eps-xi*eps-xi*eta)/(1-zeta+eps)^2
# phi_zeta[2] = -1/4*(-2*zeta+zeta^2-2*zeta*eps+eta*eps+1+2*eps+xi*eps-xi*eta)/(1-zeta+eps)^2
# phi_zeta[3] = -1/4*(-2*zeta+zeta^2-2*zeta*eps+eta*eps+1+2*eps-xi*eps+xi*eta)/(1-zeta+eps)^2
# phi_zeta[4] = 1

# And letting eps tend to zero in the numerator:
# phi_zeta[0] = 1/4*(-zeta^2+2*zeta-1+xi*eta)/(1-zeta+eps)^2
# phi_zeta[1] = 1/4*(2*zeta-zeta^2-1-xi*eta)/(1-zeta+eps)^2
# phi_zeta[2] = 1/4*(-zeta^2+2*zeta-1+xi*eta)/(1-zeta+eps)^2
# phi_zeta[3] = 1/4*(2*zeta-zeta^2-1-xi*eta)/(1-zeta+eps)^2
# phi_zeta[4] = 1


# Results - second derivatives.  In all cases, we compute the
# derivatives with general eps in place and then let eps->0 in the
# numerator only.

# phi_xixi[0] = 0
# phi_xixi[1] = 0
# phi_xixi[2] = 0
# phi_xixi[3] = 0
# phi_xixi[4] = 0
#
# phi_etaeta[0] = 0
# phi_etaeta[1] = 0
# phi_etaeta[2] = 0
# phi_etaeta[3] = 0
# phi_etaeta[4] = 0
#
# phi_zetazeta[0] = 1/2*xi*eta/(1-zeta+eps)^3
# phi_zetazeta[1] = -1/2*xi*eta/(1-zeta+eps)^3
# phi_zetazeta[2] = 1/2*xi*eta/(1-zeta+eps)^3
# phi_zetazeta[3] = -1/2*xi*eta/(1-zeta+eps)^3
# phi_zetazeta[4] = 0
#
# phi_xieta[0] = 1/4/(1-zeta+eps)
# phi_xieta[1] = -1/4/(1-zeta+eps)
# phi_xieta[2] = 1/4/(1-zeta+eps)
# phi_xieta[3] = -1/4/(1-zeta+eps)
# phi_xieta[4] = 0
#
# phi_xizeta[0] = 1/4*eta/(1-zeta+eps)^2
# phi_xizeta[1] = -1/4*eta/(1-zeta+eps)^2
# phi_xizeta[2] = 1/4*eta/(1-zeta+eps)^2
# phi_xizeta[3] = -1/4*eta/(1-zeta+eps)^2
# phi_xizeta[4] = 0
#
# phi_etazeta[0] = 1/4*xi/(1-zeta+eps)^2
# phi_etazeta[1] = -1/4*xi/(1-zeta+eps)^2
# phi_etazeta[2] = 1/4*xi/(1-zeta+eps)^2
# phi_etazeta[3] = -1/4*xi/(1-zeta+eps)^2
# phi_etazeta[4] = 0
#
# n[0] = Vector(3, [0,-1,1])
# n[1] = Vector(3, [1,0,1])
# n[2] = Vector(3, [0,1,1])
# n[3] = Vector(3, [-1,0,1])
# n[4] = Vector(3, [0,0,-1])
#
# plane[0] = zeta-eta-1
# plane[1] = zeta+xi-1
# plane[2] = zeta+eta-1
# plane[3] = zeta-xi-1
# plane[4] = -zeta