← Back to team overview

dolfin team mailing list archive

Re: Non-intuitiveness with assembling over sub_domains

 



On 15 April 2010 20:57, Anders Logg <logg@xxxxxxxxx> 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.

This makes sense and it looks nice on paper, but it the amount of code that has to be added to integrals that are, most likely, never used by the user could cause the g++ compiler to choke.
I don't see the problem in keeping the definition that we have now, after all, dx, ds, dS etc. are simply convenience tokens defined on subdomain 0.

Kristian

> 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

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.9 (GNU/Linux)

iEYEARECAAYFAkvHYZMACgkQTuwUCDsYZdHXbACcDJq3vs21IfiKTrLYHYvBO5Oq
oRkAn1uO39GGfhacshY9MqWQ5oA1eAjV
=iOwC
-----END PGP SIGNATURE-----

_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
Unsubscribe : https://launchpad.net/~dolfin
More help   : https://help.launchpad.net/ListHelp



Attachment: signature.asc
Description: OpenPGP digital signature


Follow ups

References