dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

> If not, maybe Martin can make the necessary changes in UFL, and I
> can do the same for DOLFIN?

Sounds good.

--
Anders