dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Thursday April 15 2010 12:36:36 Anders Logg wrote:
> On Thu, Apr 15, 2010 at 12:29:33PM -0700, Johan Hake wrote:
> > On Thursday April 15 2010 12:25:14 Anders Logg wrote:
> > > On Thu, Apr 15, 2010 at 12:20:40PM -0700, Johan Hake wrote:
> > > > On Thursday April 15 2010 11:57:23 Anders Logg wrote:
> > > > > On Thu, Apr 15, 2010 at 11:51:00AM -0700, Johan Hake wrote:
> > > > > > 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?
> > > > > 
> > > > > The point is that UFC should not create the -1 integral. It will be
> > > > > replaced by other integrals. The terms defined on dx(-1) will be
> > > > > added to all the other terms.
> > > > > 
> > > > > Something I realized now is what happens if there are no other
> > > > > terms, but then we must simply add it to dx(0).
> > > > > 
> > > > > And we should probably also add it to all integrals numbered before
> > > > > the highest numbered integral, even if those are missing, so if we
> > > > > have
> > > > > 
> > > > >   f*dx + g*dx(10)
> > > > > 
> > > > > then then f*dx term will be added to the dx(0), dx(1), ..., dx(9)
> > > > > integrals.
> > > > 
> > > > I think this gets complicated. Wouldn't the easiest sollution be to
> > > > let dx(0) always denote the whole domain? This would then be handed
> > > > naturally in DOLFIN, by always integrating over domain 0 if such an
> > > > integral exist. Passing a MeshFunction with domain markers of 0
> > > > could then either be ignored or we could issue a warning?
> > > 
> > > You mean like we have today?
> > 
> > No.
> > 
> > > Well, that's the easiest option. :-)
> > > (And then we don't need to remove the SubDomain possibility.)
> > > 
> > > It all boils down to this, should dx mean integral over the whole
> > > domain or integral over all cells marked by 0?
> > 
> > Eactly, and I suggest we let the 0th boundary be reserved to always
> > represent the _whole_ boundary, which isn't that unintuitively?
> 
> But what would it mean then for a cell to be marked 0?
> 
> Or should it not be allowed to mark a cell as 0?

Exactly. We might just ignore it or issue a warning?

> Perhaps this would work (or maybe that is what you suggest):
> 
> 1. dx = dx(0) = whole domain
> 
> 2. The assembler iterates over the domain and always adds the
> contributions from integral number 0 (for each cell). It then check if
> the current cell is marked something other than 0 (like 3) and in that
> case adds also the contribution from integral number 3.

Yes.

> This would lead to complications in the assembler: needing to call
> tabulate_tensor twice and add the results.

True. Not sure which one of the proposed alternative gives the smallest trade 
off. 

Johan

> --
> Anders