Mark Subdomains in 3D

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

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

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

Already reported as bug 1119803

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!