dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday May 30 2011 13:33:29 Marie E. Rognes 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 guess was that:

  dxp = dx(cell_domains=primal_cell_domains)
  L = g*dxp(1) + f*dxp(2)

should work as previously. 

But:

  dxp = dx(cell_domains=primal_cell_domains)
  dxs = dx(cell_domains=second_cell_domains)
  L = g*dxp(1) + f*dxs(2)

should not.

Johan
> >> 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
> 
> --
> Marie
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp