← Back to team overview

dolfin team mailing list archive

Re: Non-intuitiveness with assembling over sub_domains

 

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. 

Sounds like this is best managed above the UFC/UFL layer.

> > I was going to suggest this but I realize now we would need to change
> > the UFC interface (and FFC). It is currently based on the function
> > 
> >  num_foo_integrals
> 
> BTW, is there a bug in the code generation for this function? If I do:
> element = FiniteElement("Lagrange", triangle, 1)
> f = Coefficient(element)
> M = inner(f, f)*dx + inner(f, f)*dx(2)
> 
> the return value of num_cell_integrals is 3, but create_cell_integral(),
> only return an integral in case 0 and 2.

No, this is correct. But I think it is suboptimal. When create_cell_integral 
is called using 1 as argument the zero pointer is returned. Dolfin then needs 
to check this. It would be more intuitive to let ufc check it, and to expand 
ufc with foo_integral_markers(uint * markers) or something. But it might as 
well not be worth the ufc update.

Johan


> > returning the number of integrals and
> > 
> >  create_foo_integral(i)
> > 
> > returning the integral for the terms on sub domain i.
> > 
> > This would need to be extended, possibly by overloading with
> > 
> >  create_foo_integral()
> > 
> > which would return the terms that are present on all sub domains.
> 
> Should it return an array of integrals over sub domains, or one integral
> like: Integrand_0*dx(0) + Integrand_1*dx(1) --> (Integrand_0 +
> Integrand_1)*dx(-1) ? The second option will result in a lot of redundant
> code being generated.
> 
> > Is it worth the effort of updating ufc.h, the UFC manual, the
> > assembler and FFC?
> 
> No.
> 
> Kristian
> 
> > --
> > Anders
> > 
> > -----BEGIN PGP SIGNATURE-----
> > Version: GnuPG v1.4.9 (GNU/Linux)
> > 
> > iEYEARECAAYFAkvDJpkACgkQTuwUCDsYZdG11gCfTB2E+6XObaar1f1yDqYt1vzj
> > dywAn3xh49BZFoIiiYd4GnDFDYTX1eHx
> > =N5PH
> > -----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



Follow ups

References