Given PRISM15 shape functions, computes and verifies some properties of their first and second derivatives.

prism15.maple

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

prism15 := proc()

uses LinearAlgebra;

local
 x, phi, i, j, nodal_val, xi_degree, eta_degree, zeta_degree, phi_xi,
 phi_eta, phi_zeta, indexes, phi_xixi, phi_etaeta, phi_zetazeta,
 phi_xieta, phi_xizeta, phi_etazeta, phi_sum, sum_phi_xi, sum_phi_eta,
 sum_phi_zeta, d2phi_sums;

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

# The nodal points of the PRISM15 reference element:
x := Vector(15);

x[0 + 1]  := Vector(3, [0,   0,   -1]);
x[1 + 1]  := Vector(3, [1,   0,   -1]);
x[2 + 1]  := Vector(3, [0,   1,   -1]);
x[3 + 1]  := Vector(3, [0,   0,    1]);
x[4 + 1]  := Vector(3, [1,   0,    1]);
x[5 + 1]  := Vector(3, [0,   1,    1]);
x[6 + 1]  := Vector(3, [1/2, 0,   -1]);
x[7 + 1]  := Vector(3, [1/2, 1/2, -1]);
x[8 + 1]  := Vector(3, [0,   1/2, -1]);
x[9 + 1]  := Vector(3, [0,   0,   0]);
x[10 + 1] := Vector(3, [1,   0,   0]);
x[11 + 1] := Vector(3, [0,   1,   0]);
x[12 + 1] := Vector(3, [1/2, 0,   1]);
x[13 + 1] := Vector(3, [1/2, 1/2, 1]);
x[14 + 1] := Vector(3, [0,   1/2, 1]);

# The (proposed) shape functions

phi := Vector(15);

phi[0 + 1]  := (1 - zeta) * (xi + eta - 1) * (xi + eta + zeta/2);
phi[1 + 1]  := (1-zeta) * xi  * (xi  - 1 - zeta/2);
phi[2 + 1]  := (1-zeta) * eta * (eta - 1 - zeta/2); # phi1 with xi <- eta
phi[3 + 1]  := (1+zeta) * (xi + eta - 1) * (xi + eta - zeta/2); # phi0 with zeta <- (-zeta)
phi[4 + 1]  := (1+zeta) * xi  * (xi  - 1 + zeta/2); # phi1 with zeta <- (-zeta)
phi[5 + 1]  := (1+zeta) * eta * (eta - 1 + zeta/2); # phi4 with xi <- eta
phi[6 + 1]  := 2 * (1 - zeta) * xi * (1 - xi - eta);
phi[7 + 1]  := 2 * (1 - zeta) * xi * eta;
phi[8 + 1]  := 2 * (1 - zeta) * eta * (1 - xi - eta);
phi[9 + 1]  := (1 - zeta) * (1 + zeta) * (1 - xi - eta);
phi[10 + 1] := (1 - zeta) * (1 + zeta) * xi;
phi[11 + 1] := (1 - zeta) * (1 + zeta) * eta;
phi[12 + 1] := 2 * (1 + zeta) * xi * (1 - xi - eta); # phi6 with zeta <- (-zeta)
phi[13 + 1] := 2 * (1 + zeta) * xi * eta; # phi7 with zeta <- (-zeta)
phi[14 + 1] := 2 * (1 + zeta) * eta * (1 - xi - eta); # phi8 with zeta <- (-zeta)

# Check that phi_i(x_j) = delta_{ij} is satisfied.
for i from 1 to 15 do
  printf("Checking interpolary property of phi[%d]...", i-1);
  for j from 1 to 15 do
    nodal_val := subs(xi=x[j][1], eta=x[j][2], zeta=x[j][3], phi[i]);
    # printf("....................nodal_val = %a \n", nodal_val);

    if (i = j) then
        if (nodal_val <> 1) then
            printf("Error: phi[%d] does not satisfy phi_i(x_i)==1.\n", i-1);
            error;
        end if;
    end if;

    if (i <> j) then
        if (nodal_val <> 0) then
            printf("Error: phi[%d] does not satisfy phi_i(x_j)==0.\n", i-1);
            error;
        end if;
    end if;
  end do:
  printf("success!\n");
end do:

# For an arbitrary point within the reference element, the shape functions
# should sum to 1
# test_point := Vector(3, [1/3, 1/3, 0]);
# test_point := Vector(3, [1/2, 1/3, 1/2]);
# test_point := Vector(3, [1/6, 5/6, 1]);
# test_value := 0;
# for i from 1 to 15 do
#   test_value := test_value +
#     subs(xi=test_point[1], eta=test_point[2], zeta=test_point[3], phi[i]);
# end do:
# printf("test_value = %a\n", test_value);

# Maple should be able to symbolically sum all the phi values to 1 regardless of
# xi, eta, and zeta?
printf("Checking that phi functions sum to 1 symbolically...");
phi_sum := 0;
for i from 1 to 15 do
  phi_sum := phi_sum + phi[i];
end do:

phi_sum := simplify(phi_sum);

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

# Check that the maximum polynomial order of all the phi's is 2?
for i from 1 to 15 do
  printf("Checking polynomial degree of phi[%d]...", i-1);

  # Polynomial must be in collected form for degree to work
  xi_degree := degree(collect(phi[i], xi), xi);
  # printf("....................xi_degree = %a \n", xi_degree);
  if (xi_degree > 2) then
      printf("xi_degree for phi[%d] is > 2!\n", i-1);
      error;
  end if;

  eta_degree := degree(collect(phi[i], eta), eta);
  # printf("....................eta_degree = %a \n", eta_degree);
  if (eta_degree > 2) then
      printf("eta_degree for phi[%d] is > 2!\n", i-1);
      error;
  end if;

  zeta_degree := degree(collect(phi[i], zeta), zeta);
  # printf("....................zeta_degree = %a \n", zeta_degree);
  if (zeta_degree > 2) then
      printf("zeta_degree for phi[%d] is > 2!\n", i-1);
      error;
  end if;

  printf("success!\n");
end do;



# Compute first derivatives

phi_xi := Vector(15);
phi_eta := Vector(15);
phi_zeta := Vector(15);

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

# Collect in (1-zeta) for some phi_xi derivatives
indexes := Vector(3, [1, 2, 7]);
for i from 1 to Dimension(indexes) do
  phi_xi[ indexes[i] ] := subs(tmp=1-zeta, collect(subs(1-zeta=tmp, phi_xi[ indexes[i] ]), tmp));
end do:
# Collect in (1+zeta) for some phi_xi derivatives
indexes := Vector(3, [4, 5, 13]);
for i from 1 to Dimension(indexes) do
  phi_xi[ indexes[i] ] := subs(tmp=1+zeta, collect(subs(1+zeta=tmp, phi_xi[ indexes[i] ]), tmp));
end do:


# Collect in (1-zeta) for some phi_eta derivatives
indexes := Vector(3, [1, 3, 9]);
for i from 1 to Dimension(indexes) do
  phi_eta[ indexes[i] ] := subs(tmp=1-zeta, collect(subs(1-zeta=tmp, phi_eta[ indexes[i] ]), tmp));
end do:
# Collect in (1+zeta) for some phi_eta derivatives
indexes := Vector(3, [4, 6, 15]);
for i from 1 to Dimension(indexes) do
  phi_eta[ indexes[i] ] := subs(tmp=1+zeta, collect(subs(1+zeta=tmp, phi_eta[ indexes[i] ]), tmp));
end do:


# Collect in (xi+eta-1) for some phi_zeta derivatives
indexes := Vector(2, [1, 4]);
for i from 1 to Dimension(indexes) do
  phi_zeta[ indexes[i] ] := subs(tmp=xi+eta-1, collect(subs(xi+eta-1=tmp, phi_zeta[ indexes[i] ]), tmp));
end do:

phi_zeta[ 10 ] := subs(tmp=-xi-eta+1, collect(subs(-xi-eta+1=tmp, phi_zeta[ 10 ]), tmp));

# Collect in xi for some phi_zeta derivatives
indexes := Vector(2, [2, 5]);
for i from 1 to Dimension(indexes) do
  phi_zeta[ indexes[i] ] := factor(collect(phi_zeta[ indexes[i] ], xi));
end do:

# Collect in eta for some phi_zeta derivatives
indexes := Vector(2, [3, 6]);
for i from 1 to Dimension(indexes) do
  phi_zeta[ indexes[i] ] := factor(collect(phi_zeta[ indexes[i] ], eta));
end do:

# Simplify some phi_zeta derivatives
indexes := Vector(2, [11, 12]);
for i from 1 to Dimension(indexes) do
  phi_zeta[ indexes[i] ] := simplify(phi_zeta[ indexes[i] ]);
end do:


# Summing the phi_xi, phi_eta, or phi_zeta should result in zero,
# since sum(phi)==1.
printf("Checking that dphi functions sum to 0 symbolically...");
sum_phi_xi := 0;
sum_phi_eta := 0;
sum_phi_zeta := 0;
for i from 1 to 15 do
  sum_phi_xi := sum_phi_xi + phi_xi[i];
  sum_phi_eta := sum_phi_eta + phi_eta[i];
  sum_phi_zeta := sum_phi_zeta + phi_zeta[i];
end do:

sum_phi_xi := simplify(sum_phi_xi);
sum_phi_eta := simplify(sum_phi_eta);
sum_phi_zeta := simplify(sum_phi_zeta);

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

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

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

printf("success!\n");


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

printf("\n");

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

printf("\n");

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


# Compute second derivatives

phi_xixi := Vector(15);
phi_etaeta := Vector(15);
phi_zetazeta := Vector(15);

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

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


# As far as I can tell, it's not possible to get Maple to write
# 2-2*zeta as 2*(1-zeta) using simplify, collect, expand, combine,
# subs, algsubs, subsop, convert, freeze, thaw, or any other function...
# for i from 1 to 15 do
#   phi_xizeta[ i ] := factor(phi_xizeta[ i ]);
# end do:

# The first derivatives must sum to zero, so the second derivatives must as well.
printf("Checking that second derivatives sum to zero...");
d2phi_sums := Vector(6);
for i from 1 to 15 do
  d2phi_sums[1] := d2phi_sums[1] + phi_xixi[i];
  d2phi_sums[2] := d2phi_sums[2] + phi_etaeta[i];
  d2phi_sums[3] := d2phi_sums[3] + phi_zetazeta[i];
  d2phi_sums[4] := d2phi_sums[4] + phi_xieta[i];
  d2phi_sums[5] := d2phi_sums[5] + phi_xizeta[i];
  d2phi_sums[6] := d2phi_sums[6] + phi_etazeta[i];
end do:

for i from 1 to Dimension(d2phi_sums) do
  if (d2phi_sums[i] <> 0) then
      printf("Error, dphi_sums[%d] did not sum to zero!", i);
      error;
  end if;
end do:
printf("success!\n");



# Print second derivatives

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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

printf("\n");

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


return;

end proc;

# Results - first derivatives

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


# Results - second derivatives

# phi_xixi[0] = 2-2*zeta
# phi_xixi[1] = 2-2*zeta
# phi_xixi[2] = 0
# phi_xixi[3] = 2+2*zeta
# phi_xixi[4] = 2+2*zeta
# phi_xixi[5] = 0
# phi_xixi[6] = 4*zeta-4
# phi_xixi[7] = 0
# phi_xixi[8] = 0
# phi_xixi[9] = 0
# phi_xixi[10] = 0
# phi_xixi[11] = 0
# phi_xixi[12] = -4-4*zeta
# phi_xixi[13] = 0
# phi_xixi[14] = 0
#
# phi_etaeta[0] = 2-2*zeta
# phi_etaeta[1] = 0
# phi_etaeta[2] = 2-2*zeta
# phi_etaeta[3] = 2+2*zeta
# phi_etaeta[4] = 0
# phi_etaeta[5] = 2+2*zeta
# phi_etaeta[6] = 0
# phi_etaeta[7] = 0
# phi_etaeta[8] = 4*zeta-4
# phi_etaeta[9] = 0
# phi_etaeta[10] = 0
# phi_etaeta[11] = 0
# phi_etaeta[12] = 0
# phi_etaeta[13] = 0
# phi_etaeta[14] = -4-4*zeta
#
# phi_zetazeta[0] = -xi-eta+1
# phi_zetazeta[1] = xi
# phi_zetazeta[2] = eta
# phi_zetazeta[3] = -xi-eta+1
# phi_zetazeta[4] = xi
# phi_zetazeta[5] = eta
# phi_zetazeta[6] = 0
# phi_zetazeta[7] = 0
# phi_zetazeta[8] = 0
# phi_zetazeta[9] = 2*xi+2*eta-2
# phi_zetazeta[10] = -2*xi
# phi_zetazeta[11] = -2*eta
# phi_zetazeta[12] = 0
# phi_zetazeta[13] = 0
# phi_zetazeta[14] = 0
#
# phi_xieta[0] = 2-2*zeta
# phi_xieta[1] = 0
# phi_xieta[2] = 0
# phi_xieta[3] = 2+2*zeta
# phi_xieta[4] = 0
# phi_xieta[5] = 0
# phi_xieta[6] = 2*zeta-2
# phi_xieta[7] = 2-2*zeta
# phi_xieta[8] = 2*zeta-2
# phi_xieta[9] = 0
# phi_xieta[10] = 0
# phi_xieta[11] = 0
# phi_xieta[12] = -2-2*zeta
# phi_xieta[13] = 2+2*zeta
# phi_xieta[14] = -2-2*zeta
#
# phi_xizeta[0] = 3/2-zeta-2*xi-2*eta
# phi_xizeta[1] = 1/2+zeta-2*xi
# phi_xizeta[2] = 0
# phi_xizeta[3] = -3/2-zeta+2*xi+2*eta
# phi_xizeta[4] = -1/2+zeta+2*xi
# phi_xizeta[5] = 0
# phi_xizeta[6] = 4*xi+2*eta-2
# phi_xizeta[7] = -2*eta
# phi_xizeta[8] = 2*eta
# phi_xizeta[9] = 2*zeta
# phi_xizeta[10] = -2*zeta
# phi_xizeta[11] = 0
# phi_xizeta[12] = -4*xi-2*eta+2
# phi_xizeta[13] = 2*eta
# phi_xizeta[14] = -2*eta
#
# phi_etazeta[0] = 3/2-zeta-2*xi-2*eta
# phi_etazeta[1] = 0
# phi_etazeta[2] = 1/2+zeta-2*eta
# phi_etazeta[3] = -3/2-zeta+2*xi+2*eta
# phi_etazeta[4] = 0
# phi_etazeta[5] = -1/2+zeta+2*eta
# phi_etazeta[6] = 2*xi
# phi_etazeta[7] = -2*xi
# phi_etazeta[8] = 2*xi+4*eta-2
# phi_etazeta[9] = 2*zeta
# phi_etazeta[10] = 0
# phi_etazeta[11] = -2*zeta
# phi_etazeta[12] = -2*xi
# phi_etazeta[13] = 2*xi
# phi_etazeta[14] = -2*xi-4*eta+2