dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Thursday April 15 2010 11:39:28 Anders Logg wrote:
> On Mon, Apr 12, 2010 at 10:29:18AM -0700, Johan Hake wrote:
> > On Monday April 12 2010 08:43:54 Anders Logg wrote:
> > > On Mon, Apr 12, 2010 at 08:32:41AM -0700, Johan Hake wrote:
> > > > On Monday April 12 2010 08:14:09 Kristian Oelgaard wrote:
> > > > > On 12 April 2010 15:56, Anders Logg <logg@xxxxxxxxx> wrote:
> > > > > > On Thu, Apr 08, 2010 at 08:44:21AM -0700, Johan Hake wrote:
> > > > > >> 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.
> > > > > > 
> > > > > > Yes, this is the case.
> > > > > > 
> > > > > > I think the current behavior is correct, if one knows that ds =
> > > > > > ds(0).
> > > > > > 
> > > > > > But I understand it can be confusing. Marie and I discussed this
> > > > > > last week and I suggested that we let
> > > > > > 
> > > > > >  ds = ds(-1)
> > > > > > 
> > > > > > and that should denote integration over all sub domains, and
> > > > > > ds(0) would mean integration specifically over sub domain 0.
> > > > 
> > > > What would the default subdomain then be which we integrate over when
> > > > a SubDomain is passed?
> > > > 
> > > > I don't know if I really get what you are trying to accomplish.
> > > > 
> > > > Whould:
> > > >   f*ds + g*ds(0) + h*ds(1)
> > > > 
> > > > create three integrals one for "subdomain" -1 which returns the
> > > > integral over all subdomains, and one for each of 0 and 1.
> > > 
> > > Yes.
> > > 
> > > > Sounds like this is best managed above the UFC/UFL layer.
> > > 
> > > Yes, but I don't see how that should be done. Somehow we need to pass
> > > the information from
> > > 
> > >   f*ds + g*ds(0) + h*ds(1)
> > > 
> > > to the assembler.
> > > 
> > > One option would be to add f*ds to all the sub integrals. So we don't
> > > change the UFC interface (or DOLFIN) but do the following:
> > > 
> > > 1. ds = ds(-1)
> > > 
> > > 2. FFC adds ds(-1) to all its sub domain integrals.
> > 
> > But how would one integrate over just the ds domain?
> 
> Sorry for my late reply on this one...
> 
> Anyway, integrating over just the ds domain would just be
> 
>   f*ds
> 
> I think this should work out fine. Here's a summary (now with dx
> instead of ds).
> 
> 1. Omega = Omega_0 U Omega_1 U Omega_2 ...
> 2. dx means integral over Omega
> 3. dx(i) means integral over Omega_i
> 
> So if one writes
> 
>   f*dx + g*dx(0)
> 
> that means that the integral is a sum of two integrals, one over the
> whole domain and one just over the domain marked with 0.
> 
> This can alternatively be expressed as
> 
>   (f+g)*dx(0) + f*dx(not marked 0)

Ok.

> This can be handled without changing the UFC interface by making just
> one change in UFL (or FFC), which is to add the terms involving dx
> (which we can represent as dx(-1)) to all the other integrals.

It will be difficult to ask an UFC form to create the -1 integral as the 
interface now use unsigned int. But another int might be as god?

> Then DOLFIN will happily iterate over all cells and on each cell check
> how the cell is marked (=i) and evaluate integral number i.

Why does this sounds intuitive now and not when it was first introduced? Time 
has its effects! :)

But then we might remove the possibility (in PyDOLFIN) to pass a SubDomain as 
foo_domains, because it will not be clear what subdomain (marker) it should 
represent.

Johan

> --
> Anders