VBA code for transforming a 4th rank tensor represented by a 6x6 matrix

[JonMark13]

I have written a code in VBA (derived from python) that takes a 6x6 matrix input and rotates it by 45 degrees. Unfortunately the code returns nothing and I can't see where I have gone wrong:

Function transform_elastic_tensor_2(elastic_tensor As Variant)

Dim transform(1 To 3, 1 To 3) As Double
Dim C_ijkl(1 To 3, 1 To 3, 1 To 3, 1 To 3) As Double
Dim new_C_ijkl(1 To 3, 1 To 3, 1 To 3, 1 To 3) As Double
Dim indices(1 To 9, 1 To 2) As Integer
Dim index_matrix(1 To 3, 1 To 3) As Integer
Dim new_elastic_tensor(1 To 6, 1 To 6) As Double
Dim i As Integer, j As Integer, k As Integer, l As Integer, m As Integer, n As Integer, ij(1 To 2) As Integer, kl(1 To 2) As Integer

'Rotation Matrix
transform(1, 1) = 0.7071068
transform(1, 2) = -0.7071068
transform(1, 3) = 0
transform(2, 1) = 0.7071068
transform(2, 2) = 0.7071068
transform(2, 3) = 0
transform(3, 1) = 0
transform(3, 2) = 0
transform(3, 3) = 1

'Defining 4th rank tensor
For i = 1 To 3
For j = 1 To 3
For k = 1 To 3
For l = 1 To 3
C_ijkl(i, j, k, l) = 0
Next l
Next k
Next j
Next i

'Relation of 4th rank tensor to 6x6 2D representation
indices(1, 1) = 0: indices(1, 2) = 0
indices(2, 1) = 1: indices(2, 2) = 1
indices(3, 1) = 2: indices(3, 2) = 2
indices(4, 1) = 1: indices(4, 2) = 2
indices(5, 1) = 2: indices(5, 2) = 0
indices(6, 1) = 0: indices(6, 2) = 1
indices(7, 1) = 2: indices(7, 2) = 1
indices(8, 1) = 0: indices(8, 2) = 2
indices(9, 1) = 1: indices(9, 2) = 0

index_matrix(1, 1) = 1: index_matrix(1, 2) = 6: index_matrix(1, 3) = 5
index_matrix(2, 1) = 6: index_matrix(2, 2) = 2: index_matrix(2, 3) = 4
index_matrix(3, 1) = 5: index_matrix(3, 2) = 4: index_matrix(3, 3) = 3

For i = 1 To 3
For j = 1 To 3
For k = 1 To 3
For l = 1 To 3
m = index_matrix(i, j) - 1
n = index_matrix(k, l) - 1
C_ijkl(i, j, k, l) = elastic_tensor(m, n)
Next l
Next k
Next j
Next i

'Transformation of 4th rank tensor
For i = 1 To 3
For j = 1 To 3
For k = 1 To 3
For l = 1 To 3
For m = 1 To 3
For n = 1 To 3
new_C_ijkl(i, j, k, l) = new_C_ijkl(i, j, k, l) + transform(m, i) * transform(n, j) * C_ijkl(m, n, k, l)
Next n
Next m
Next l
Next k
Next j
Next i

'Relation of transformed 4th rank tensor to 6x6 2D representation
For i = 1 To 9
ij(1) = indices(i, 1) + 1
ij(2) = indices(i, 2) + 1

For j = 1 To 9
kl(1) = indices(j, 1) + 1
kl(2) = indices(j, 2) + 1

For k = 1 To 3
For l = 1 To 3
For m = 1 To 6
For n = 1 To 6
new_elastic_tensor(index_matrix(ij(1), ij(2)), index_matrix(kl(1), kl(2))) = new_C_ijkl(ij(1), ij(2), k, l)
Next n
Next m
Next l
Next k
Next j
Next i

transform_elastic_tensor_2 = new_elastic_tensor
End Function

The analogous python code is:

def transform_elastic_tensor(elastic_tensor):

#Rotation Matrix
transform = np.array([[0.7071068, -0.7071068, 0],
[0.7071068, 0.7071068, 0],
[0, 0, 1]])




#Defining 4th rank tensor.
C_ijkl = np.zeros((3, 3, 3, 3))

indices = [(0,0), (1,1), (2,2), (1,2), (2,0), (0,1), (2,1), (0,2), (1,0)]

index_matrix = np.array([[1, 6, 5],
[6, 2, 4],
[5, 4, 3]])
#Relation of 4th rank tensor to 6x6 2D representation.
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
m = index_matrix[i, j] - 1
n = index_matrix[k, l] - 1
C_ijkl[i,j,k,l] = elastic_tensor[m,n]

new_C_ijkl = np.einsum('ip,jq,kr,ls,pqrs->ijkl', transform, transform, transform, transform, C_ijkl)
new_elastic_tensor = np.zeros((6,6))
for m in range(elastic_tensor.shape[0]):
for n in range(elastic_tensor.shape[1]):
ij = indices[m]
kl = indices[n]
new_elastic_tensor[m, n] = new_C_ijkl[ij + kl]

return new_elastic_tensor

#Elastic tensor to be rotated, input tensor values
elastic_tensor = np.array([[200, 120, 120, 0, 0, 0],
[120, 200, 120, 0, 0, 0],
[120, 120, 200, 0, 0, 0],
[0, 0, 0, 90, 0, 0],
[0, 0, 0, 0, 90, 0],
[0, 0, 0, 0, 0, 90]])



new_elastic_tensor = transform_elastic_tensor(elastic_tensor)
print(np.round_(new_elastic_tensor, 2))


The correct answer is:
[[250. 70. 120. 0. 0. 0.]
[ 70. 250. 120. 0. 0. 0.]
[120. 120. 200. 0. 0. 0.]
[ 0. 0. 0. 90. 0. 0.]
[ 0. 0. 0. 0. 90. 0.]
[ 0. 0. 0. 0. 0. 40.]]

Any help would be greatly appreciated.

[beyond Excel]

Hello. You don't need to complicate your life because Excel already has matrix multiplication.

You multiply the rotation matrix by the original tensor and you're done.

[JonMark13]

Unfortunately for this matrix multiplication you would have to first convert the 6x6 matrix into a 3x3x3x3 tensor. Then rotate each individual component of the tensor and convert back to the new 6x6 matrix.