dolfin team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

On 30 May 2011 22:33, Marie E. Rognes <meg@xxxxxxxxx> wrote:
> On 05/30/2011 07:36 PM, Anders Logg wrote:
>>
>> On Mon, May 30, 2011 at 10:24:48AM -0700, Johan Hake wrote:
>>>
>>> On Monday May 30 2011 04:33:29 Martin Sandve Alnæs wrote:
>>>>
>>>> On 30 May 2011 12:51, Marie E. Rognes<meg@xxxxxxxxx>  wrote:
>>>>>
>>>>> On 05/30/2011 11:26 AM, Martin Sandve Alnæs wrote:
>>>>>>
>>>>>> Since this feature implementation relies on modifying immutable
>>>>>> objects,
>>>>>> I'm not the least surprised you're getting problems. The bug is not
>>>>>> that
>>>>>> dolfin subdomains are not passed with forms, but that they are allowed
>>>>>> to be attached in the first place on an existing and assumed immutable
>>>>>> form object.
>>>>>
>>>>> Yes...
>>>>>
>>>>>> The short term solution to this bug is to revert pydolfin back to
>>>>>> providing subdomains as arguments to assemble and variationalproblem
>>>>>> where they belong, instead of attaching them to forms. I think this
>>>>>> should be done for fenics 1.0 if this bug is a problem.
>>>>>>
>>>>>> Improvements to the language for expressing subdomains of various
>>>>>> kinds
>>>>>> is in the design stage, but that won't happen before the summer.
>>>>>
>>>>> Specifications of subdomain does not belong as arguments to assemble
>>>>> and
>>>>> variational problem. If you have a form
>>>>>
>>>>>        L = g*v*dG
>>>>>
>>>>> where G is a part of a boundary (In semi-math, semi-UFL notation), G
>>>>> should be related to the form. Not to the matrix resulting from the
>>>>> assembly of the form.
>>>>>
>>>>> (cc to DOLFIN since the below involves DOLFIN mainly)
>>>>>
>>>>> The interface to VariationalProblem
>>>>>
>>>>>        VariationalProblem(., ., bcs, exterior_facet_domains,
>>>>> interior_facet_domains, cell_facet_domains)
>>>>>
>>>>> was rather suboptimal because it assumed implicitly that the bilinear
>>>>> and
>>>>> the linear form were defined over the same subdomains. That in,
>>>>> combination with dx = dx(0) etc, is highly bugprone.
>>>>>
>>>>> I care of course because if you want to use the same patent for an
>>>>> variational problem with automatic adaptivity, and take care of the
>>>>> above, the required  input will look something like this
>>>>>
>>>>>        VariationalProblem(., ., bcs,
>>>>>                        primal_bilinear_exterior_facet_domains,
>>>>>                        primal_bilinear_interior_facet_domains,
>>>>>                        primal_bilinear_cell_domains,
>>>>>                        primal_linear_exterior_facet_domains,
>>>>>                        primal_linear_interior_facet_domains,
>>>>>                        primal_linear_cell_domains,
>>>>>                        goal_exterior_facet_domains,
>>>>>                        goal_interior_facet_domains,
>>>>>                        goal_cell_domains)
>>>>>
>>>>> which I can't live with.
>>>>>
>>>>> The Coefficient/Function magic must involve some of the same issues as
>>>>> this. I imagine that a similar way of fixing it should be possible.
>>>>
>>>> I'm not saying it can't be implemented, only that it can't
>>>> be implemented well within the FEniCS 1.0 timeframe, and
>>>> that fixing the current solution will lead down a bad path.
>>>> I'm hoping for a much better solution later this year, but
>>>> that will require some design work first.
>>>>
>>>> An alternative short term approach could be to attach the data to the
>>>> measures.
>>>>
>>>> dxp = dx(cell_domains=primal_cell_domains)
>>>> dsp = ds(exterior_facet_domains=primal_exterior_facet_domains)
>>>> L = g*dxp(1) + f*dsp(0)
>>>>
>>>> and making sure that this data follows measure objects around.
>>>> They can then be collected in ufl preprocess just like functions
>>>> and function spaces. Then the connection between the
>>>> meshfunction and the dx(i) index looks more explicit as well.
>>>
>>> I like that syntax. I guess we only allow one type of cell integral
>>> within one
>>> form?
>
> I don't see how this last guess would follow from the above. Why?
>
>
>>>
>>> My conserns regards changing the subdomains after creating the form, as
>>> this
>>> can be a convinient way of reusing a compiled form. But I guess letting
>>> domain
>>> arguments in assemble override domains in the form should fix that.
>>
>> I like the current syntax and would not like to change it again. An
>> important point is that it allows the same syntax to be used in both
>> Python and C++:
>>
>> a.exterior_facet_domains = exterior_facet_domains
>> a.exterior_facet_domains = exterior_facet_domains;
>>
>> This won't work with dx etc.
>>
>
> The same syntax is not used in C++ and Python for coefficients. With that
> argument, we should not allow
>
>        f = Function(V)
>        a = f*ds
>
> but force
>
>        a.f = f
>
> in Python also. (Which I of course do not advocate, but also makes me not
> buy the argument.)
>
>>
>>
>> I suggest we keep the current syntax and find an improved
>> implementation later.
>>
>
>
> I'll throw in the following suggestion in the mix, which is based on same
> idea as Martin's, but a bit terser:
>
> -------------------------------------------------------------------
> Defining sub-domains by mesh functions using boundary domains as an
> example
> -------------------------------------------------------------------
>
> Python
> ********************************************************************
>
> boundaries = FacetFunction("uint", mesh)
> ds = ds(boundaries)
> a = u*v*ds(0) + f*v*ds(1)
>
> C++
> ********************************************************************
>
> UFL:
> --------------------------------------------------------------------
> a = u*v*ds(0) + f*v*ds(1)
>
> main:
> --------------------------------------------------------------------
>
> Form a(V);
> FacetFunction<uint> boundaries(mesh);
> a.ds = boundaries;
>
> Advantages:
>
> * Backwards compatible (except last 2 releases or so)
> * Very similar construction to Functions/Coefficients
> * ds is way shorter than exterior_facet_domains
> * Retains desired link between form and domains of integration
> * Does not mess with UFL immutability
> * Can implement with ufl.Measure.metadata or revamp Measure further

I think this looks great :)

Johan: preprocess can fetch the metadata from measures,
shouldn't be too hard. Much easier than the coefficients
and function spaces anyway, since the integrals and
thus the measures are available directly as a list in the form.

The only issue I see is that:
- L = g*dxp(1) + f*dxs(1) will have 1 meaning two different subdomains
- ufc uses one unique numbering of all integrals
so to support this we need to renumber subdomains between ufl and ufc code.
But that would be additional functionality we do not have today, right?
So we don't have to handle that now.

Martin