Question: How do I optimize algebraic computations of derivatives, matrices inverse and determinants?

Dear users,

I have a Maple file to evaluate the shape functions on hexahedrons and to save the results in a .cpp file. However, the execution time is very long and when I do not loose the connection with the server, I ran out of memory (I have 32GB!!). So, I very keen to hear some tips on how to optimize it! I am sure this computation is simple enough for the Maple powerful and I quite sure there are some tricks to optimize it! The code is below.

Thanks in advance!

# =====================
# reference 2nd order hexadron (VTK numbering)
# =====================

pc := Matrix(20,3,[
[-1, -1, -1], # 1st coordinate.
[1, -1, -1],
[1, 1, -1],
[-1, 1, -1],
[-1, -1, 1], # 5th coordinate.
[1, -1, 1],
[1, 1, 1],
[-1, 1, 1],  # 8th coordinate
[0, -1, -1], # 9th coordinate
[1, 0, -1],
[0, 1, -1],
[-1, 0, -1], # 12th coordinate
[0, -1, 1], # 13th coordinate
[1, 0, 1],
[0, 1, 1],
[-1, 0, 1], # 16th coordinate
[-1, -1, 0], # 17th coordinate.
[1, -1, 0],
[1, 1, 0],
[-1, 1, 0]  # 20th coordinate
] );

# ==============================
# definition of the 2nd order shape functions
# ==============================

for k from 1 to 8 do
N_2[k] := 1/8 * (1+x*pc[k,1]) * (1+ y*pc[k,2]) * (1 + z*pc[k,3]) * (x*pc[k,1] + y*pc[k,2] + z*pc[k,3] - 2);
end do:

N_2[9] := 1/4 * (1-x*x) * (1+ y*pc[9,2]) * (1 + z*pc[9,3]):
N_2[10] := 1/4 * (1-y*y) * (1+ x*pc[10,1]) * (1 + z*pc[10,3]):
N_2[11] := 1/4 * (1-x*x) * (1+ y*pc[11,2]) * (1 + z*pc[11,3]):
N_2[12] := 1/4 * (1-y*y) * (1+ x*pc[12,1]) * (1 + z*pc[12,3]):
N_2[13] := 1/4 * (1-x*x) * (1+ y*pc[13,2]) * (1 + z*pc[13,3]):
N_2[14] := 1/4 * (1-y*y) * (1+ x*pc[14,1]) * (1 + z*pc[14,3]):
N_2[15] := 1/4 * (1-x*x) * (1+ y*pc[15,2]) * (1 + z*pc[15,3]):
N_2[16] := 1/4 * (1-y*y) * (1+ x*pc[16,1]) * (1 + z*pc[16,3]):

for k from 17 to 20 do
N_2[k] := 1/4 * (1 - z*z) * (1 + x*pc[k,1]) * (1 + y*pc[k,2]);
end do:

# TESTS
#for k from 1 to 20 do
#eval( N_2[k], [x=pc[k,1],y=pc[k,2],z=pc[k,3]] );
#end do:

# ====================================
# Tranformation from reference (2nd order) element
# ====================================

F_2[1] := simplify( add(ax[k]*N_2[k], k = 1..20), size ):
F_2[2] := simplify( add(ay[k]*N_2[k], k = 1..20), size ):
F_2[3] := simplify( add(az[k]*N_2[k], k = 1..20), size ):

# TESTS
#eval(F_2[1], [x=-1,y=1,z=1]);
#eval(F_2[2], [x=-1,y=1,z=1]);
#eval(F_2[3], [x=-1,y=1,z=1]);

# =========================
# Compute jacobian and determinant
# =========================

J_2 := Matrix( [
[ simplify( diff(F_2[1], x), size ) , simplify( diff(F_2[2],x), size ) , simplify( diff(F_2[3],x), size ) ],
[ simplify( diff(F_2[1], y), size ) , simplify( diff(F_2[2],y), size ) , simplify( diff(F_2[3],y), size ) ],
[ simplify( diff(F_2[1], z), size ) , simplify( diff(F_2[2],z), size ) , simplify( diff(F_2[3],z), size ) ]
] ):

M_2 := MatrixInverse(J_2):

# TESTS
#eval( M_2, [x=1,y=1,z=1] );

det_2 = abs( simplify( Determinant(M_2), size ) ):

# ====================================
# First order derivatives of 2nd order shape functions
# ====================================

for k from 1 to 20 do
dN_dx_2[k] := simplify( diff(N_2[k], x), size ):
dN_dy_2[k] := simplify( diff(N_2[k], y), size ):
dN_dz_2[k] := simplify( diff(N_2[k], z), size ):
end do:

for k from 1 to 20 do
dn_dx_2[k] := simplify( M_2[1,1] * dN_dx_2[k] + M_2[1,2] * dN_dy_2[k] + M_2[1,3] * dN_dz_2[k], size );
dn_dy_2[k] := simplify( M_2[2,1] * dN_dx_2[k] + M_2[2,2] * dN_dy_2[k] + M_2[2,3] * dN_dz_2[k], size );
dn_dz_2[k] := simplify( M_2[3,1] * dN_dx_2[k] + M_2[3,2] * dN_dy_2[k] + M_2[3,3] * dN_dz_2[k], size );
end do: