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

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')
cells_f = MeshFunction('size_t', mesh, 3)
facets_f = MeshFunction('size_t', mesh, 2)
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

``````

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