Function space model in python layer of dolfin is flawed

Issue #576 resolved
Martin Sandve Alnæs created an issue

In UFL there are various element types, e.g. FiniteElement and VectorElement.

In UFC all elements are represented behind a single ufc::finite_element/dofmap interface.

In DOLFIN/C++ there is only one FunctionSpace (ignoring MultiMeshFunctionSpace), no MixedFunctionSpace or VectorFunctionSpace etc.

In DOLFIN/Python, there is FunctionSpace, MixedFunctionSpace, VectorFunctionSpace etc. This is wrong. The python layer should wrap the C++ library, which has the concept that a FunctionSpace has a Mesh and a finite element. This also matches the new ufl.FunctionSpace(mesh,element). The element can be a vector element or mixed element, and dolfin.FunctionSpace should handle that.

In fact, you can create a function space in C++ using a VectorElement, and when passed to python it will not be a dolfin.VectorFunctionSpace but a dolfin.FunctionSpaceFromCpp. This means the dolfin.VectorFunctionSpace is only usable as a construction method, it cannot be relied upon for type checking etc. since it's not the only way to represent a vector element based function space.

I propose the following steps:

  • Introduce a constructor dolfin.FunctionSpace(mesh, element)

  • Change dolfin.FooFunctionSpace into factory functions and add deprecation warnings.

  • Merge classes dolfin.FunctionSpace and dolfin.FunctionSpaceFromCpp

Comments (25)

  1. Jan Blechta

    As a side note, I'd also like to see a deprecation of FunctionSpace(mesh, family, degree). It creates an overhead of constructing an unneeded dofmap in user codes (and our demo codes!):

    P2 = VectorFunctionSpace(mesh, 'Lagrange', 2) # its dofmap not needed anywhere
    P1 = FunctionSpace(mesh, 'Lagrange', 1) # its dofmap not needed anywhere
    TH = P2 * P1
    

    Moreover it tends to support a wrong idea that TH has "something" in common with P1 and P2. Construction from UFL element should be preferred.

    Sorry for polluting the thread but if there is a time to interface change, this can be solved alongside.

    It this interface change is not welcome, other solution would be a lazy construction of dofmap while keeping the FunctionSpace(mesh, family, degree) interface.

  2. Prof Garth Wells

    Changes proposed by @martinal sound good to me.

    @blechta :

    1. Your first point is a common cause of confusion, and would be good to fix (but the necessary interface change could cause some consternation).

    2. Lazy construction of dofmaps would be an easy fix for the efficiency issue.

  3. Jan Blechta

    The principal problem is that

    1. ufl.FunctionSpace uses __slots__
    2. We need to dynamically change __class__ of cpp.FunctionSpace instance to FunctionSpace (which should inherit from ufl.FunctionSpace).

    Presence of both 1. and 2. results in

    Traceback (most recent call last):
      File "space.py", line 5, in <module>
        W.sub(0)
      File "/home/jan/dev/fenics-master/src/dolfin/local.jan.fix-issue-576/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 259, in sub
        return FunctionSpaceFromCpp(cpp.FunctionSpace.sub(self, i))
      File "/home/jan/dev/fenics-master/src/dolfin/local.jan.fix-issue-576/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 367, in __new__
        cppV.__class__ = FunctionSpace
    TypeError: __class__ assignment: 'FunctionSpace' object layout differs from 'FunctionSpace'
    

    when running

    from dolfin import *
    mesh = UnitSquareMesh(3, 3)
    V = FunctionSpace(mesh, 'CG', 1)
    W = V*V
    W.sub(0)
    

    with jan/fix-issue-576. Any ideas, @martinal ? Do you really need slots? If yes, I'll try looking what can be done.

  4. Jan Blechta

    @martinal, anyway is it intended that ufl.FunctionSpace has a __dict__ (probably because of inheritance of nonslotted ufl.AbstractFunctionSpace)?

    import ufl
    e = ufl.FiniteElement('Lagrange', ufl.triangle, 1)
    V = ufl.FunctionSpace(ufl.triangle, e)
    V.__dict__
    
  5. Jan Blechta

    Ok, now I see that you've disable __slots__ for Coefficient so I've done the same in UFL jan/fix-dolfin-issue-576.

  6. Jan Blechta

    How the design should correspond to UFL where is special class MixedFunctionSpace (and TensorProductFunctionSpace). Or did you mean to keep MixedFunctionSpace class in PyDOLFIN?

  7. Martin Sandve Alnæs reporter

    Disabling slots for functionspace is fine.

    But is the design with dynamic class change necessary? From what I remember we don't do that with Function: dolfin.Function inherits from cpp.Function and ufl.Coefficient. It would be nice if we had only one pattern for how to join ufl and dolfin classes.

  8. Martin Sandve Alnæs reporter

    The classes ufl.MixedFunctionSpace and ufl.TensorProductFunctionSpace are not used yet. We might have to use a different name than MixedFunctionSpace to avoid confusion. The idea is that MixedElement is local and MixedFunctionSpace is global and can contain different meshes or meshviews in its subspaces. Function spaces on different meshes could be useful for the MultiMesh functionality that Anders is working on, but function spaces on different meshviews of the same mesh is something we will need, typically for coupled multiphysics problems.

    mesh = ...
    element1 = ...
    element2 = ...
    
    # Space of locally mixed elements
    M1 = FunctionSpace(mesh, MixedElement(element1, element2))
    
    # Space of globally mixed spaces over different parts of same mesh:
    meshview1 = MeshView(mesh, ...)
    meshview2 = MeshView(mesh, ...)
    V1 = FunctionSpace(meshview1, element1)
    V2 = FunctionSpace(meshview2, element2)
    M2 = MixedFunctionSpace(V1, V2)
    
    # Space of globally mixed spaces over different meshes:
    mesh2 = ...
    V1 = FunctionSpace(mesh, element1)
    V2 = FunctionSpace(mesh2, element2)
    M3 = MixedFunctionSpace(V1, V2)
    

    But we need to keep the MixedFunctionSpace notation (not the class) in dolfin at least for a while.

    One naming alternative is to follow the FiniteElement/MixedElement pattern and use FunctionSpace/MixedSpace for the new globally mixed class. However while VectorElement works ok as a name, VectorSpace is too loaded with other meaning.

  9. Jan Blechta

    But is the design with dynamic class change necessary? From what I remember we don't do that with Function.

    Function.__init_cpp_copy_constructor does that and there is the same purpose for it: wrap cpp.Function into Function. Thanks for pointing me to Function - now I can follow the pattern more closely.

    The name of the method is misleading - it does not behave like a copy constructor:

    from dolfin import *
    mesh = UnitSquareMesh(3, 3)
    V = FunctionSpace(mesh, 'CG', 1)
    u = cpp.Function(V)
    v = Function(u) # "copy" of cpp.Function
    w = Function(v) # copy of Function
    print u.vector().id(), v.vector().id(), w.vector().id()
    

    yields 5 5 7. I guess we need construction (not copy) from cpp.Function but the difference between two is quite subtle and prone to bugs.

  10. Jan Blechta

    As a side note, we could probably move Function[Space] to SWIG layer, merging with cpp.Function[Space] and subclassing ufl.(Coefficient|FunctionSpace). Then Function[Space] objects returned in C++ would already have correct type and no dynamic retyping would be needed. But this like a big rewrite so let's avoid it while not necessary.

  11. Jan Blechta

    Let's simply deprecate Function(Function) copy constructor and keep only Function(cpp.Function) wrapper constructor. There is Function.copy(deepcopy=False) for copying.

  12. Martin Sandve Alnæs reporter

    Yes, please fix that bug and deprecate the copy constructor. Is the wrapper constructor necessary? Maybe it can be a stand-alone function if it's used internally. Changing the way construction works tends to be hard to follow for developers other than the original author.

    But is it possible to inject the ufl base class in a .i file? I don't think I've seen that anywhere. If so it might not be much more work than moving the 'class Function' body to a .i file. It could have unforeseen side effects though.

    Currently the python interface of dolfin classes use some subclassing in .py files and some %extend in .i files.

    Another pattern is in meshes.py, where cpp.Mesh.init is 'patched' with a python override. That way UnitSquareMesh etc. actually use that same init function. I'm also not sure this is possible to do in a .i file. I don't particularly like this solution though.

  13. Martin Sandve Alnæs reporter

    Actually, you can just remove slots from all FunctionSpace classes, no need for the ufl_noslots marker outside of the Expr class hierarchy.

  14. Martin Sandve Alnæs reporter

    I suggest not enabling a deprecation warning for {Vector|Tensor|Mixed}FunctionSpace at least yet. The fallout of that operation is rather large and should be discussed more broadly. It's good that they're not classes but we need to keep the syntax to not break the world.

  15. Jan Blechta

    Don't we want to support count kwarg for ufl.FunctionSpace.__init__ as for ufl.Coefficient for a finer coupling between UFL count and DOLFIN id?

  16. Jan Blechta
    1. I'll put the deprecation of {Mixed|Enriched}FunctionSpace into separate pull request so that we can discuss it on next meeting with the world.

    2. I'll revert deprecation of {convenience_constructor|Vector|Tensor}FunctionSpace. I went really overboard here.

    Note 1. and 2. are really different story. TH = P2v * P1 is not semantically correct. We, in fact, do not support something like creation of TH from P2v and P1. It actually means TH = FunctionSpace(mesh, P2v.ufl_element() * P1.ufl_element()).

  17. Martin Sandve Alnæs reporter

    No, the FunctionSpace doesn't have a count in ufl. I don't see the need for one. The Mesh has a ufl_id (same thing, more descriptive) however, and mesh id + element should be sufficiently unique for the ufl/ffc/ufc side of things.

    1., 2. good.

    Yes, and ideally P2v * P1 should later become an actual MixedFunctionSpace where each subspace can have a different mesh(view), e.g. for including Lagrange multiplier on boundary as part of a coupled system.

  18. Log in to comment