dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Thursday April 8 2010 08:36:24 Marie Rognes wrote:
> Johan Hake wrote:
> > On Thursday April 8 2010 04:51:24 Marie Rognes wrote:
> >> There is something suboptimal with regard to assembly over sub_domains
> >> specified
> >> by domains (in contrast to markers)
> >> 
> >> Say I have some functional
> >> 
> >> (*)    M = f*ds
> >> 
> >> where f is some function. Let 'domain' be some sub-domain.
> >> 
> >> Now, it is not clear to me what
> >> 
> >>     v = assemble(M, mesh=mesh, exterior_facet_domains=domain)
> >> 
> >> gives. It would be really convenient (for my purposes) if it gave
> >> the same as
> >> 
> >>     v = assemble(M, mesh=mesh)
> > 
> > Considering that ds defaults to ds(0) I think it is a logical behaviour.
> 
> So
> 
>     ds == ds(0)
> 
> ?

In ufl/objects.py 

  ds = Measure(Measure.EXTERIOR_FACET, 0)

So I guess yes.

> > When
> > you want the whole boundary you do not pass any exterior domain. If you
> > want the ds form integrated over a certain boundary you pass that
> > boundary.
> 
> Right, and I guess that most people actually know which boundary that
> they want to integrate over ;)
> 
> However, I just have some abstract form M and some subdomains that I try
> to assemble correctly.

But wouldn't you expect that a passed exterior domain would define the 
integration domain for your exterior form? Maybe I have to think more 
abstract...

> Is there any way of checking whether a form is *ds or *ds(0)? If indeed
> ds == ds(0), I guess not... Hm.

I guess not...

Johan


> --
> Marie
> 
> >> It would also make some sense, considering that no special
> >> tag has been defined in (*), and so the exterior_facet_domains
> >> argument could be ignored.
> >> 
> >> It does not at the moment though. At the moment, it seems
> >> as though ds -> ds(0) and the sub-domain is converted
> >> to markers, marking anything inside 'domain' as 0, and
> >> everything else as 1...?  This is really not obvious to me.
> > 
> > I do not know if changing the default representation of ds, or
> > documenting (better) that ds defaults to ds(0) could help. The latter
> > whould be the easiest fix.
> > 
> >> Could "we" change the interface to something slightly
> >> more transparent?
> >> 
> >> --
> >> Marie
> >> 
> >> 
> >> PS: Example attached.
> >> 
> >> 
> >> 
> >> from dolfin import *
> >> 
> >> class Domain(SubDomain):
> >>     def inside(self, x, on_boundary):
> >>         return x[0] < 0.5 and on_boundary
> >> 
> >> mesh = UnitSquare(8, 8)
> >> f = Expression("sin(pi*x[0])")
> >> M = f*ds
> >> 
> >> domain = Domain()
> >> 
> >> markers = MeshFunction("uint", mesh, 1)
> >> markers.set_all(0)
> >> domain.mark(markers, 1)
> >> 
> >> v1 = assemble(M, mesh=mesh)
> >> v2 = assemble(M, mesh=mesh, exterior_facet_domains=domain)
> >> v3 = assemble(M, mesh=mesh, exterior_facet_domains=markers)
> >> 
> >> print "It is not intuitive what these become ..."
> >> print "v1 = ", v1
> >> print "v2 = ", v2
> >> print "v3 = ", v3
> > 
> > At least v2 and v3 add up to v1 :)
> > 
> > If you changed
> > 
> >   markers.set_all(0)
> >   domain.mark(markers, 1)
> > 
> > to
> > 
> >   markers.set_all(1)
> >   domain.mark(markers, 0)
> > 
> > You would get the same.