BUG: += does not assign when combining multiple Cofunctions

[pbrubeck]

Describe the bug

The Cofunction.__iadd__ method is not using assign when the update consisits of a (weighted) sum of Cofunctions. Instead it returns a FormSum, as opposed to the same operation on Functions, where we numerically compute the linear combination and accumulate the result in the LHS.

Steps to Reproduce

from firedrake import *

mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh, "CG", 1)

V = V.dual()
b = Function(V)
clist = [Function(V) for k in range(3)]
btype = type(b)

b.dat.data_wo[...] = 1
for c in clist:
    c.dat.data_wo[...] = 1

# this uses assign with Functions and Cofunctions
b += clist[0]
assert type(b) == btype, "b is now a %s" % type(b).__name__

# this only uses assign with Functions
b += clist[0] + clist[1]
assert type(b) == btype, "b is now a %s" % type(b).__name__

Expected behavior

+= should be consistent with Functions and Cofunctions. As these are well defined vector spaces, we have closure of linear combinations. I would expect += to use assign if I'm adding several Cofunctions

Error message

Traceback (most recent call last):
  File "/home/brubeckmarti/git/dphil_thesis/static_condensation/snippets/bug_cofunction.py", line 21, in <module>
    assert type(b) == btype, "b is now a %s" % type(b).__name__
AssertionError: b is now a FormSum

[Ig-dolci]

I am interested in this issue. I will dig into this.