The use of assemble_system with multiple subdomains

Question

Hello!

I was trying to use assemble_system at the same time as using multiple subdomains and I was getting the following error:

*** Warning: Bilinear and linear forms do not have same cell_domains subdomains in SystemAssembler. Taking cell_domains subdomains from bilinear form

so I took the demo_subdomains-poisson.py that FEniCS offers, modifying just the way of solving the system, to do it now with assemble_system,

# Define new measures associated with the interior domains and exterior boundaries
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]

# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
 - g_L*v*ds(1) - g_R*v*ds(3)
 - f*v*dx(0) - f*v*dx(1))

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
A = assemble(a)
b = assemble(L)
A, b = assemble_system(a, L, bcs)
u = Function(V)
solve(A, u.vector(), b)

and I still get the same message. Why is that not possible?

Thank you in advance!

Comments

Hi, I assume you are using DOLFIN version 1.3.0. There, the way to call assemble_system which yields the same results as the original demo is

A, b = assemble_system(a, L, bcs, cell_domains=domains,
                       exterior_facet_domains=boundaries)

There is no such warning with version 1.4.0+ and the assembler picks up correct cell/facet domains.

---

I am using DOLFIN version 1.4.0, is it the correct way to call assemble_system in such a version as I did?

---

I get that warning you described with 1.4.0. If you solve $AU=b$ with $A, b$ from assemble_system, how does the solution compare with the one from solve(a == L, ...)?

---

Using solve(A, u.vector(), b) I get exactly the same warning as using solve(a == L, u, bcs), both with A, b = assemble_system(a, L, bcs).

---

Is this an error from FEniCS that is not possible to use assemble_system at the same time as multiple subdomains?

---

It is a warning and checking the difference of solution u from solve(a == L, u, bcs) and
uu from solve(A, uu.vector(), b) where $A, b$ come from assemble_system you could see how serious the warning is. If they differ, assemble_system is probably not picking up cell/facet domains correctly.

---

With both u from solve(a == L, u, bcs) and uu from solve(A, uu.vector(), b) I am getting the same solution, 0, so I suppose the problem is not with assemble_system. Any idea what it could be?

Answer 1

I think this is fixed in dev. Here's a possible workaround:

a, L = ...
# replace 99 with any unused boundary id
a += u*v*dx(99) + u*v*ds(99)
L += v*ds(99) + v*ds(99)

This makes both the lhs and rhs form contain the same cell and facet functions. It shouldn't contribute to the result or the assembly computation time.

Let me know if this works.