dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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

--
Anders

Attachment: signature.asc
Description: Digital signature