Creating many subdomains iteratively

Asked by Alex Evans on 2013-02-14

Hi,

I am wanting to solve a problem very similar to that presented in http://fenicsproject.org/documentation/tutorial/materials.html
however I would like to be able to create many subdomains (e.g. on the order of 256) without having to manually define each one as described in that demo. Is there a way to do this?

ie. I would like to have the effects of this:

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

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

but for many such subdomains.

Thanks in advance for your time and for your patience over what is undoubtedly a stupid question :).

Alex

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Alex Evans
Solved:
2013-02-15
Last query:
2013-02-15
Last reply:
2013-02-14
Garth Wells (garth-wells) said : #1

On 14 February 2013 10:51, Alex Evans
<email address hidden> wrote:
> New question #221862 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/221862
>
> Hi,
>
> I am wanting to solve a problem very similar to that presented in http://fenicsproject.org/documentation/tutorial/materials.html
> however I would like to be able to create many subdomains (e.g. on the order of 256) without having to manually define each one as described in that demo. Is there a way to do this?
>

You could probably use some Python trickery to generate the classes.

I would like to change SubDomain::inside to return an integer which
would make it possible to create multiple subdomains with just one
loop over the mesh.

Garth

> ie. I would like to have the effects of this:
>
> class Omega0(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] <= 0.5 else False
>
> class Omega1(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] >= 0.5 else False
>
> but for many such subdomains.
>
> Thanks in advance for your time and for your patience over what is undoubtedly a stupid question :).
>
> Alex
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

This is really a python question, keep in mind that it's just a regular
class:

class Omega0(SubDomain):
    def __init__(self, x_max):
        self.x_max = x_max
    def inside(self, x, on_boundary):
        return x[1] <= self.x_max

Otherwise the answer depends on in what way the subdomains differ. Just
different values? Use variations of the above. Different functions? Then
there's no good way around defining different functions (i.e. classes).

Martin

On 14 February 2013 12:06, Garth Wells <<email address hidden>
> wrote:

> Question #221862 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/221862
>
> Status: Open => Answered
>
> Garth Wells proposed the following answer:
> On 14 February 2013 10:51, Alex Evans
> <email address hidden> wrote:
> > New question #221862 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/221862
> >
> > Hi,
> >
> > I am wanting to solve a problem very similar to that presented in
> http://fenicsproject.org/documentation/tutorial/materials.html
> > however I would like to be able to create many subdomains (e.g. on the
> order of 256) without having to manually define each one as described in
> that demo. Is there a way to do this?
> >
>
> You could probably use some Python trickery to generate the classes.
>
> I would like to change SubDomain::inside to return an integer which
> would make it possible to create multiple subdomains with just one
> loop over the mesh.
>
> Garth
>
> > ie. I would like to have the effects of this:
> >
> > class Omega0(SubDomain):
> > def inside(self, x, on_boundary):
> > return True if x[1] <= 0.5 else False
> >
> > class Omega1(SubDomain):
> > def inside(self, x, on_boundary):
> > return True if x[1] >= 0.5 else False
> >
> > but for many such subdomains.
> >
> > Thanks in advance for your time and for your patience over what is
> undoubtedly a stupid question :).
> >
> > Alex
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Anders Logg (logg) said : #3

On Thu, Feb 14, 2013 at 12:01:22PM -0000, Martin Sandve Alnæs wrote:
> Question #221862 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/221862
> Martin Sandve Alnæs proposed the following answer:
> This is really a python question, keep in mind that it's just a regular
> class:
> class Omega0(SubDomain):
> def __init__(self, x_max):
> self.x_max = x_max
> def inside(self, x, on_boundary):
> return x[1] <= self.x_max
> Otherwise the answer depends on in what way the subdomains differ. Just
> different values? Use variations of the above. Different functions? Then
> there's no good way around defining different functions (i.e. classes).
> Martin

Since you have so many subdomains, I assume you can express the
numbering somehow based on the element numbers or coordinates. In that
case, just create a CellFunction on the mesh, then iterate over the
cells and set the markers based on the coordinates or indices.

markers = CellFunction("size_t", mesh)
markers.set_all(0)

for cell in cells(mesh):
    i = cell.index()
    x = cell.midpoint()
    markers.set(i, f(i, x))

(untested code)

--
Anders

Alex Evans (evansalex) said : #4

Thanks for the suggestions all, this may not be relevant but I think the elephant in the room here is that (with my imported mesh) the cell indexing isn't the same for the mesh function and the function of my vector space that I am trying to define:

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

for cell_no in range(len(subdomains.array())):
            subdomain_no = subdomains.array()[cell_no]
            kxx.vector()[cell_no] = noisexx[subdomain_no]

the 'kxx.vector()' is seemingly ordered in a different fashion to that of the cell function array.

A representing of what I have tried to do is below, due to this issue of different orderings the cells aren't given the correct value that should correspond to their position in space.

________________________________

from dolfin import *
from random import gauss

mesh = Mesh("meshfiniteheight10.xml")
mesh.order()

V0 = FunctionSpace(mesh, 'DG', 0)
kxx = Function(V0)
kxy = Function(V0)
kyy = Function(V0)

noisexx = range(40)
noisexy = range(40)
noiseyy = range(40)

subdomains = MeshFunction('size_t', mesh, 2)
subdomains.set_all(0)

for i in range(20):
    for j in range(2):
        class Omega(SubDomain):
            def inside(self, x, on_boundary):
                return True if x[0] >= 0.5*i-5.0 and x[0] <= 0.5*i-4.5 and x[1] <= 0.5*j and x[1] >= 0.5*j-0.5 else False

# Mark subdomain cells with corresponding number
        subdomain = Omega()
        subdomain.mark(subdomains, 2*i+j)

for I in range(40):
    noisexx[I] = gauss(0,1)
    noisexy[I] = gauss(0,1)
    noiseyy[I] = gauss(0,1)

for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    kxx.vector()[cell_no] = noisexx[subdomain_no]
    kxy.vector()[cell_no] = noisexx[subdomain_no]
    kyy.vector()[cell_no] = noisexx[subdomain_no]
___________________________________________

Cheers
Alex

Jan Blechta (blechta) said : #5

Don't rely on arrays indexing.

Use Cell object to index CellFunction:
   cell_func[cell]

Use DofMap.cell_dofs(cell_index) to get dof number to coresponding cell_index.

Alex Evans (evansalex) said : #6

Jan, thanks for your comment but as I'm unfamiliar with the DofMap I'm not sure I follow - are these changes to be made to all of the code I have or is it a better way of doing:

for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    kxx.vector()[cell_no] = noisexx[subdomain_no]
    kxy.vector()[cell_no] = noisexx[subdomain_no]
    kyy.vector()[cell_no] = noisexx[subdomain_no]

if so could you clarify?

Thanks
Alex

Johan Hake (johan-hake) said : #7

On 02/14/2013 05:45 PM, Alex Evans wrote:
> Question #221862 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/221862
>
> Alex Evans posted a new comment:
> Jan, thanks for your comment but as I'm unfamiliar with the DofMap I'm
> not sure I follow - are these changes to be made to all of the code I
> have or is it a better way of doing:

Try:

for cell in cells(mesh):
  subdomain_no = subdomains[cell] # Not array index
  dof_ind = V0.dofmap().cell_dofs(cell.index())[0]
  kxx.vector()[dof_ind] = noisexx[subdomain_no]

I have not tested the code, and I have no clue what noisexx is.

Johan

>
> for cell_no in range(len(subdomains.array())):
> subdomain_no = subdomains.array()[cell_no]
> kxx.vector()[cell_no] = noisexx[subdomain_no]
> kxy.vector()[cell_no] = noisexx[subdomain_no]
> kyy.vector()[cell_no] = noisexx[subdomain_no]
>
> if so could you clarify?
>
> Thanks
> Alex
>

Alex Evans (evansalex) said : #8

Thanks Johan/Jan, that's cracked it.