Division operator does not work to substitute a function

[Bourceanu Cristian]

I am solving a diffusion equation that was initially not linear and I have done a substitution to get a linear form. The solver throws results that I was expecting, but when I want to back subtitute to get the form of the original distribution, I get the following error:

‌

<code>

File "/home/cristi/Documents/Scripts and Projects of IDEs/Python/Internship/FEniCS_myModels/Jing Liu Paper/dimensionless_model_V3.py", line 239, in solver
vtkfile << (elem_div(u*(1+alpha),(1+alpha*u)), t)
File "/home/cristi/.conda/envs/fenicsproject/lib/python3.7/site-packages/dolfin/io/init.py", line 19, in lshift
self.write(u[0], u[1])
AttributeError: 'Division' object has no attribute '_cpp_object'

</code>

‌

The code:

# Create mesh and define function space
mesh = UnitIntervalMesh(num_intervals)
V = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement(V, R))
VS = FunctionSpace(mesh, V)

(u, c) = TrialFunction(W)
(v, d) = TestFunction(W)

# Initial Condition
u_past = interpolate(Constant(u_init_val),VS)
u_avg = Constant(u_init_val)

#Boundary Conditions
(bcs, integrals_N, integrals_R, bcE) = boundaryConditions(
    u, v, V=VS, mesh=mesh)
integBound = sum(integrals_N)+sum(integrals_R)

def tDeriv():
    return (u-u_past)*v/dt
def Laplacian():
    return -(inner(grad(u),grad(v))+c*v+(u-u_avg)*d)
def FaradaicConv():
    return 2*beta*i_electrodes*grad(u)[0]*v

# Define variational problem
g = Expression("x[0]<0.5 ? 1 : -1", degree=1)
F = (-Laplacian()+FaradaicConv())*dx+integBound
w = Function(W)
limiting_current(F,w,bcs,max_conc=10)
F = F+tDeriv()*dx
a,L=lhs(F),rhs(F)


# Make VTK format file to save the solution
vtkfile = File('diffusion_toyProblem/solution.pvd')

# Create time series
timeseries_u = TimeSeries('diffusion_reaction/concentration')
timeseries_u.store(u_past.vector(), 0)

# Save mesh to file
File('diffusion_reaction/mesh.xml.gz') << mesh

t = 0
for n in range(time_steps):
    t += dt

    solve(a==L, w, bcs)

    (u, c) = w.split(True)

    u_past.assign(u)

    updateTimeBoundary(bcE, t)

    vtkfile << (elem_div(u*(1+alpha),(1+alpha*u)), t)

    timeseries_u.store(u_past.vector(), t)

‌

I have also tried the following expression, but the same error occurs:

<code class=”python”>

vtkfile<<(u*(1+alpha)/(1+alpha*u),t)

</code>