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.
