# Mark Subdomains in 3D

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(

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(

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.

subdomain1 = Omega1()

subdomain1.

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(

subdomain_no = subdomains.

k.vector(

## 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(

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(

plot(mesh3d, "3D mesh");

V = FunctionSpace(

subdomains = MeshFunction(

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.

subdomain0 = Omega0()

subdomain0.

subdomain1 = Omega1()

subdomain1.

plot(subdomains)

V0 = FunctionSpace(

k = Function(V0)

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

help = numpy.asarray(

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/

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_

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!