dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Jun 01, 2011 at 12:07:58AM +0200, Marie E. Rognes wrote:
> On 06/01/2011 12:05 AM, Anders Logg wrote:
> >On Wed, Jun 01, 2011 at 12:01:50AM +0200, Marie E. Rognes wrote:
> >>On 05/30/2011 11:51 PM, Johan Hake wrote:
> >>>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 :)
> >>>
> >>>+
> >>>
> >>
> >>Ok, so we are three in favor. Anders, do you want to veto?
> >
> >I think I was the first one to say it looks good.
> >
>
> Ok, great! I must have misunderstood your "I don't want to change
> the current syntax again" ;-)

I wrote that before you suggested your solution. Then I wrote

  "I like a.ds = exterior_facet_domains!"

in response to that suggestion. :-)

--
Anders