Hello everyone,

I am using FEniCS to solve the following 1D system equation:

\begin{array}{l} \frac{{\partial \xi }}{{\partial t}} = - {{\hat L}_\sigma }\hat g\left( \xi \right) + \hat \kappa \frac{{{\partial ^2}\xi }}{{\partial t}} - {{\hat L}_\eta }\dot h\left( \xi \right)\left( {{e^{\left( {1 - \alpha } \right)\hat \phi }} - \hat c{e^{ - \alpha \hat \phi }}} \right)\\ \frac{\partial }{{\partial x}}\left[ {\hat K\left( \xi \right)\frac{{\partial \hat \phi }}{{\partial x}} + {{\hat K}_D}\left( \xi \right)\frac{1}{{\hat c}}\frac{{\partial \hat c}}{{\partial x}}} \right] = - \frac{{\partial \xi }}{{\partial t}}\\ \hat \varepsilon \left( \xi \right)\frac{{\partial \hat c}}{{\partial t}} - \frac{\partial }{{\partial x}}\left( {\hat D\left( \xi \right)\frac{{\partial \hat c}}{{\partial x}}} \right) = - \frac{{\partial \xi }}{{\partial t}} \end{array}

With the following boundaries conditions:

\begin{array}{l} \xi = 1,\,\,\,\hat \phi = - 0.05,\,\,\,\frac{{\partial \hat c}}{{\partial x}} = 0\,\,\,\,\,\,at\,\,\hat x = 0\\ \xi = 0,\,\,\,0.3\frac{{\partial \hat \phi }}{{\partial x}} + 0.4\frac{{\partial \hat c}}{{\partial x}} = 1,\,\,\,\hat c = 1\,\,\,at\,\,\hat x = 20 \end{array}

I´m trying to implement the inhomogeneus Neumann boundary condition as the manual explains and I´m not sure what I´m doing wrong.

The code regarding to the boundaries conditions is the next :

```
# %%Define boundary condition
lox = 20
x = sym.symbols('x[0]')
xi, c, phi = u.split()
# g = dphi/dx = 1/0.3 - (0.4/0.3)*dc/dx
g = Constant(1/0.3) - Constant(0.4/0.3)*sym.diff(c, x).subs(x,lox)
# Define boundary conditions
boundary_conditions = {1: {'Neumann': g}} # x = 1
# Define boundary subdomains
tol = 1e-14
class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], lox, tol)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim())
#boundary_markers.set_all(9999)
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)
#Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v_2*ds(i))
def boundary0(x, on_boundary):
return on_boundary and near(x[0], 0)
def boundaryL(x, on_boundary):
return on_boundary and near(x[0], lox)
bc_xi1 = DirichletBC(V.sub(0), Constant(chi_ini), boundary0)
bc_xi2 = DirichletBC(V.sub(0), Constant(chi_fin), boundaryL)
bc_c2 = DirichletBC(V.sub(1), Constant(c_fin), boundaryL)
bc_phi1 = DirichletBC(V.sub(2), Constant(phi_ini), boundary0)
bcs = [bc_xi1, bc_xi2, bc_c2, bc_phi1] # Condiciones dirichlet
################################################################
```

The message error is the following:

File “”, line 1, in

runfile(’/home/josemanuel/Documents/Fenics/Newmann1D.py’, wdir=’/home/josemanuel/Documents/Fenics’)

File “/usr/lib/python3/dist-packages/spyder_kernels/customize/spydercustomize.py”, line 827, in runfile

execfile(filename, namespace)

File “/usr/lib/python3/dist-packages/spyder_kernels/customize/spydercustomize.py”, line 110, in execfile

exec(compile(f.read(), filename, ‘exec’), namespace)

File “/home/josemanuel/Documents/Fenics/Newmann1D.py”, line 154, in

g = Constant(1/0.3) - Constant(0.4/0.3)*sym.diff(c, x).subs(x,lox)

File “/usr/lib/python3/dist-packages/sympy/core/function.py”, line 2448, in diff

return Derivative(f, *symbols, **kwargs)

File “/usr/lib/python3/dist-packages/sympy/core/function.py”, line 1243, in **new**

expr = sympify(expr)

File “/usr/lib/python3/dist-packages/sympy/core/sympify.py”, line 328, in sympify

coerced = coerce(a)

File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 468, in **float**

raise RuntimeError(“Cannot convert spatially varying function to float.”)

RuntimeError: Cannot convert spatially varying function to float.

Anyone who can help me?

Thank you in advance.