← Back to team overview

dolfin team mailing list archive

Re: Non-intuitiveness with assembling over sub_domains

 

Anders Logg wrote:
On Fri, Apr 16, 2010 at 10:39:05AM +0200, Kristian Oelgaard wrote:
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

That's the easiest option... What does Marie think? (She wanted to
change the behavior in the first place.) I can live with both options.

So, conclusion is to leave as before? If so, Marie will just have to be creative.

--
Marie



Follow ups

References