Mark Subdomains in 3D

Asked by Michael Trogdon on 2013-02-07

Ok so I know there are a lot of posts on this but I can't seem to find one to address the error I am getting. I am simply trying to make two subdomains in a 3D mesh but I keep getting the following error:

Traceback (most recent call last):
  File "BulkDiffusion.py", line 45, in <module>
    k.vector()[cell_no] = k_values[subdomain_no]
IndexError: list index out of range

Here is a snippet of the relevant code, if someone can notice something I am doing wrong it would be greatly appreciated. Let me know if you have any questions. (The .xml file is just a 3D cylinder with spherical caps).

#Create mesh from Coli problem
mesh = Mesh('coli.xml')
V = FunctionSpace(mesh, "Lagrange", 1)

subdomains = MeshFunction("uint", mesh, 3)

class Omega0(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] <= 1 else False
class Omega1(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] >= 1 else False

# Mark subdomains with numbers 0 and 1
subdomain0 = Omega0()
subdomain0.mark(subdomains, 0)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 1)

V0 = FunctionSpace(mesh, 'DG', 0)
k = Function(V0)

k_values = [1.5, 50] # values of k in the two subdomains
for cell_no in range(len(subdomains.array())):
 subdomain_no = subdomains.array()[cell_no]
 k.vector()[cell_no] = k_values[subdomain_no]

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Michael Trogdon
Solved:
2013-02-08
Last query:
2013-02-08
Last reply:
2013-02-08
Jan Blechta (blechta) said : #1

Consider that not all cells belongs to Omega0 or Omega1. I'm not sure what's the exact algorithm but maybe if whole cell is not completely x[1]<=1 or x[1]>=1 then it is not marked by any domain.

Since MeshFunction is not initialized - it has random values (I think it changed since dolfin 1.0.0 but not sure), for these unmarked cells subdomain_no gets theses random values and k_values has no corresponding key.

I reproduced this behaviour with:

mesh = UnitSquare(3, 3)
subdomains = MeshFunction("uint", mesh, 2)

Consider also that DG0 dof numbering may not correspond to CellFunction numbering in newer versions of DOLFIN. There are such renumbering issues but I'm not sure wheter it occurs for DG0.

Jan

M Waste (michael-waste) said : #2

hi,

strange enough. It was not possible to create two subdomains with different k values on a cylinder with caps.
Generate always a crash.
Also with a cgs cylinder I created only messy orderd k_values.

But it is possible to create different k_values in domains with a BoxMesh.

Maybe a hint for you to do this with your mesh.
Here the code with BoxMesh.

from dolfin import *
import numpy # needed for k_values

mesh3d = BoxMesh(0,0,0,2,3,3,10,10,10)
plot(mesh3d, "3D mesh");

V = FunctionSpace(mesh3d, "Lagrange", 1)

subdomains = MeshFunction("size_t", mesh3d, 3)

class Omega0(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] <= 1 else False
class Omega1(SubDomain):
 def inside(self, x, on_boundary):
  return True if x[1] >= 1 else False

# Mark subdomains with numbers 0 and 1
subdomains.set_all(0) # to avoid random numbering
subdomain0 = Omega0()
subdomain0.mark(subdomains, 0)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 1)
plot(subdomains)
V0 = FunctionSpace(mesh3d, 'DG', 0)
k = Function(V0)

k_values = [1.5, 50.] # values of k in the two subdomains

help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
k.vector()[:] = numpy.choose(help, k_values)

plot(project(k))

interactive()

good luck

Michael

Jan Blechta (blechta) said : #3

Well, it's not strange at all. I checked dolfin/mesh/SubDomain.cpp to
realize how SubDomain::mark() works. For entity to be marked, it is
needed that all incident vertices of that entity plus midpoint when
passed to SubDomain's inside method to yield true.

So you have to construct appropriate SubDomains. For example: Omega0 =
{x<=1.2}, Omega1 = {x>=0.8} with difference in values decreasing with
cell size... Or you can simply initialize values of subdomains
MeshFunction before marked to some default value.

You can also think of the situation as flaw of SubDomain::mark() method
and report a bug or improve it. It could take an additional voluntary
argument bool midpoint_only = false, which would ensure that only
midpoint is checked if midpoint_only==true. And/or another voluntary
bool argument could be added, that would marked entity if at least on
vertex of (incident vertices + midpoint) was inside...

Jan

Jan Blechta (blechta) said : #4

Already reported as bug 1119803

Michael Trogdon (mindgame252) said : #5

Thank you both for the reply. The idea of initializing the subdomains to zero to avoid random numbering seems to be working correctly. Appreciate the help!