population of elasticity tensor

I am trying to create a code to populate a tensor with the symmetries of an elasticity tensor which should have 36 independent components and the symmetries C[i,j,k,l] = C[j,i,k,l] = C[i,j,l,k].

I am doing:

C = sym.MutableDenseNDimArray.zeros(3,3,3,3)
var = sym.symbols('C1:%d'%82)

ct = 0
for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                if C[i,j,k,l] == 0:
                    C[i,j,k,l] = var[ct]
                    C[j,i,k,l] = var[ct]
                    C[i,j,l,k] = var[ct]
                    ct += 1

This is giving me 45 different entries. What am I doing wrong?

You seem to be missing a symmetry. If the 2 symmetries you mentioned:

C[i,j,k,l] = C[j,i,k,l] (symmetry 1)
C[i,j,k,l] = C[i,j,l,k] (symmetry 2)

hold for every element that means that:

C[a,b,c,d] = C[b,a,c,d] (symmetry 1)
C[b,a,c,d] = C[b,a,d,c] (symmetry 2)

and therefore:

C[a,b,c,d] = C[b,a,d,c] or
C[i,j,k,l] = C[j,i,l,k] (implied symmetry)

If you include that in your code you'll get:

C = sym.MutableDenseNDimArray.zeros(3,3,3,3)
var = sym.symbols('C1:%d'%82)

ct = 0
for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                if C[i,j,k,l] == 0:
                    C[i,j,k,l] = var[ct]
                    C[j,i,k,l] = var[ct]
                    C[i,j,l,k] = var[ct]
                    C[j,i,l,k] = var[ct]
                    ct += 1

And that should give you 36 different entries.