Subfunction assignment

[Johan Hake]

Given the following set of Functions:

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)
Vv = VectorFunctionSpace(mesh, "CG", 1)
W = V*Vv
u0 = Function(V)
u1 = Function(Vv)
u2 = Function(W)

It would be very nice if all of these assignments were allowed:

u0.assign(u1.split()[0])
u1.assign(0, u0)
u2.assign([1,0], u1.split()[0]))

Right now the first case, where all values of the first sub function of u1 is assigned to the whole of u0, is supported, by collapsing the first sub function of u1. In the second case all values of u0 are assigned to the first sub function of u1. In the third case all values of the first sub function of u1 is assigned to the first sub function of the second sub function of u2.

To accomplish the second and third case I suggest extending the C++ interface with:

Function::assign(const std::vector<std::size_t>& component, const Function& v);

which in its body iterate over all cells and tabulates dofs and call GenericVector::get and GenericVector::set with the proper indices and values. The correct indices are acquired by calling cell_dofs of the correct sub_space acquired from the component argument.

As we iterate over cells one could extend the interface to also take a MeshFunction and a marker, making it possible to assign values associated to a subdomain too.

[Patrick Farrell]

Yess!! A big upvote for this.

The syntax I had in my head was

u2.assign((u0, u1))

for the simple case where you have all the subfunctions that you want to pack, and

u2.assign({0: u0, 1: u1})

for the more complicated one where you might only want to pack some subfunctions. But so long as it is pretty, I don't mind ..

[Johan Hake]

You suggest a fourth option, which I think should be possible to cover if we change the syntax to:

Function::assign(std::map<const std::vector<std::size_t>,  boost::shared_ptr<const GenericFunction> >);

Note that it now includes a GenericFunction instead of just Function, as I figured one can just call GenericFunction::restrict to get the needed cell dofs.

[Martin Sandve Alnæs]

From python I'd like to do things that are easy to implement with a standalone function on top of a sufficiently generic C++ interface:

# Instead of u0.assign(u1.split()[0]):
assign(u0, u1[0])
# Instead of u1.assign(0, u0):
assign(u1[0], u0)
# Instead of u2.assign([1,0], u1.split()[0])):
assign(u2, (u1[1], u1[0]))

# And maybe something like this:
assign(u2[0], u1[1], domain=(meshfunction, domainid))

Might be possible in C++ as well, but at least in python this should be easy if the C++ side you describe is there.

[Johan Hake]

Why would you want that as a standalone function? Why not just:

u0.assign(u1[0])
u1.assign({0:u0})
u2.assign(u1[1], u1[0])

which could easily be handled by the python layer of Function.assign. I guess your approach is cleaner? Especially #2 would is nicer in your example I think, as one would not be able to ever do:

u1[0].assign(u0)

right?

[Martin Sandve Alnæs]

I want to treat these the same:

u0 = Function(V)
u1 = Function(V)
u = as_vector((u0,u1))

v = Function(Vv)

I.e. make these two possible:

assign(u, v)
assign(v, u)

This has applications in our segregated Navier Stokes solvers where the unknowns live in separate scalar spaces.

[Prof Garth Wells]

Some of Johan's suggestions are straightforward, and some are not. Assignment of a constant and assignment of Functions that share the same FunctionSpace are easy. For

V = FunctionSpace(mesh, "CG", 1)
Vv = VectorFunctionSpace(mesh, "CG", 1)
W = V*Vv

Vv.sub(0) != V, which makes u1.assign(0, u0) non-trivial. It will require a re-mapping of dofs.

[Johan Hake]

For this to be robust I envision that the C++ function would implement something like:

# Ensure that the functionspaces has the same element.
for cell in cells(mesh):
    self_dofs = self_dofmap.cell_dofs(cell.index())        
    other_dofs = other_dofmap.cell_dofs(cell.index())
    other_u.vector().get(values, other_dofs)
    self_u.vector().set(values, self_dofs)

To cover Martins and Patricks suggestions of having several other functions, we just need to grab the correct subdofmaps for the receiving function.

It might not be the most effective way of assigning values but it should be robust, using abstractions which are used in the assemble algorithm.

[Prof Garth Wells]

The above operations are potentially very expensive.There is no guarantee that the FunctionSpaces are partitioned in the same way across processes.

[Prof Garth Wells]

I would suggest in the first instance to not support assignment with non-matching FunctionSpaces and encourage users to use subspaces properly, ie. do

Vv = VectorFunctionSpace(mesh, "CG", 1)
V = Vv.sub(0)
u0 = Function(V)
u1 = Function(Vv)

instead of

V = FunctionSpace(mesh, "CG", 1)
Vv = VectorFunctionSpace(mesh, "CG", 1)
u0 = Function(V)
u1 = Function(Vv)

[Jan Blechta]

Is not this protocol suggested by Garth sufficiently robust to cover all the relevant use cases (except loading external data) ?

[Johan Hake]

It would be nice to support both. A typically use case is that you have a split solver where one solution should be interwoven into another. Then one would start with non matching FunctionSpaces. Also your code above raises an error:

Error:   Unable to create function.
Reason:  Cannot be created from subspace. Consider collapsing the function space.

Also a slow but robust solution is better than nothing.

Rather than calling get/set for each cell one could instead build a mapping between all dofs. Then call get and set only onces for each (sub)function assignment. But I am not sure how one would deal with duplicated dofs.

Such mapping might also be expensive to build. But once it is built an assignment would be just a call to get/set. Maybe create a dedicated assignment structure that cached such mapping?

[Patrick Farrell]

The mapping [index of sub-function-space] -> [dofs of the big function space associated with that sub-function-space] is also exactly what you need to implement FIELDSPLIT efficiently, right?

[Martin Sandve Alnæs]

Surely this cannot be more expensive than a projection, which is the current workaround for situations like this today.

Does collapsed function spaces preserve dof ownership? Is there a way to check if a space is originally a subspace of another after it is collapsed?

[Prof Garth Wells]

Whether the operation is more or less expensive than a projection is not the point. The point is that it should be clear to a user if a potentially expensive operation is being performed. We should also, when possible, have a design that helps/encourages users to write code that will be efficient.

It is not presently possible to check which space a collapsed dofmap came from. We could use a boost::weak_ptr to hold a pointer to the parent space. This would in most cases allow the parent to be checked without introducing the possibly large memory overhead of using a shared_ptr.

[Prof Garth Wells]

Patrick: the functionality in the difficult case is not the same as what is required for field split because in the case of field split we do not work with different FunctionSpaces, but with subspaces that all come from a common parent. The challenging case is with non-matching functions spaces.

[Chaffra Affouda]

In the mean time is there a work around for this? I was thinking something like,

Q = FunctionSpace(mesh,'CG',1)
M = MixedFunctionSpace([Q,Q,Q])
u_mixed = Function(M)
s0 = M.sub(0).dim()
u_mixed.vector()[0:s0] = 1.0 # to assign 1.0 to u_mixed[0]

I tried it but it doesn't seem to work.

[Nico Schlömer]

Is there already a feature branch for this? At the moment, I'm going great lengths to accomplish something like

u2.assign(0, u2[0] + u1)

aka

u2[0] += u1

[Johan Hake]

@chaffra : No, the dofs are no longer ordered by subspaces. @nschloe : No

[Johan Hake]

Should now be fixed in: 42cb6f88c9034deab27bb5dc24a6a894c403377d

[Simon Funke]

Thumbs up for this very useful feature! Thanks.

[Prof Garth Wells]

Removing milestone: 1.3 (automated comment)