dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday May 30 2011 14:47:32 Martin Sandve Alnæs wrote:
> 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.

Yes, I assumed that.

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

I suggest we restrict this in UFL. Maybe as part of Form.__add__ or something. 
It checks that the same meassures are used for all integrals of the same 
domain.

Johan