dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

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?

> > > > 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.
> 
> It's a very easy check for DOLFIN and it avoids adding an extra
> function to UFC

Ok.

Johan

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