double tensorial product for stiffness tensor from Lamé constants

doubletensorproduct.py

import numpy as np

# 2nd order identity tensor
I = np.eye(3)

# Lamé Constants
lam_0 = 456
mu_0 = 123

print(f"G = {mu_0}")
print(f"lambda + 2G = {lam_0 + 2*mu_0}")

def double_tensor(A, B):
    # Double tensorial product, used to build stiffness tensor from Lamé constants
    # Skipping the bases... This assummes e_1 = [1,0,0], e_2 = [0,1,0], e_3 = [0,0,1]
    return 0.5 * (np.einsum('ik,jl', A, B) + np.einsum('il,jk', A, B)) 

# Using the double tensorial product to generate a real 4th order tensor (3x3x3x3 elements)!
S = lam_0 * np.tensordot(I, I, axes=0) + 2 * mu_0 * double_tensor(I, I)

def MtoV(S):
    """Quick and Dirty method to convert a 3x3x3x3 array to Voigt notation (6x6)"""
    x = np.empty((6, 6))
    # This is really quick and dirty :D
    # And no check if the tensor is really symmetric where it should be.
    packing = {  # Indices that are used in Voigt Notation as 6×6 matrix
    11: 1111, 12: 1122, 13: 1133, 14: 1123, 15: 1113, 16: 1112,
    21: 2211, 22: 2222, 23: 2233, 24: 2223, 25: 2213, 26: 2212,
    31: 3311, 32: 3322, 33: 3333, 34: 3323, 35: 3313, 36: 3312,
    41: 2311, 42: 2322, 43: 2333, 44: 2323, 45: 2313, 46: 2312,
    51: 1311, 52: 1322, 53: 1333, 54: 1323, 55: 1313, 56: 1312,
    61: 1211, 62: 1222, 63: 1233, 64: 1223, 65: 1213, 66: 1212,
    }
    assert len(set(packing.keys())) == 36
    assert len(set(packing.values())) == 36
    
    for k, v in packing.items():
        x0, x1, i0, i1, i2, i3 = [int(x) - 1 for x in str(k) + str(v)]
        x[x0, x1] = S[i0, i1, i2, i3]
    return x
    
print("C^-1:")
print(MtoV(S))

print("-------")

print(MtoV(np.tensordot(I, I, axes=0)))

print("-------")

print(MtoV(double_tensor(I, I)))