dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

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)

?


> 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.

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

--
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.
>