Subdomains with different coefficient dependent on solution

Hi, guys.
I am currently stuck at the following problem.

I want to define a coefficient dependent on solution and subdomains. But it produced an error: setting an array element with a sequence.
If i don’t use different coefficient on subdomains, it works. Also, if the coefficient doesn’t dependent on the solution, it works.

Here is part of my code related with the error:

from fenics import *


mesh = UnitSquareMesh(5,5);
U = FunctionSpace(mesh, “Lagrange”, 1);

class Air(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0 + DOLFIN_EPS

class Glass(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= 0 - DOLFIN_EPS

subdomains = MeshFunction(“size_t”, mesh, mesh.topology().dim(),0)

air = Air().mark(subdomains,0)
glass = Glass().mark(subdomains, 1)

u = Function(U);
v =TestFunction(U);

class KK(UserExpression):
def init(self, subdomains, k_0, k_1, **kwargs):
self.subdomains = subdomains
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 0:
values[0] = self.k_0
values[0] = self.k_1
def value_shape(self):
return ()

Wpii_g=sigmaK_g * inner(u,u)**2
Wpii_a=sigmaK_a * inner(u,u)**7
Wpii = KK(subdomains=subdomains, k_0 = Wpii_a, k_1 = Wpii_g, degree=1)

F=Wpii * v * u * dx;

solve(F == 0, u)

And Here is my error code:

ValueErrorTraceback (most recent call last)
46 F=Wpii * v * u * dx;
—> 48 solve(F == 0, u)

/usr/local/lib/python3.6/dist-packages/dolfin/fem/ in solve(*args, **kwargs)
218 # tolerance)
219 elif isinstance(args[0], ufl.classes.Equation):
–> 220 _solve_varproblem(*args, **kwargs)
222 # Default case, just call the wrapped C++ solve function

/usr/local/lib/python3.6/dist-packages/dolfin/fem/ in _solve_varproblem(*args, **kwargs)
264 solver = NonlinearVariationalSolver(problem)
265 solver.parameters.update(solver_parameters)
–> 266 solver.solve()

/usr/local/lib/python3.6/dist-packages/dolfin/function/ in wrapped_eval_cell(self, values, x, cell)
55 def wrapped_eval_cell(self, values, x, cell):
—> 56 self.user_expression.eval_cell(values, x, cell)
58 # Attach user-provided Python eval functions (if they exist in

in eval_cell(self, values, x, cell)
36 values[0] = self.k_0
37 else:
—> 38 values[0] = self.k_1
39 def value_shape(self):
40 return ()

ValueError: setting an array element with a sequence.

I found this error is related with shape of np array. But i can’t understand why it happened.
Is there any way to define the different nonlinear coefficient in different subdomains?

Thank you for taking the time to help me out!


I assume that k0 and k1 are just fenics Functions. What you can do, to make it work, is to actually evaluate these functions, i.e., in your eval_cell routine of KK, instead of

values[0] = self.k_0


values[0] = self.k_0(x)