← Back to team overview

dolfin team mailing list archive

Re: [Bug 734527] Re: New style sub-domains do not carry across form transformation

 

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

--
Marie


Follow ups

References