Index SubDomains with Array in dx

Hey,

I have a quick problem: Is there an elegant solution to index two or more subdomains at the same time in dx?

I basically want to write something like dx([0,1]), instead of dx(0)+dx(1)

Minimal (not) working examaple:

from __future__ import print_function
from fenics import *


mesh = UnitSquareMesh(64, 64)
tol  = 1e-14
# Create Subdomain
class Omega_0(SubDomain): # bottom subdomain
    def inside(self, x, on_boundary):
        return x[1] <= 0.25 + tol

# Create Subdomain
class Omega_1(SubDomain): # middle subdomain
    def inside(self, x, on_boundary):
        return x[1] <= 0.75 + tol and x[1] >= 0.25 - tol

# Create Subdomain
class Omega_2(SubDomain): # top subdomain
    def inside(self, x, on_boundary):
        return x[1] >= 0.75 - tol


materials   = MeshFunction("size_t", mesh, 2)
subdomain0  = Omega_0()
subdomain1  = Omega_1()
subdomain2  = Omega_2()

subdomain0.mark(materials, 0)
subdomain1.mark(materials, 1)
subdomain2.mark(materials, 2)

V = FunctionSpace(mesh, 'Lagrange',1)
u = TrialFunction(V)
v = TestFunction(V)

dx = Measure("dx", domain=V, subdomain_data = materials)
# this works
a = Constant(1.)*dot(grad(u), grad(v))*(dx(0) + dx(1)) + Constant(2.)*dot(grad(u), grad(v))*dx(2)
# this does not work, I find it however more elegant
a = Constant(1.)*dot(grad(u), grad(v))* dx([0,1])+ Constant(2.)*dot(grad(u), grad(v))*dx(2)
L = v*dx

bc = DirichletBC(V, Constant(1.), 'on_boundary')

u = Function(V)
solve(a==L, u, bc)

You can use tuples, i.e.

a = Constant(1.)*dot(grad(u), grad(v))* dx(tuple([0,1]))+ Constant(2.)*dot(grad(u), grad(v))*dx(2)
2 Likes

works like a charm, thank you :slight_smile:!