Hello

First of all, I installed FEniCS on a Windows 10 subsystem for Linux with an Ubuntu terminal and I am working with version 2019.1.0. I am attempting to test the error of a fully discretised mixed finite element scheme for a transformed and linearised form of Richards Equation. The exact solution against which the error is measured is:

u_{\text{ex}}(x, y, t) = \begin{cases}-\frac{2(e^{s}-1)}{e^{s}+1},\qquad\text{for } s\geq0,\\ -s,\qquad\ \qquad\text{for } s<0,\end{cases}

with s=x-y-t.

The variational form used in the scheme is as follows, find \left(u^{n,i}_h, \vec{q^{n,i}_h}\right)\in W_h\times V_h such that:

k\left(u^{n,i}_h, w_h \right)_{L^2} + \tau\left(\nabla\cdot\vec{q^{n,i}_h}, w_h\right)_{L^2}=\left(b\left(u^{n-1}_h\right)-b\left(u^{n, i-1}_h\right)+ ku^{n, i-1}_h, w_h\right)_{L^2},\quad \forall w_h\in{W_h},

\left(\vec{q^{n, i}_h}, \vec{v_h}\right)_{L^2} - \left(u^{n,i}_h,\nabla\cdot \vec{v_h} \right)_{L^2} = 0, \qquad \forall \vec{v_h}\in V_h.

Here W_h and V_h are the test spaces which are finite subspaces of L^2\left(\Omega\right) and H\left(\text{div}, \Omega\right) respectively. The superscript n denotes the time step which is being considered and the superscript i indicates the stage within the iterative Newton linearisation procedure performed at each time step n. The value k is a constant involved in the Newton linearisation. Also (\cdot,\cdot)_{L^2} is the L^2 inner product and \Omega=[0,1]^2.

The issue I am having relates to the fact that the function b is a conditional function of u, specifically:

b\left(u\right)=\begin{cases}\frac{\pi^2}{2}-\frac{u^2}{2}, \qquad\text{where }u\leq 0,\\\frac{\pi^2}{2},\qquad\qquad \ \text{where } u>0.\end{cases}

I have attempted to express this function using the `conditional`

command but I either generate `nan`

solutions at each time step or I receive the error message:

`NotImplementedError: Cannot take length of non-vector expression`

.

I provide the working examples for each case below.

# Generation of nan solutions:

Running this will show that the numerical solution `u`

has value `nan`

at each vertex on the mesh.

```
from fenics import *
from dolfin import *
from mshr import *
import numpy as np
#*****************range of mesh refinements and step totals********************
dofs = [4, 8, 16, 32]
steps = [25, 50, 100, 200]
#********solving numerically at each time step and mesh size level************
for j in range(len(dofs)):
#******************* Setting up unit square mesh ***************************
mesh = UnitSquareMesh(dofs[j], dofs[j])
#***************Defining variables and corresponding function spaces**********
W = FiniteElement('DG', mesh.ufl_cell(), 0) #Trial/test space of Discontinuous Lagrange elements
V = FiniteElement('RT', mesh.ufl_cell(), 1) #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element) #Creating mixed finite element space
W_0 = M.sub(0).collapse() #Rebuilding trial/test space of constant elements
V_0 = M.sub(1).collapse()
#*******************Exact solution and boundary condition**********************
u_Ex = Expression('x[0]-x[1]-t >= 0 ? -2*(exp(x[0]-x[1]-t)-1)/(exp(x[0]-x[1]-t)-1): -(x[0]-x[1]-t)',
degree = 3, t = 0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(M.sub(0), u_Ex, boundary)
#*******************Initial Condition******************************************
u_ = project(u_Ex, W_0)
#*****************Setting values of certain parameters************************
k = 2.0 # Linearisation constant
#************** Fixing duration of evolution and timestep size****************
T = 1 #Final time
N = steps[j] #Number of steps
tau = T/N #Time step size
I = 8 #Number of interim steps in linearization method
#*****************Defining variational form************************************
u_i_ = u_ #Taking initial value within linearisation iteration as value from previous timestep
pi = np.pi
u_i, q_i = TrialFunctions(M) #Defining trial and test functions of the mixed formulation
w, v = TestFunctions(M)
bn_ = conditional(le(0, u_), pow(pi,2)/2 - pow(u_,2)/2, pow(pi,2)/2) #The function b evaluated at the value of u from the previous time step n-1
bni_ = conditional(le(0, u_i_), pow(pi,2)/2 - pow(u_i_,2)/2, pow(pi,2)/2) #The function b evaluated at the value of u from the previous step of the nested iteration (n, i-1)
a1 = k*(u_i*w*dx) + tau*((div(q_i))*w*dx)
L1 = (bn_- bni_ + k*u_i_)*w*dx
a2 = dot(q_i, v)*dx - u_i*(div(v))*dx
L2 = 0
a = a1 + a2
L = L1 + L2
#****************************Advancing scheme in time*************************
t = 0 #Initial time
uq = Function(M) #Solution
for n in range(N):
t += tau #Updating current time
u_Ex.t = t #Updating boundary condition and exact solution
u_Ex_nodes = u_Ex.compute_vertex_values(mesh) #Extracting vertex values of the exact solution
print(u_Ex_nodes)
for i in range(I): #Linearisation iteration
print(i)
solve(a==L, uq, bc)
u_i, q_i = uq.split(True)
u_i_.assign(u_i) #assigning the current solution to be the previous solution in the next run of the loop
u_.assign(u_i_) #once iteration is finished the value of u at this time is recorded and entered as the previous value in the next run of the time loop.
u_nodes = u_.compute_vertex_values(mesh) #storing vertex values of numerical solution at this time step
print(u_nodes)
```

# NotImplementedError: Cannot take length of non-vector expression:

In this case the error message appears upon definition of the variational form.

```
from fenics import *
from dolfin import *
from mshr import *
import numpy as np
#*****************range of mesh refinements and step totals********************
dofs = [4, 8, 16, 32]
steps = [25, 50, 100, 200]
#********solving numerically at each time step and mesh size level************
for j in range(len(dofs)):
#******************* Setting up unit square mesh ***************************
mesh = UnitSquareMesh(dofs[j], dofs[j])
#***************Defining variables and corresponding function spaces**********
W = FiniteElement('DG', mesh.ufl_cell(), 0) #Trial/test space of Discontinuous Lagrange elements
V = FiniteElement('RT', mesh.ufl_cell(), 1) #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element) #Creating mixed finite element space
W_0 = M.sub(0).collapse() #Rebuilding trial/test space of constant elements
V_0 = M.sub(1).collapse()
#*******************Exact solution and boundary condition**********************
u_Ex = Expression('x[0]-x[1]-t >= 0 ? -2*(exp(x[0]-x[1]-t)-1)/(exp(x[0]-x[1]-t)-1): -(x[0]-x[1]-t)',
degree = 3, t = 0)
print(u_Ex)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(M.sub(0), u_Ex, boundary)
#*******************Initial Condition******************************************
u_ = project(u_Ex, W_0)
#*****************Setting values of certain parameters************************
k = 2.0 # Linearisation constant
#************** Fixing duration of evolution and timestep size****************
T = 1 #Final time
N = steps[j] #Number of steps
tau = T/N #Time step size
I = 8 #Number of interim steps in linearization method
#*****************Defining variational form************************************
u_i_ = u_ #Taking initial value within linearisation iteration as value from previous timestep
pi = np.pi
u_i, q_i = TrialFunctions(M) #Defining trial and test functions of the mixed formulation
w, v = TestFunctions(M)
bn_ = conditional(le(0, u_), pow(pi,2)/2 - pow(u_,2)/2, pow(pi,2)/2) #The function b evaluated at the value of u from the previous time step n-1
bni_ = conditional(le(0, u_i_), pow(pi,2)/2 - pow(u_i_,2)/2, pow(pi,2)/2) #The function b evaluated at the value of u from the previous step of the nested iteration (n, i-1)
a1 = k*(u_i*w*dx) + tau*((div(q_i))*w*dx)
L1 = (bn_(u_)- bni_(u_i_) + k*u_i_)*w*dx
a2 = dot(q_i, v)*dx - u_i*(div(v))*dx
L2 = 0
a = a1 + a2
L = L1 + L2
```

This produces the error:

```
Traceback (most recent call last):
File "r_k_m_c_conditional_mwe_notimplerr.py", line 84, in <module>
L1 = (bn_(u_)- bni_(u_i_) + k*u_i_)*w*dx
File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 341, in _call
return _eval(self, arg, mapping, component)
File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 331, in _eval
return f.evaluate(coord, mapping, component, index_values)
File "/usr/lib/python3/dist-packages/ufl/conditional.py", line 250, in evaluate
c = self.ufl_operands[0].evaluate(x, mapping, component, index_values)
File "/usr/lib/python3/dist-packages/ufl/conditional.py", line 134, in evaluate
b = self.ufl_operands[1].evaluate(x, mapping, component, index_values)
File "/usr/lib/python3/dist-packages/ufl/core/terminal.py", line 71, in evaluate
return self.ufl_evaluate(x, component, derivatives)
File "/usr/lib/python3/dist-packages/dolfin/function/function.py", line 270, in ufl_evaluate
return self(*x)
File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 406, in __len__
raise NotImplementedError("Cannot take length of non-vector expression.") NotImplementedError: Cannot take length of non-vector expression.
```

I accept this is a large post but I would really appreciate any help at all so that I can finally incorporate this conditional expression into my variational form.

Thank you