dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

> > 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! :)

I guess...

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

That's also possible in C++. It's always been unclear what that does
so it should be removed.

--
Anders

Attachment: signature.asc
Description: Digital signature