How to derive 4th order tensor component by component at each node and assemble it to form a tensorfunction in the mesh?

Hi all,
I am now trying to solve a nonlinear hyperelasticity problem with complicated constitutive relations. So I have to derive the Jacobian of F myself, and take derivative of stress S with respect to strain E manually to find each component of the 4th order tensor C at each node and. I have shown how I compute each component of C at each node in the following code. The question is how to assemble each component of C to a 4th order tensor and make it a tensor function throughout the whole mesh?
Any help would be appreciated!

#u is the displacement field

#deformation gradient 
def Fmat(u):
  d = u.geometric_dimension()
  I = Identity(d)             
  F = I + grad(u)            
  return F

# strain tensor E
def Emat(u):
  d = u.geometric_dimension()
  I = Identity(d)             
  F = Fmat(u) 
  return 0.5*(F.T*F - I)


# fourth order tensor Cijkl
def C_fourth_order(u, mesh):
    E = Emat(u)
    E11 = project(E[0,0],FunctionSpace(mesh, "DG", 0))
    E11 = E11.vector().array()

    E22=project(E[1,1],FunctionSpace(mesh, "DG", 0))
    E22 = E22.vector().array()

    E33=project(E[2,2],FunctionSpace(mesh, "DG", 0))
    E33=E33.vector().array()

    E12 = project(E[0,1],FunctionSpace(mesh, "DG", 0))
    E12 = E12.vector().array()

    E13 = project(E[0,2],FunctionSpace(mesh, "DG", 0))
    E13 = E13.vector().array()

    E23 = project(E[1,2],FunctionSpace(mesh, "DG", 0))
    E23 = E23.vector().array()

    # weight and points are used for gauss integration
    weight = [0.5688888888888889, 0.4786286704993665,0.4786286704993665,0.2369268850561891,0.3478548451374538]
    points = [0.0000000000000000,-0.5384693101056831,0.5384693101056831,-0.9061798459386640,0.9061798459386640]


    # example of how I compute component C1111= dS11/dE11 of 4th order tensor Cijkl.
    # Final11() is a function of integration, it needs components of strain Eij at each node to compute Cijkl at that node 
    # I am not clever, so I just use a list to store the component C1111 at all the nodes  
    C_1111 = []
    for n  in range(0,len(E11)):
        total = 0
        for i in range(0,len(weight)):
            for j in range(0,len(weight)):
                for k in range(0,len(weight)):
                       total = total + pi/24*pi*pi/2*weight[i]*weight[j]*weight[k]*Final11(points[i],points[j],points[k],E11[n],E22[n],E33[n],E12[n],E13[n],E23[n],'11')
        C_1111.append(total)

So I just stuck here, I can compute the other components of Cijkl correspondingly, but how to form a tensor function is a problem for me.
Any help would be greatly appreciated!!!