hex20.maple

# Compute second derivatives of Hex20 basis functions using Maple

hex20 := proc()

uses LinearAlgebra;

local N, i, phi, phi_x, phi_xx, phi_y, phi_yy, phi_z, phi_zz, phi_xy, phi_xz, phi_yz, indexes;

N:=20;

phi := Vector(N);
phi_x := Vector(N);
phi_xx := Vector(N);

phi_y := Vector(N);
phi_yy := Vector(N);

phi_z := Vector(N);
phi_zz := Vector(N);

phi_xy := Vector(N);
phi_xz := Vector(N);
phi_yz := Vector(N);

# Note: 1-based numbering... these functions are copied directy from
# libmesh code.
phi[0 + 1]  := (1 - x)*(1 - y)*(1 - z)*(1 - 2*x - 2*y - 2*z);
phi[1 + 1]  := x*(1 - y)*(1 - z)*(2*x - 2*y - 2*z - 1);
phi[2 + 1]  := x*y*(1 - z)*(2*x + 2*y - 2*z - 3);
phi[3 + 1]  := (1 - x)*y*(1 - z)*(2*y - 2*x - 2*z - 1);
phi[4 + 1]  := (1 - x)*(1 - y)*z*(2*z - 2*x - 2*y - 1);
phi[5 + 1]  := x*(1 - y)*z*(2*x - 2*y + 2*z - 3);
phi[6 + 1]  := x*y*z*(2*x + 2*y + 2*z - 5);
phi[7 + 1]  := (1 - x)*y*z*(2*y - 2*x + 2*z - 3);
phi[8 + 1]  := 4*x*(1 - x)*(1 - y)*(1 - z);
phi[9 + 1]  := 4*x*y*(1 - y)*(1 - z);
phi[10 + 1] := 4*x*(1 - x)*y*(1 - z);
phi[11 + 1] := 4*(1 - x)*y*(1 - y)*(1 - z);
phi[12 + 1] := 4*(1 - x)*(1 - y)*z*(1 - z);
phi[13 + 1] := 4*x*(1 - y)*z*(1 - z);
phi[14 + 1] := 4*x*y*z*(1 - z);
phi[15 + 1] := 4*(1 - x)*y*z*(1 - z);
phi[16 + 1] := 4*x*(1 - x)*(1 - y)*z;
phi[17 + 1] := 4*x*y*(1 - y)*z;
phi[18 + 1] := 4*x*(1 - x)*y*z;
phi[19 + 1] := 4*(1 - x)*y*(1 - y)*z;

# # Print basis functions 0-based
# for i from 1 to N do
#   printf("....................phi[%d] = \n", i-1);
#   print( phi[i] );
# end do:

# Compute first derivatives
for i from 1 to N do
  # Note: extra 1/2 comes from xi-x transformation of variables
  phi_x[i] := (1/2) * diff(phi[i], x);
  phi_y[i] := (1/2) * diff(phi[i], y);
  phi_z[i] := (1/2) * diff(phi[i], z);

  # Print derivatives 0-based
  # printf("....................phi_x[%d] = \n", i-1);
  # print( phi_x[i] );
end do:


# Compute second derivatives.  Another 1/2 comes from differentiating
# again.
for i from 1 to N do
  phi_xx[i] := (1/2) * diff(phi_x[i], x);
  phi_yy[i] := (1/2) * diff(phi_y[i], y);
  phi_zz[i] := (1/2) * diff(phi_z[i], z);

  phi_xy[i] := (1/2) * diff(phi_x[i], y);
  phi_xz[i] := (1/2) * diff(phi_x[i], z);
  phi_yz[i] := (1/2) * diff(phi_y[i], z);
end do:

# Do specific simplifications (collect in 1-z by jumping through some
# hoops) on certain terms
indexes := Vector(8, [1, 2, 3, 4, 9, 10, 11, 12]);
for i from 1 to 8 do
  phi_xy[ indexes[i] ] := subs(tmp=1-z, collect(subs(1-z=tmp, phi_xy[ indexes[i] ]), tmp));
end do:

# Do simplifications by collecting in 1-y
indexes := Vector(8, [1, 2, 5, 6, 9, 13, 14, 17]);
for i from 1 to 8 do
  phi_xz[ indexes[i] ] := subs(tmp=1-y, collect(subs(1-y=tmp, phi_xz[ indexes[i] ]), tmp));
end do:

# Do simplifications by collecting in 1-x
indexes := Vector(8, [1, 4, 5, 8, 12, 13, 16, 20]);
for i from 1 to 8 do
  phi_yz[ indexes[i] ] := subs(tmp=1-x, collect(subs(1-x=tmp, phi_yz[ indexes[i] ]), tmp));
end do:

# Make a few other straightforward simplifications and re-collect
indexes := Vector(4, [17, 18, 19, 20]);
for i from 1 to 4 do
  phi_xy[ indexes[i] ] := collect(simplify(phi_xy[ indexes[i] ]), z);
end do:

indexes := Vector(4, [11, 15, 16, 19]);
for i from 1 to 4 do
  phi_xz[ indexes[i] ] := collect(simplify(phi_xz[ indexes[i] ]), y);
end do:

indexes := Vector(4, [10, 14, 15, 18]);
for i from 1 to 4 do
  phi_yz[ indexes[i] ] := collect(simplify(phi_yz[ indexes[i] ]), x);
end do:

# And a few other simplifications based on factoring...
indexes := Vector(4, [5, 6, 7, 8]);
for i from 1 to 4 do
  phi_xy[ indexes[i] ] := factor(simplify(phi_xy[ indexes[i] ]));
end do:

indexes := Vector(4, [3, 4, 7, 8]);
for i from 1 to 4 do
  phi_xz[ indexes[i] ] := factor(simplify(phi_xz[ indexes[i] ]));
end do:

indexes := Vector(4, [2, 3, 6, 7]);
for i from 1 to 4 do
  phi_yz[ indexes[i] ] := factor(simplify(phi_yz[ indexes[i] ]));
end do:



# Print all xx derivatives
for i from 1 to N do
  printf("....................phi_xx[%d] = %a \n", i-1, phi_xx[i]);
end do:

printf("\n");

# Print all yy derivatives
for i from 1 to N do
  printf("....................phi_yy[%d] = %a \n", i-1, phi_yy[i]);
end do:

printf("\n");

# Print all zz derivatives
for i from 1 to N do
  printf("....................phi_zz[%d] = %a \n", i-1, phi_zz[i]);
end do:

printf("\n");

# Print all xy derivatives
for i from 1 to N do
  printf("....................phi_xy[%d] = %a \n", i-1, phi_xy[i]);
end do:

printf("\n");

# Print all xz derivatives
for i from 1 to N do
  printf("....................phi_xz[%d] = %a \n", i-1, phi_xz[i]);
end do:

printf("\n");

# Print all yz derivatives
for i from 1 to N do
  printf("....................phi_yz[%d] = %a \n", i-1, phi_yz[i]);
end do:

return;

end proc;

################################################################################
# Results: xx derivative, re-ordered
# phi_xx[0] = (1-y)*(1-z)
# phi_xx[1] = (1-y)*(1-z)

# phi_xx[2] = y*(1-z)
# phi_xx[3] = y*(1-z)

# phi_xx[4] = (1-y)*z
# phi_xx[5] = (1-y)*z

# phi_xx[6] = y*z
# phi_xx[7] = y*z

# phi_xx[8] = -2*(1-y)*(1-z)

# phi_xx[10] = -2*y*(1-z)

# phi_xx[16] = -2*(1-y)*z

# phi_xx[18] = -2*y*z

# phi_xx[9] = 0
# phi_xx[11] = 0
# phi_xx[12] = 0
# phi_xx[13] = 0
# phi_xx[14] = 0
# phi_xx[15] = 0
# phi_xx[17] = 0
# phi_xx[19] = 0


################################################################################
# Results: yy derivative, re-ordered
# phi_yy[0] = (1-x)*(1-z)
# phi_yy[3] = (1-x)*(1-z)

# phi_yy[1] = x*(1-z)
# phi_yy[2] = x*(1-z)

# phi_yy[4] = (1-x)*z
# phi_yy[7] = (1-x)*z

# phi_yy[5] = x*z
# phi_yy[6] = x*z

# phi_yy[9] = -2*x*(1-z)

# phi_yy[11] = -2*(1-x)*(1-z)

# phi_yy[17] = -2*x*z

# phi_yy[19] = -2*(1-x)*z

# phi_yy[8] = 0
# phi_yy[10] = 0
# phi_yy[12] = 0
# phi_yy[13] = 0
# phi_yy[14] = 0
# phi_yy[15] = 0
# phi_yy[16] = 0
# phi_yy[18] = 0


################################################################################
# Results: zz derivative, re-ordered
# phi_zz[0] = (1-x)*(1-y)
# phi_zz[4] = (1-x)*(1-y)

# phi_zz[1] = x*(1-y)
# phi_zz[5] = x*(1-y)

# phi_zz[2] = x*y
# phi_zz[6] = x*y

# phi_zz[3] = (1-x)*y
# phi_zz[7] = (1-x)*y

# phi_zz[12] = -2*(1-x)*(1-y)

# phi_zz[13] = -2*x*(1-y)

# phi_zz[14] = -2*x*y

# phi_zz[15] = -2*(1-x)*y

# phi_zz[8] = 0
# phi_zz[9] = 0
# phi_zz[10] = 0
# phi_zz[11] = 0
# phi_zz[16] = 0
# phi_zz[17] = 0
# phi_zz[18] = 0
# phi_zz[19] = 0


################################################################################
# Results: xy derivatives
# phi_xy[0] = (5/4-x-y-1/2*z)*(1-z)
# phi_xy[1] = (-x+y+1/2*z-1/4)*(1-z)
# phi_xy[2] = (x+y-1/2*z-3/4)*(1-z)
# phi_xy[3] = (-y+x+1/2*z-1/4)*(1-z)
# phi_xy[4] = -1/4*z*(4*x+4*y-2*z-3)
# phi_xy[5] = -1/4*z*(-4*y+4*x+2*z-1)
# phi_xy[6] = 1/4*z*(-5+4*x+4*y+2*z)
# phi_xy[7] = 1/4*z*(4*x-4*y-2*z+1)
# phi_xy[8] = (-1+2*x)*(1-z)
# phi_xy[9] = (1-2*y)*(1-z)
# phi_xy[10] = (1-2*x)*(1-z)
# phi_xy[11] = (-1+2*y)*(1-z)
# phi_xy[12] = z*(1-z)
# phi_xy[13] = -z*(1-z)
# phi_xy[14] = z*(1-z)
# phi_xy[15] = -z*(1-z)
# phi_xy[16] = (-1+2*x)*z
# phi_xy[17] = (1-2*y)*z
# phi_xy[18] = (1-2*x)*z
# phi_xy[19] = (-1+2*y)*z


################################################################################
# Results: xz derivatives
# phi_xz[0] = (5/4-x-1/2*y-z)*(1-y)
# phi_xz[1] = (-x+1/2*y+z-1/4)*(1-y)
# phi_xz[2] = -1/4*y*(2*y+4*x-4*z-1)
# phi_xz[3] = -1/4*y*(-2*y+4*x+4*z-3)
# phi_xz[4] = (-z+x+1/2*y-1/4)*(1-y)
# phi_xz[5] = (x-1/2*y+z-3/4)*(1-y)
# phi_xz[6] = 1/4*y*(2*y+4*x+4*z-5)
# phi_xz[7] = 1/4*y*(-2*y+4*x-4*z+1)
# phi_xz[8] = (-1+2*x)*(1-y)
# phi_xz[9] = -y*(1-y)
# phi_xz[10] = (-1+2*x)*y
# phi_xz[11] = y*(1-y)
# phi_xz[12] = (-1+2*z)*(1-y)
# phi_xz[13] = (1-2*z)*(1-y)
# phi_xz[14] = (1-2*z)*y
# phi_xz[15] = (-1+2*z)*y
# phi_xz[16] = (1-2*x)*(1-y)
# phi_xz[17] = y*(1-y)
# phi_xz[18] = (1-2*x)*y
# phi_xz[19] = -y*(1-y)


################################################################################
# Results: yz derivatives
# phi_yz[0] = (5/4-1/2*x-y-z)*(1-x)
# phi_yz[1] = 1/4*x*(2*x-4*y-4*z+3)
# phi_yz[2] = -1/4*x*(2*x+4*y-4*z-1)
# phi_yz[3] = (-y+1/2*x+z-1/4)*(1-x)
# phi_yz[4] = (-z+1/2*x+y-1/4)*(1-x)
# phi_yz[5] = -1/4*x*(2*x-4*y+4*z-1)
# phi_yz[6] = 1/4*x*(2*x+4*y+4*z-5)
# phi_yz[7] = (y-1/2*x+z-3/4)*(1-x)
# phi_yz[8] = x*(1-x)
# phi_yz[9] = (-1+2*y)*x
# phi_yz[10] = -x*(1-x)
# phi_yz[11] = (-1+2*y)*(1-x)
# phi_yz[12] = (-1+2*z)*(1-x)
# phi_yz[13] = (-1+2*z)*x
# phi_yz[14] = (1-2*z)*x
# phi_yz[15] = (1-2*z)*(1-x)
# phi_yz[16] = -x*(1-x)
# phi_yz[17] = (1-2*y)*x
# phi_yz[18] = x*(1-x)
# phi_yz[19] = (1-2*y)*(1-x)