Use of vector value in final form to be solved for solving a 3D equation

variational formulation
1

I am currently trying so solve a simple time dependent 3D heat equation, however I am currently encountering some issues using a vector quantity in the final form to solve. My final form to solve looks like this:

    G = as_tensor(T.dx(i), (i))  #gradient of T
    G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step 
        
    q = -kappa*G

    F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv

But there seems to be an issue with using the vector quantity q as I get the following error:

Traceback (most recent call last):
  File "/Users/parameterspace-test.py", line 79, in <module>
    F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

How do I resolve this issue? My code is below.

from __future__ import print_function
import numpy as np
import scipy as scipy
import matplotlib.pyplot as plt
from dolfin import *
import ufl
from ufl import as_tensor
from ufl import Index
import math

parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
set_log_level(50)

#----------------
t_start = 0.0
t_end = 200
nstep = 400
dtime = (t_end - t_start)/ nstep  #time step
p = 280 #W
#---------------------------------------

kappamin = 0
kappamax = 1000
npstep = 100

for b in np.linspace(kappamin, kappamax, npstep):
    kappa = b

    data = 'mesh_3mm'
    mesh = Mesh()
    hdf = HDF5File(mesh.mpi_comm(), data+'.h5', 'r')
    hdf.read(mesh, '/mesh', False)
    cells_f = MeshFunction('size_t', mesh, 3)
    hdf.read(cells_f, '/cells')
    facets_f = MeshFunction('size_t', mesh, 2)
    hdf.read(facets_f, '/facets')
    boundaries = facets_f

    all_cells = np.array(list(cells(mesh)))

    tdim = mesh.topology().dim()
    fdim = tdim - 1
    mesh.init(fdim, tdim)
    f_to_c = mesh.topology()(fdim, tdim)
    c_to_f = mesh.topology()(tdim, fdim)

    #------------
    i,j,k = ufl.indices(3)
    n = FacetNormal(mesh)

    #------------------------------------------------

    VectorSpace = VectorFunctionSpace(mesh, 'P', 1)
    da = Measure('ds', domain=mesh, subdomain_data = facets_f, metadata = {'quadrature_degree': 2})  #area element

    dv = Measure('dx', domain=mesh, subdomain_data = cells_f, metadata = {'quadrature_degree': 2})   #volume element

    Space = FunctionSpace(mesh, 'P', 1) 
    T = Function(Space)
    T0 = Function(Space)
    T_init = Expression('Tambient', degree=1, Tambient=600.)
    T = project(T_init, Space)
    assign(T0,T)

#-------------------------------------------

    V = TestFunction(Space)     # Test Function used for FEA
    dT = TrialFunction(Space)   # Temperature Derivative
    q0 = Function(VectorSpace)  # heat flux at previous time step
    i = Index()
    G = as_tensor(T.dx(i), (i))  #gradient of T
    G0 = as_tensor(T0.dx(i), (i)) # gradient of T at previous time step 
        
    q = -kappa*G

    F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv 
    
    Gain = derivative(F, T, dT)    # Gain will be used as the Jacobian required to determine the evolution of a dynamic system 


    problem = NonlinearVariationalProblem(F, T, [], Gain)
    solver  = NonlinearVariationalSolver(problem)

    solver.parameters["newton_solver"]["relative_tolerance"] = 1E-4
    solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-3
    solver.parameters["newton_solver"]["convergence_criterion"] = "residual"
    solver.parameters["newton_solver"]["error_on_nonconvergence"] = True
    solver.parameters["newton_solver"]["linear_solver"] = "cg"
    solver.parameters["newton_solver"]["maximum_iterations"] = 10000
    solver.parameters["newton_solver"]["preconditioner"] = "hypre_euclid"
    solver.parameters["newton_solver"]["report"] = True

    solver.parameters["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = False
    solver.parameters["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-4
    solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-3
    solver.parameters["newton_solver"]["krylov_solver"]["monitor_convergence"] = False
    solver.parameters["newton_solver"]["krylov_solver"]["maximum_iterations"] = 100000

    t = 0

    #file_T << (T,t)
    is_looping = True
    for t in np.arange(t_start + dtime,t_end + dtime,dtime):
        print( "Time", t)
        solver.solve()

        assign(T0, T)

    else:
        is_looping = False
        break
    if not is_looping:
        break
2

Please note that your code cannot reproduce the error, as you have not supplied the mesh.
A simple fix to make it reproducible would be to use:

mesh = UnitCubeMesh(10, 10, 10)
cells_f = MeshFunction('size_t', mesh, 3)
facets_f = MeshFunction('size_t', mesh, 2)

Secondly, I would strongly advice you to read in meshes outside your for-loop, as this is very expensive (and unnecessary).

Finally, to resolve your issue simply doing:

    kappa = Constant(b)
    # .. other code
    q = -kappa*G
    F = (1/dtime*(T-T0)*V - q[i]*V.dx(i)) * dv

resolves the issue