Error: cannot replace unknown element domain without unique common domain

Question

Hi there,
I followed the Cahn-Hilliard tutorial to solve a system of two PDEs, but met some problem. I don't quite understand this error message. Could you please help me find which part of the code is wrong? The main difference from the tutorial is that I use two assign commands to initialise function. Is it correct here?
I post the complete code here:

import random
from dolfin import *

# Class representing the intial conditions
#  initial conditon for unknown u
class InitialConditions_u(Expression):
    def eval(self, values, x):
        values[0] = 2*x[0]
        values[1] = 2*x[1]
    def value_shape(self):
        return (2,)
 # initial condition for unknown p
class InitialConditions_p(Expression):
    def eval(self, values, x):
        values[0] = 1 + x[0]*x[0] + x[1]*x[1]

# Class for interfacing with the Newton solver
class MyNonlinear(NonlinearProblem):
    def __init__(self, L, a, bc):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc= bc
    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bc.apply(b,x)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        self.bc.apply(A)

# Model parameters
dt     = 5.0e-06  # time step
t      = 0

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"

# Create mesh and define function spaces
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
Q = FunctionSpace(mesh, "DG", 0)
W = V*Q

# Define trial and test functions
dw    = TrialFunction(W)
v, q  = TestFunctions(W)


# functions to initialise w and w0
ui = Function(V)
pi = Function(Q)

# Create intial conditions and interpolate
u_init = InitialConditions_u()
ui.interpolate(u_init)

# Create intial conditions and interpolate
p_init = InitialConditions_p()
pi.interpolate(p_init)

# Define functions
w   = Function(W)  # current solution
w0  = Function(W)  # solution from previous converged step

# assign and split
assign(w,[ui,pi])
assign(w0,[ui,pi])
du, dq  = split(dw)
(u,p) = w.split()
(u0,p0) = w0.split()

# plot subclass for test
plot(p0,title ="Initial pressure")
plot(u0,title ="Initial velocity")
interactive()

# define boudary
def boundaries(x, on_boundary):
    return on_boundary
# define normal vector on boundary
normal = FacetNormal(mesh)

# define p expression on boundary (essential boundary condtion)
u_boundary_expression = Expression("1 + x[0]*x[0] + x[1]*x[1] + t", t = t)

# define boundary condition
bc = DirichletBC(W.sub(1), u_boundary_expression, boundaries)

# define f expression
f = Expression("5 + 8*x[0]*x[0] + 8*x[1]*x[1] + 4*t", t = t)

# Weak statement of the equations
L0 = p*q*dx - p0*q*dx - dt*p*inner(u, nabla_grad(q))*dx - f*q*dx
L1 = inner(u, v)*dx + p*inner(v, normal)*ds - p*div(v)*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, w, dw)

# Create nonlinear problem and Newton solver
problem = MyNonlinear(a, L, bc)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("output.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    f.t = t
    w0.vector()[:] = w.vector()
    solver.solve(problem, w.vector())
    file << (w.split()[0], t)

plot(w.split()[0])
interactive()

And the error message I got is:

Cannot replace unknown element domain without unique common domain in form.
Traceback (most recent call last):
File "MyProblem_updated.py", line 119, in <module>
  solver.solve(problem, w.vector())
File "MyProblem_updated.py", line 25, in F
  assemble(self.L, tensor=b)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 230, in assemble
  cell_domains, exterior_facet_domains, interior_facet_domains)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 80, in _create_dolfin_form
  mesh=mesh)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 98, in __init__
  = jit(form, form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 62, in mpi_jit
  return local_jit(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 126, in jit
  return jit_compile(form, parameters=p)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 77, in jit
  return jit_form(ufl_object, parameters)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 103, in jit_form
  form_data = form.compute_form_data()
File "/usr/lib/python2.7/dist-packages/ufl/form.py", line 169, in compute_form_data
  self._form_data = preprocess(self, object_names=object_names)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/preprocess.py", line 188, in preprocess
  element_mapping = _compute_element_mapping(extract_elements(original_form), common_domain)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/preprocess.py", line 348, in _compute_element_mapping
  "Cannot replace unknown element domain without unique common domain in form.")
File "/usr/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert
  if not condition: error(*message)
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 151, in error
  raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot replace unknown element domain without unique common domain in form.

Answer 1

Seems like you have mixed a and L in MyNonlinear. And try

u, p = split(w)    # not u, p = w.split()
u0, p0 = split(w0) # not u0, p0 = w0.split()

Comments

Yes. It's a mixed function space of one vector and one scalar. Does this influence the way to define the nonlinear problem?

---

I mean, MyNonlinear takes (L, a, bc) in that order in the signature and you are creating in with (a, L, bc). Could be the problem...

---

Thank you very much. It works now.