How to create a Form without using deprecated free variable notation?

variational formulation
1

I am reading through Computational Reality by B. Emek Abali and a code sample contains the following script (with minor modifications on my part in the form of adding a degree argument to each of the Expressions):



from fenics import *
import numpy
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
set_log_level(False)

rho = 7860.0
c = 624.0
kappa = 30.1
h = 18.0
Ta = 300.0
l = 0.1
thickness = 0.02
P = 3.0e6
speed = 0.02

# Create mesh and function space
mesh = BoxMesh(Point(0.0,0.0,0.0), Point(1,1,thickness), 200, 200, 2)
Space = FunctionSpace(mesh, 'P', 1)
cells = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
da = Measure('ds', domain=mesh, subdomain_data=facets)
dv = Measure('dx', domain=mesh, subdomain_data=cells)
t = 0.0
t_end = 50.0
Dt = 0.1
initial_T = Expression("Tini", Tini=Ta, degree=1)
T0 = interpolate(initial_T, Space)
Laser = Expression("P*exp(-50000.0*(pow(x[0]-0.5*l*(1+0.5*sin(2*pi*t/t_e)), 2)+pow(x[1]-velo*t, 2)))",\
                    P=P, t=0, t_e=t_end/10.,l=l, velo=speed, degree=2)
T = TrialFunction(Space)
del_T = TestFunction(Space)
Form = (rho*c/Dt*(T-T0)*del_T\
        +kappa*T.dx(i)*del_T.dx(i)\
        -rho*Laser*del_T)*dv\
        +h*(T-Ta)*del_T*da

left = lhs(Form)
right = rhs(Form)

A = assemble(left)
b = None
T = Function(Space)
file_T = File("Solutions/laser_on_steel.pvd")

for t in numpy.arange(0, t_end, Dt):
    print("Time ", t)
    Laser.t = t
    b = assemble(right, tensor=b)
    solve(A, T.vector(), b, 'cg')
    if t == int(t): file_T << (T, t)
    T0.assign(T)

Unfortunately, I’m returned the following with regards to the definition of the Form variable ~ line 35:

Ordering mesh.
Elapsed wall, usr, sys time: 0.0295507, 0.02, 0 (Build BoxMesh)
Elapsed wall, usr, sys time: 1.353e-06, 0, 0 (Number distributed mesh entities)
Elapsed wall, usr, sys time: 0.067164, 0.06, 0.01 (Init dofmap from UFC dofmap)
Computing mesh entities of dimension 2
Elapsed wall, usr, sys time: 0.370391, 0.34, 0.03 (Compute entities dim = 2)
Requesting connectivity 2 - 3.
Requesting connectivity 3 - 2.
Computing mesh connectivity 2 - 3 from transpose.
Elapsed wall, usr, sys time: 0.0257503, 0.02, 0 (Compute connectivity 2-3)
Determining node ownership for parallel dof map
Finished determining dof ownership for parallel dof map
Elapsed wall, usr, sys time: 0.000112316, 0, 0 (SCOTCH: call SCOTCH_graphBuild)
Elapsed wall, usr, sys time: 0.0130756, 0.01, 0 (SCOTCH: call SCOTCH_graphOrder)
Elapsed wall, usr, sys time: 0.0159682, 0.02, 0 (Compute SCOTCH graph re-ordering)
Elapsed wall, usr, sys time: 0.703672, 0.66, 0.04 (Init dofmap)
Elapsed wall, usr, sys time: 1.8163e-05, 0, 0 (Apply (PETScVector))
Elapsed wall, usr, sys time: 1.958e-06, 0, 0 (Apply (PETScVector))
Elapsed wall, usr, sys time: 0.00445111, 0.01, 0 (Init dof vector)
Elapsed wall, usr, sys time: 4.072e-06, 0, 0 (Apply (PETScVector))
Elapsed wall, usr, sys time: 1.0877e-05, 0, 0 (Apply (PETScVector))
Traceback (most recent call last):
  File "/mnt/leroux/p_laufwerk/A133_HIWIS/Simulationen/Fenics/fenics_tests.py", line 36, in <module>
    +kappa*T.dx(i)*del_T.dx(i)\
NameError: name 'i' is not defined. Did you mean: 'io'?

It seems the notation of using Method.dx(i) as a free variable is no longer supported, however, I can’t seem to find a solution online that keeps the functionality that this version of the code would have, unless I am misunderstanding.

I appreciate any help on this issue, thanks.

2

You seem to be missing a definition i = Indices(1). You cannot have an undefined variable in Python in general, so I have no clue where you got this code from.

3

Hi Dokken, thanks for responding. As I said, it’s directly from Computational Reality, on page 119, to be precise. I’m brand new to fenics and have been trying to parse it out, though the solution you gave isn’t a defined name in itself. Where does Indices() come from?

NameError: name 'Indices' is not defined

Thanks

4

I dont have that book. Indices can be imported from ufl, ie.

from ufl import indices
i = indices(1)