solving pde on disconnected meshes

Question

Hi!
I would like to solve a diffusion pde over two disconnected domains that communicate through two surfaces (e.g. two cubes with a flux applied on one side of each cube).

I created a geometry using gmsh, for which i marked the cells as 29 and 229, and two faces as 27 and 227 (available here: http://pkh.as.uky.edu/html/files/fenics.tgz )

I am attempting to define a weak form that is defined on domains dx(29) and dx(229) with fluxes applied to ds(27) and ds(227) [similar to http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/subdomains-poisson/python/documentation.html] , but I get a segmentation fault on the first solve. If I instead solve the same problem over dx() and ds(), the solver is fine, so I suspect I am missing something stupid.

Question: is there a better way to define the weak form/solver so that I can solve the PDE over two disconnected domains?

My code is below:

from dolfin import *

name ="twocube"
subdomainMarkers=[29,229] # from 'Physical Volume(29)' in .geo file
boundaryMarkers = [27,227] # from 'Physical Surface(28)' in .geo file

mesh = Mesh(name+".xml")

## load subdomains, boundaries
subdomains = MeshFunction("size_t", mesh, name+"_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, name+"_facet_region.xml")

dx = Measure("dx",domain=mesh)
ds = Measure("ds",domain=mesh,subdomain_data=boundaries)

# Surfaces
for i,marker in enumerate(boundaryMarkers):
  print "area marker ", marker, " ", assemble(Constant(1.)*ds(marker))

#print "vol", assemble(Constant(1.)*dx)
for i,marker in enumerate(subdomainMarkers):
  submesh = SubMesh(mesh, subdomains, marker)
  print "vol marker ", marker, " ", assemble(Constant(1.)*dx(submesh))

# Apply these markings to our surface/distance measures
ds = ds(domain=mesh,subdomain_data=boundaries)
dx = dx(domain=mesh)
V = FunctionSpace(mesh,"CG",1)

v = TestFunction(V)
u_n = Function(V)
u_0 = Function(V)

Du = Constant(1)
dt = 1.
#
bcs=[]
#bcs = [
# DirichletBC(V,Constant(2.),boundaries,boundaryMarkers[0]),
# DirichletBC(V,Constant(0.),boundaries,boundaryMarkers[1])]

# LHS
# works
#F = ((u_n-u_0)*v/dt + Du*inner(grad(u_n), grad(v)))*dx()
#F+= -Constant(0.1)*v*ds()
# doesn't work
F = ((u_n-u_0)*v/dt + Du*inner(grad(u_n), grad(v)))*dx( subdomainMarkers[0] )
F+= ((u_n-u_0)*v/dt + Du*inner(grad(u_n), grad(v)))*dx( subdomainMarkers[1] )
F+= -Constant(0.1)*v*ds( boundaryMarkers[0])
F+= Constant(0.1)*v*ds( boundaryMarkers[1])

u_0.interpolate(Constant(2.))
u_n.vector()[:] = u_0.vector()[:]

t=0
tstop=5
while t < tstop+DOLFIN_EPS:
    solve(F==0, u_n,bcs)
    print "############### ", t

    submesh = SubMesh(mesh, subdomains, subdomainMarkers[0])
    #print "vol", assemble(Constant(1.)*dx)
    for i,marker in enumerate(subdomainMarkers):
      submesh = SubMesh(mesh, subdomains, marker)
      volSubMesh = assemble(Constant(1.)*dx(submesh))
      concSubMesh = assemble(u_n*dx(submesh))/volSubMesh
      print "conc marker ", marker, " ", concSubMesh

    u_0.assign(u_n)
    t += float(dt)

huskeypm · Answer 1 · Dec 10, 2015

From a Fenics enthusiast:

The problem generating the errors in your script was that you tried to
extract a submesh. Since you already have MeshFunction with
appropriate markers, you may use them to compute measures on each
subdomain (each box). Unless for some reasons you really want/need to
extract submeshes but usually it is not a good idea.

from dolfin import *

name ="twocube"
subdomainMarkers=[29,229]  # from 'Physical Volume(29)' in .geo file  
boundaryMarkers = [27,227] # from 'Physical Surface(28)' in .geo file 

mesh = Mesh(name+".xml")

## load subdomains, boundaries 
subdomains = MeshFunction("size_t", mesh, name+"_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, name+"_facet_region.xml")

# COMMENT: Provide mesh function to the spatial measure
dx = Measure("dx",domain=mesh, subdomain_data=subdomains)
ds = Measure("ds",domain=mesh,subdomain_data=boundaries)

# Surfaces 
for i,marker in enumerate(boundaryMarkers):
  print "area marker ", marker, " ", assemble(Constant(1.)*ds(marker))

#print "vol", assemble(Constant(1.)*dx)
for i,marker in enumerate(subdomainMarkers):
    # COMMENT: No need to extract a submesh. Use the markers.
#  submesh = SubMesh(mesh, subdomains, marker)
  print "vol marker ", marker, " ", assemble(Constant(1.)*dx(marker))


V = FunctionSpace(mesh,"CG",1)
v = TestFunction(V)
u_n = Function(V)
u_0 = Function(V)

Du = Constant(1)
Dub = Constant(0.1)
dt = 1.
#
bcs=[]
# LHS 
F = ((u_n-u_0)*v/dt + Du*inner(grad(u_n), grad(v)))*dx( subdomainMarkers[0] )
F += ((u_n-u_0)*v/dt + Dub*inner(grad(u_n), grad(v)))*dx( subdomainMarkers[1] )
F -= Constant(0.1)*v*ds( boundaryMarkers[0])
F += Constant(0.1)*v*ds( boundaryMarkers[1])
#
u_0.interpolate(Constant(2.))
u_n.vector()[:] = u_0.vector()[:]

t=0
tstop=5
while t < tstop+DOLFIN_EPS:
    solve(F==0, u_n,bcs)
    print "############### ", t

    submesh = SubMesh(mesh, subdomains, subdomainMarkers[0])
    #print "vol", assemble(Constant(1.)*dx)
    for i,marker in enumerate(subdomainMarkers):
        # COMMENT: Again use markers instead of extracting submeshes.
#      submesh = SubMesh(mesh, subdomains, marker)
      volSubMesh = assemble(Constant(1.)*dx(marker))
      concSubMesh = assemble(u_n*dx(marker))/volSubMesh
      print "conc marker ", marker, " ", concSubMesh

    u_0.assign(u_n)
    t += float(dt)

solving pde on disconnected meshes

1 Answer