dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2008/10/22 Anders Logg <logg@xxxxxxxxx>:
> On Wed, Oct 22, 2008 at 11:55:10AM +0200, Martin Sandve Alnæs wrote:
>> 2008/10/22 Johan Hake <hake@xxxxxxxxx>:
>> > On Wednesday 22 October 2008 11:12:02 Martin Sandve Alnæs wrote:
>> >> 2008/10/22 Johan Hake <hake@xxxxxxxxx>:
>> >> > On Wednesday 22 October 2008 10:17:43 Martin Sandve Alnæs wrote:
>> >> >> 2008/10/22 Johan Hake <hake@xxxxxxxxx>:
>> >> >> > On Wednesday 22 October 2008 09:32:31 Martin Sandve Alnæs wrote:
>> >> >> >> 2008/10/22 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> > On Tuesday 21 October 2008 23:23:27 Martin Sandve Alnæs wrote:
>> >> >> >> >> 2008/10/21 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> >> > On Tuesday 21 October 2008 22:34:04 Martin Sandve Alnæs wrote:
>> >> >> >> >> >> 2008/10/21 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> >> >> > On Tuesday 21 October 2008 21:37:13 Martin Sandve Alnæs wrote:
>> >> >> >> >> >> >> 2008/10/21 Anders Logg <logg@xxxxxxxxx>:
>> >> >> >> >> >> >> > On Tue, Oct 21, 2008 at 06:01:53PM +0100, Garth N. Wells
>> >> >
>> >> > wrote:
>> >> >> >> >> >> >> >> Anders Logg wrote:
>> >> >> >> >> >> >> >> > On Tue, Oct 21, 2008 at 04:45:01PM +0100, Garth N.
>> >> >> >> >> >> >> >> > Wells
>> >> >> >
>> >> >> > wrote:
>> >> >> >> >> >> >> >> >> I have a few questions and thoughts regarding the new
>> >> >> >> >> >> >> >> >> Function design
>> >> >> >> >> >> >> >> >>
>> >> >> >> >> >> >> >> >> * It's not clear to me what the intention is with
>> >> >> >> >> >> >> >> >> user-defined functions. The functions
>> >> >> >> >> >> >> >> >> Function::interpolate(...) never call eval(..), so
>> >> >> >> >> >> >> >> >> they can't pick up user-defined values. Should
>> >> >> >> >> >> >> >> >> Function::interpolate test for the presence of a
>> >> >> >> >> >> >> >> >> GenericVector to decide whether or not the Function is
>> >> >> >> >> >> >> >> >> discrete or user-defined?
>> >> >> >> >> >> >> >> >
>> >> >> >> >> >> >> >> > Yes, sorry. I've missed this. I'll fix it.
>> >> >> >> >> >> >> >> >
>> >> >> >> >> >> >> >> >> * It would be useful to declare user-defined functions
>> >> >> >> >> >> >> >> >> without associating a FunctionSpace. If we want to
>> >> >> >> >> >> >> >> >> interpolate the function, a FunctionSpace must then be
>> >> >> >> >> >> >> >> >> provided. Anyone see any problems with this?
>> >> >> >> >> >> >> >> >
>> >> >> >> >> >> >> >> > The reasoning here is that all Functions must always be
>> >> >> >> >> >> >> >> > associated with a FunctionSpace so that they may be
>> >> >> >> >> >> >> >> > correctly interpreted in forms and correctly plotted.
>> >> >> >> >> >> >> >> > When a Function is created in PyDOLFIN, it must always
>> >> >> >> >> >> >> >> > be associated with a certain FiniteElement (and in a
>> >> >> >> >> >> >> >> > while FunctionSpace). It would simplify the handling of
>> >> >> >> >> >> >> >> > Functions if they are always associated with a
>> >> >> >> >> >> >> >> > FunctionSpace.
>> >> >> >> >> >> >> >>
>> >> >> >> >> >> >> >> I agree that is makes life simple if every function has a
>> >> >> >> >> >> >> >> space, but it is a bit clunky for declaring user-defined
>> >> >> >> >> >> >> >> functions. The forms must be declared first to extract
>> >> >> >> >> >> >> >> the finite element to create the function space. Could
>> >> >> >> >> >> >> >> look nasty when a lot of functions are involved.
>> >> >> >> >> >> >> >>
>> >> >> >> >> >> >> >> We have a function Function::interpolate which takes a
>> >> >> >> >> >> >> >> function space V as an argument and it interpolates the
>> >> >> >> >> >> >> >> function u in V. What if we permit undefined function
>> >> >> >> >> >> >> >> spaces (which perhaps only have a domain)? We would then
>> >> >> >> >> >> >> >> interpolate the user defined function u in the provided
>> >> >> >> >> >> >> >> space V.
>> >> >> >> >> >> >> >>
>> >> >> >> >> >> >> >> Garth
>> >> >> >> >> >> >> >
>> >> >> >> >> >> >> > Are user-defined functions ever used without being related
>> >> >> >> >> >> >> > to a particular element/function space?
>> >> >> >> >> >> >> >
>> >> >> >> >> >> >> > It don't think it will be very clumsy. The clumsy thing
>> >> >> >> >> >> >> > will be to (in C++) get from something compiled by a form
>> >> >> >> >> >> >> > compiler to a FunctionSpace.
>> >> >> >> >> >> >> >
>> >> >> >> >> >> >> > If we can make that operation smooth, then creating
>> >> >> >> >> >> >> > (user-defined) functions will be very simple and
>> >> >> >> >> >> >> > convenient. One just needs to supply the variable V
>> >> >> >> >> >> >> > holding the function space.
>> >> >> >> >> >> >> >
>> >> >> >> >> >> >> > The current way of extracting function space data from the
>> >> >> >> >> >> >> > form is not very nice (in C++). What would be the optimal
>> >> >> >> >> >> >> > way to initialize a FunctionSpace in C++? We could think
>> >> >> >> >> >> >> > of extending the code generation to generate code that
>> >> >> >> >> >> >> > makes this convenient.
>> >> >> >> >> >> >> >
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> The current way of extracting function space data from the
>> >> >> >> >> >> >> form is not very nice in Python either, since it doesn't
>> >> >> >> >> >> >> work with compiled functions. (Never mind that the current
>> >> >> >> >> >> >> code is FFC-specific, this will be the same with UFL).
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> Using Python functors can easily make the assembly slower
>> >> >> >> >> >> >> than solving the linear system, so it's not really
>> >> >> >> >> >> >> interesting to do in real applications...
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> To make a function object that is both of a C++ subclass of
>> >> >> >> >> >> >> dolfin::Function and of the Python class ufl.Function, we
>> >> >> >> >> >> >> can't use the fixed multiple inheritance
>> >> >> >> >> >> >> solution in the current PyDOLFIN.
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> We would have to define a new class dynamically in python,
>> >> >> >> >> >> >> inheriting from both ufl.Function and the freshly compiled
>> >> >> >> >> >> >> C++ Function subclass. After all this work cleaning up the
>> >> >> >> >> >> >> Function class hierarchy, is that really something you want?
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> I'm not sure if that is even possible to do while
>> >> >> >> >> >> >> maintaining efficiency, with cross-language inheritance and
>> >> >> >> >> >> >> SWIG directors and all that.
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> If anyone has another solution, I'm very interested in
>> >> >> >> >> >> >> hearing it! Otherwise, I'm all for keeping the ufl.Function
>> >> >> >> >> >> >> objects used in form definition separated from
>> >> >> >> >> >> >> dolfin.Function objects used in assembly.
>> >> >> >> >> >> >
>> >> >> >> >> >> > I agree with Martin that we need to have a solution for
>> >> >> >> >> >> > PyDOLFIN users that does not depend on using python functors,
>> >> >> >> >> >> > as it will take forever for a complex form together with a
>> >> >> >> >> >> > moderate mesh to just assemble the form.
>> >> >> >> >> >> >
>> >> >> >> >> >> > Is it possible to let compile_functions compile a cpp
>> >> >> >> >> >> > function, with a FunctionSpace and all, instead of a mesh as
>> >> >> >> >> >> > it is today. Then after doing
>> >> >> >> >> >>
>> >> >> >> >> >> If you have a dolfin::FunctionSpace object already, there's no
>> >> >> >> >> >> reason compile_functions can't take this instead of
>> >> >> >> >> >> dolfin::Mesh. That's exactly the same and no problem at all.
>> >> >> >> >> >>
>> >> >> >> >> >> > this compile_function extract the element, and instantiate a
>> >> >> >> >> >> > UFL/FFC/PyFunction-function, and "attach" the compiled
>> >> >> >> >> >> > version to it. This
>> >> >> >> >> >>
>> >> >> >> >> >> What I state above is that this "attachment" must be done with
>> >> >> >> >> >> dynamic creation of a new class with multiple inheritance.
>> >> >> >> >> >> And I am unsure whether this will work out properly with SWIG
>> >> >> >> >> >> directors etc. I believe it _may_ work, but I don't dare to
>> >> >> >> >> >> keep my hopes up :-)
>> >> >> >> >> >
>> >> >> >> >> > Ok, I get it. For a moment I thought we could get away by
>> >> >> >> >> > defineing our own PyDOLFIN::Function class that could inherit
>> >> >> >> >> > from UFL/FFC, and then have a cpp_Function, but I realise this
>> >> >> >> >> > will not work.
>> >> >> >> >> >
>> >> >> >> >> >> See the attached python file for a prototype of dynamic class
>> >> >> >> >> >> creation with multiple inheritance using pure python classes.
>> >> >> >> >> >> (I think this is called "aspect oriented programming" by some
>> >> >> >> >> >> people)
>> >> >> >> >> >>
>> >> >> >> >> >> > can be used to define forms, but more important it can be
>> >> >> >> >> >> > handed to the python assembly that check if the function has
>> >> >> >> >> >> > a compiled version attached to it and send this to the
>> >> >> >> >> >> > cpp_assembler?
>> >> >> >> >> >>
>> >> >> >> >> >> If the "attachment" is anything other than inheritance, it will
>> >> >> >> >> >> have to be checked with manually written python code
>> >> >> >> >> >> _everywhere_ a dolfin::Function is expected... We can't have
>> >> >> >> >> >> one kind of functions for assembly and one for other stuff.
>> >> >> >> >> >
>> >> >> >> >> > Ok, I guess we have three different cases:
>> >> >> >> >> >
>> >> >> >> >> >  1) PyFunctions inherting from both UFL/FFC and cpp_Function as
>> >> >> >> >> > today, taking a functionsspace in its constructor. This will
>> >> >> >> >> > work with both user defined and discrete functions, more or less
>> >> >> >> >> > as we have it today.
>> >> >> >> >> >
>> >> >> >> >> >  2) The special functions, MeshSize, etc, can also be defined in
>> >> >> >> >> > the same way as now, right?
>> >> >> >> >> >
>> >> >> >> >> >  3) Using compile_functions, that creates a multi inheritance
>> >> >> >> >> > object that can be sent to any function expecting a
>> >> >> >> >> > cpp_Function, without manually extending the python interface.
>> >> >> >> >>
>> >> >> >> >> I'm with you up to this point.
>> >> >> >> >>
>> >> >> >> >> > Could the last be done by letting compile_function create a
>> >> >> >> >> > muliti inheritance Function. Instantiate the cpp_one with the
>> >> >> >> >> > function space and by that creating a dummy cpp_function. Then
>> >> >> >> >> > "attach" the compiled function to a protected attribute and
>> >> >> >> >> > define eval, by overloading it in python. This will then just
>> >> >> >> >> > call the attached and compiled cpp_functions eval.
>> >> >> >> >>
>> >> >> >> >> What you describe here sounds like the envelope-letter design
>> >> >> >> >> that was just _removed_ from dolfin.
>> >> >> >> >
>> >> >> >> > Yes, but only for compiled functions in Python. No other places.
>> >> >> >> >
>> >> >> >> >> What I'm suggesting is that
>> >> >> >> >> compile_functions dynamically creates a Python class that inherits
>> >> >> >> >> from ufl.Function and the freshly compiled C++ class, which is
>> >> >> >> >> a dolfin::Function subclass. Then it can construct an object of
>> >> >> >> >> this new class, passing a FunctionSpace object given by the user
>> >> >> >> >> to the dolfin::Function constructor, and an ufl.FiniteElement to
>> >> >> >> >> the ufl.Function constructor.
>> >> >> >> >
>> >> >> >> > This sounds doable. I realize now that this was what you were
>> >> >> >> > talking about in your previous emails, but I did not get it until
>> >> >> >> > now ;)
>> >> >> >> >
>> >> >> >> >> This of course requires that dolfin.FunctionSpace
>> >> >> >> >> is a Python subclass of dolfin::FunctionSpace with an additional
>> >> >> >> >> ufl.FiniteElement member variable. Using jit, dolfin.FunctionSpace
>> >> >> >> >> can compile the ufc::finite_element and ufc::dof_map classes it
>> >> >> >> >> needs from an ufl.FiniteElement. And then there's the issue of
>> >> >> >> >> reusing dofmaps, where DofMapSet enters the play...
>> >> >> >> >
>> >> >> >> > Do we need to jit compile ufc::finite_elements and ufc::dof_maps
>> >> >> >> > from the created ufl.FiniteElement? What about the one that follows
>> >> >> >> > from the FunctionSpace?
>> >> >> >>
>> >> >> >> I was thinking about when _constructing_ the FunctionSpace.
>> >> >> >> Just like PyDOLFIN uses jit in Function.__init__ today.
>> >> >> >
>> >> >> > Ok, something like:
>> >> >> >
>> >> >> > # Note pseudo code...
>> >> >> > class FunctionSpace(cpp_FunctionSpace):
>> >> >> >    def __init__(self,ufl_finite_element,mesh):
>> >> >> >        ufc_finit_element = jit(ufl_finite_element)
>> >> >> >        form = ufl.FiniteElement*ufl.TestFunction*ufl.dx
>> >> >> >        dof_map = jit(form)
>> >> >> >        cpp_FucntionSpace.__init__(mesh,ufc_FinitElement,dof_map)
>> >> >> >        self._UFL_FiniteElement = ufl_finite_element
>> >> >> >
>> >> >> >    def UFL_FiniteElement(self):
>> >> >> >        return self._UFL_FiniteElement
>> >> >> >
>> >> >> > By this the the ufc_element, ufl_element, the dofmaps and the mesh,
>> >> >> > are cached in the FunctionSpace.
>> >> >> >
>> >> >> > The Function would then be something like:
>> >> >> >
>> >> >> > class Function(cpp_Function,ufl.Function):
>> >> >> >    def __init__(self,function_space):
>> >> >> >        cpp_Function.__init__(function_space):
>> >> >> >        ufl.Function.__init__(function_space.UFL_FiniteElement())
>> >> >> >
>> >> >> > and dynamical created code in compile_functions()
>> >> >> >
>> >> >> > class MyFunction(MyCompiledFunction,ufl.Function):
>> >> >> >    def __init__(self,function_space):
>> >> >> >        MyCompiledFunction.__init__(function_space):
>> >> >> >        ufl.Function.__init__(function_space.UFL_FiniteElement())
>> >> >>
>> >> >> Something like that, yes. This is close to the current PyDOLFIN.
>> >> >>
>> >> >> But FunctionSpace might become a subclass of ufl.FunctionSpace
>> >> >> if we introduce that in UFL, and it should be possible to get
>> >> >> cached initialized and renumbered DofMaps from a DofMapSet.
>> >> >>
>> >> >> Since a DofMapSet will typically be initialized with a Form,
>> >> >> a Form depends on a Function, and a Function depends on
>> >> >> a FunctionSpace which should be initialized by the DofMapSet,
>> >> >> we have a cirular dependency right there.
>> >> >
>> >> > But won't you have this circular dependency in UFL already?
>> >>
>> >> In UFL this is simple:
>> >>
>> >>   FiniteElement depends on nothing
>> >>   Function depends on FiniteElement
>> >>   Form depends on Function
>> >
>> > I ment for a potential FunctionSpace class in UFL.
>>
>> Then it is simply:
>>
>>    FiniteElement depends on nothing
>>    FunctionSpace depends on FiniteElement
>>    Function depends on FunctionSpace
>>    Form depends on Function
>>
>>
>> >> The relation between UFL and UFC code is (at a certain level) simple:
>> >>
>> >>   ufc::* is generated by form compilers from ufl.* (equivalently ffc.*)
>> >>   once generated, ufc::* depends on nothing
>> >>
>> >> In DOLFIN (C++) it is also simple:
>> >>
>> >>   dolfin::* does not depend directly on ufl.* (equivalently ffc.*)
>> >>   dolfin::* depends on ufc::*
>> >>
>> >> PyDOLFIN could (should!) be kept simple. I've been using PyDOLFIN
>> >> without the FFC-dependent multiple inheritance stuff all the time
>> >> with the dolfin.cpp_* classes.
>> >
>> > Me too.
>> >
>> >> I believe that's also how it must be
>> >> done if you have external UFC code, e.g. precompiled in a library
>> >> or manually written.
>> >
>> > Agree.
>> >
>> >> So no matter what, it should at least be possible to avoid
>> >> relying on the multiple inheritance and JIT in PyDOLFIN.
>> >
>> > and more or less have the possibilities we have today. I agree.
>> >
>> > Johan
>>
>> Seems at least we're on the same page here.
>
> I've been in meetings the whole day and it took me a while to get
> through the discussion... :-)
>
> I agree that life would be easier if we did not need to worry about
> integrating UFL, UFC, and DOLFIN, but I think we should. I don't think
> a user should need to worry about the relations between all the
> different flavours of for example the finite element abstractions. We
> have quite a few:
>
>  ufl.FiniteElement
>  ufc::finite_element
>  dolfin::FiniteElement
>  dolfin.FiniteElement
>  FIAT.FiniteElement
>  (ffc.FiniteElement)
>  ...

There are different kinds of "worry about" though. Beginners
might only want to copy literally what the manual and demos do.
But I believe semi-technical users who want to debug their
applications by looking into the underlying code would find
clear separation between these types easier to comprehend.

In reality, they would only need to worry about:
  ufl.FiniteElement
  ufc::finite_element
  dolfin::FiniteElement == dolfin.FiniteElement


> DOLFIN plays two roles here:
>
> 1. It functions as a user interface/problem-solving environment and as
>   such should present a clean and simple interface where only one of the
>   above flavours should be visible.
>
> 2. It provides a general assembler, and as such should provide an
>   assembly algorithm that just works with ufc:: types (+ mesh and
>   linear algebra).
>
> So, we need to provide (1) for general use (and make it efficient) but
> also allow (2) for experts.
>
> So far, I like Martin's idea with dynamic inheritance to solve (1).
>
> --
> Anders

Do you have a solution for the circular dependency problem when reusing dofmaps?

I have a possible suggestion:

# First initialize Mesh and a matching DofMapSet with no dofmaps in it
mesh = dolfin.Mesh(...)
dofmapset = dolfin.DofMapSet(mesh)

# Define elements
element = ufl.FiniteElement(...) # or dolfin.FiniteElement?

# Define function spaces based on all the above
U = dolfin.FunctionSpace(mesh, element, dofmapset)
V = dolfin.FunctionSpace(mesh, element, dofmapset)

# dolfin.FunctionSpace.__init__ calls "dm = dofmapset.update(element)"
# Thus the first line (U = ...) results in a jit(element) call
# and the second line (V = ...) gets the previously constructed dolfin.DofMap.
# Is this ok with later renumbering, distribution etc of the dolfin.DofMap?

# Define forms based on function spaces
f = dolfin.Function(V) # both a dolfin::Function and ufl.Function
g = dolfin.Function(U)
a = f*dx # ufl.Form
b = f*g*dx

# Get all dofmaps used by forms
a_dofmaps = dofmapset.update(a)
b_dofmaps = dofmapset.update(b)


The backside here is that elements need to be compiled separately,
that is, they're not taken from the form. But there are relatively few
elements so the cache should hide that cost fairly well.

--
Martin