evaluate drive term f on submesh

Question

Hello all,

I need to solve a PDE (let us say a Poisson equation) on a domain $\Omega$ with a forcing term f that is non-zero only in a subdomain $\Omega_0\subset \Omega$. (Both $\Omega$ and) $\Omega_0$ cannot be described by a simple mathematical function. I have therefore meshed $\Omega$ with gmsh and made sure that the mesh is aligned with $\Omega_0$.
I thought of extracting a submesh from the mesh

mesh = Mesh("mymesh.xml")
subdomains = MeshFunction("size_t", mesh, "mymesh_physical_region.xml")
submesh = SubMesh(mesh, subdomains, 1)
V = FunctionSpace(mesh, 'Lagrange', 2)
u, v = TrialFunction(V),  TestFunction(V)
Vsub = FunctionSpace(submesh, 'Lagrange', 2)

and then of defining a function on $\Omega_0$ but the following is not a valid syntax

f = Function(Vsub, "pow(x[0], 2)")

(taken from the undocumented example dg-demo.py).

I have seen the example in the tutorial Handling Domains with Different Materials but the use of MeshFunction seems clunky compared to defining the forcing term directly on a submesh where it is straightforward to access the coordinates directly.

For as much as I have searched for examples, I could not find anything. The closer that I found was this unanswered question. I think it would be very useful to include such an example somewhere. Or am I the only one dealing with a forcing term that is space dependent?

Thanks

Answer 1

You could include your forcing term in a subdomain integral:

mesh = Mesh("mymesh.xml")
subdomains = MeshFunction("size_t", mesh, "mymesh_physical_region.xml");
dxx = dx[subdomains]
V = FunctionSpace(mesh, 'Lagrange', 2)
u, v = TrialFunction(V),  TestFunction(V)
f = Expression("pow(x[0], 2)")
L = f*v*dxx(1)

Now it will only be integrated over subdomain 1.

Comments

Now that you mention it, something similar is mentioned in Section 1.5.3 of the FEniCS book, but your answer is clearer.

Related to this, suppose in the example of Section 1.5.1 of the FEniCS book that the function k(x,y) was defined as a space dependent k0(x,y) is subdomain $\Omega_1$ and k1(x,y) in $\Omega_2$ with $\Omega_1$, $\Omega_2$ some complex subdomains. Would one use again the approach you described?

---

Of course. You don't need to bother where is k0 and k1 defined when defining them. You just involve appropriate subdomain integrals:

L = k0*v*dxx(0) + k1*v*dxx(1)