Defining a function on a subdomain

Question

Hi,

I have built an external cube mesh file using gmsh and dolfin-convert. There are two subdomains with labels 1 and 2 that make up the volume of the cube.

I have succesfully imported all relevant information using the following three lines.

mesh = Mesh("cube.xml")                                                  
subdomains = MeshFunction("size_t", mesh, "cube_physical_region.xml")    
boundaries = MeshFunction("size_t", mesh, "cube_facet_region.xml")

I am unsure of how to manipulate the subdomains (I plan on defining poisson equation source terms on them). E.g. Instead of

V = FunctionSpace(mesh,'Lagrange',1)
func = project(1.0,V)

I would like to define func1 and func2 on subdomains 1 and 2.

Thanks

Edit: The manual suggests the following lines

V = FunctionSpace(mesh,'Lagrange',1)
func = Function(V)

func_values = [0,1.0,2.0]
for cell_no in range(len(subdomains.array())):                          
    subdomain_no = subdomains.array()[cell_no]                          
    func.vector()[cell_no] = func_values[subdomain_no]                

but I get an error message saying the index cell_no is out of range.

Comments

I also had a similar problem to yours and here is the answer I got:

https://bitbucket.org/fenics-project/dolfin/issue/382/bug-when-using-expression-and-subdomains

Answer 1

Hi, I believe it is more elegant to define functions over the entire domain and include them in the source term by integrating over the desired subdomains. Consider

from dolfin import *

mesh = UnitSquareMesh(10, 10)

subdomains = CellFunction('size_t', mesh, 0)
for cell in cells(mesh):
    M = cell.midpoint()
    if M.y() > M.x():
        subdomains[cell] = 1

plot(subdomains)

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

f0 = Constant(0)  # source to be used in subdomain 0
f1 = Constant(1)  # source to be usde in subdomain 1

dx = Measure('dx')[subdomains]
a = inner(grad(u), grad(v))*dx('everywhere')
# Inluclude sources
L = inner(f0, v)*dx(0) + inner(f1, v)*dx(1)
bc = DirichletBC(V, Constant(0), DomainBoundary())
uh = Function(V)

solve(a == L, uh, bc)
plot(uh)
interactive()

Comments

Thanks for the reply. Would I be right in thinking this would also work if I wanted f0 and f1 to be permitivity functions? I.e. If I'm solving

$$ \nabla \cdot (f \nabla u) = -\rho$$

I could replace

a = inner(grad(u), grad(v))*dx('everywhere')

with

a = inner(f0*grad(u),grad(v))*dx(0) + inner(f1*grad(u),grad(v))*dx(1)

?

As an aside, here is how I get my original solution working

W = FunctionSpace(mesh,"Lagrange",1) # For solution function
V = FunctionSpace(mesh,"DG",0) # For source and permittivity
k = Function(V)

dm = V.dofmap()                                                        
helper = numpy.asarray(subdomains.array(),dtype=numpy.int32)             
for cell in cells(mesh):                                                 
    helper[dm.cell_dofs(cell.index())] = subdomains[cell]                

k_values=[0.0,3.0,4.0]                                                   
k.vector()[:] = numpy.choose(helper,k_values)  

---

Hi, the bilinear form with dx(0) and dx(1) will work. +1 for using dofmap in your original solution.