UFL: unsupported operand type(s) for *: 'Form' and 'Form'

Question

In my weak formulation, I have a term of the following form:

(v * f * ds(1)) * (u * f * ds(1))

where, as usual, u is the trialfunction and v the testfunction. ds(1) refers to the left boundary of a rectangular domain and f is a known function defined on this boundary, but I don't think that matters. When I try to run this code, I get the complaint:

TypeError: unsupported operand type(s) for *: 'Form' and 'Form'

Any suggestions of how to proceed?

Edit
Below is a small working example that produces the error. I plan to rewrite my own code so that it uses FEniCS for the finite element part, and this test is the last in a series to see if everything I would need is supported (everything else is), so any help would be much appreciated.

from dolfin import *

mesh = UnitSquareMesh(20, 20)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
leftboundary = LeftBoundary()

boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
leftboundary.mark(boundaries, 1)
ds = Measure('ds')[boundaries]

V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('x[1]*x[1]')
a = inner(grad(u), grad(v))
a += (v * f * ds(1)) * (u * f * ds(1))

Answer 1

You can reformulate the problem as a saddle point problem:

V = FunctionSpace(mesh, 'Lagrange', 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
u, r = TrialFunctions(W)
v, s = TestFunctions(W)
f = Expression('x[1]*x[1]')
a = inner(grad(u), grad(v))* dx  +  v * f * r *ds(1) \
  +      u * f * s * ds(1)       -   r * s * ds(1)

This results in a variational formulation that you can assemble, but it is symmetric indefinite so you have to adapt the linear solver accordingly.

Alternatively, if you want to preserve the SPD structure of your variational form, you can form the operator without explicitly constructing the matrix:

a = inner(grad(u), grad(v))* dx
b = f * v * ds(1)
p = a + u * v * dx

A = assemble(a)
B = assemble(b)
P = assemble(p)

class CompositeOperator(object):
    def __init__(self, A, B):
        self.A_mat = as_backend_type(A).mat()
        self.B_vec = as_backend_type(B).vec()
    def mult(self, mat, x, y):
        self.A_mat.mult(x, y)
        a = self.B_vec.dot(x)
        y.axpy(a, self.B_vec)

    def as_petscmatrix(self):
        from petsc4py import PETSc
        mat = PETSc.Mat().createPython(self.A_mat.getSizes(), 
                                       comm = mpi_comm_world())

        mat.setPythonContext(self)
        return PETScMatrix(mat)


A_B = CompositeOperator(A, B).as_petscmatrix()

solver = PETScKrylovSolver("cg", "hypre_amg")
solver.set_operators(A_B, P)

This solution requires that you have compiled dolfin with petsc4py.

Answer 2

In my weak formulation, I have a term of the following form:

(v * f * ds(1)) * (u * f * ds(1))

What are you trying to do with this line? I read it as
$$\int_{\partial \Omega_1} vf \, dx \int_{\partial \Omega_1} u f \, dx $$
which does not seem natural at all in a weak formulation.
Maybe you miswrote what you meant?

Comments

It might seem unnatural, but it is what I meant. The boundary condition to be fullfilled at $\partial \Omega_1$ is

$$\nabla u (\mathbf{r}) \cdot \hat{\mathbf{n}} = f(\mathbf{r}) \int_{\partial \Omega_1} f(\mathbf{r'}) u(\mathbf{r'}) d\mathbf{r'}$$

which, coming from Schroedinger's equation, results in this problematic form. I don't see an a priori reason why it should be problematic, as it is bilinear. From what I understand from UFL language, one should eventually obtain a summation of forms. So I tried:

 u * f * (v * f * ds(1)) * ds(1)

but that just changes the error to

ufl.log.UFLException: An Integral without a Domain is now illegal.

Or using version 1.4.0

ufl.log.UFLException: Calling is_cellwise_constant on a utility type is an error.

Answer 3

It was pointed out to me that there is no support for such multiplication of forms, but a possible workaround would be something along the lines of:

t1 = assemble(v * f * ds(1))
t2 = assemble(v * f * ds(1))
to_add = numpy.outer(t1.array(), t2.array())
a = a + to_add # won't work

Now I am trying to make this work, if possible avoiding any dense vectors/matrices like t1, t2, and to_add. So some questions I have are:

• How can I assemble sparse vectors directly, and save their outer product as a sparse matrix?
• How to add this second sparse matrix to_add to a (both with a different sparsity pattern)?

I will have to solve many systems, all with the same a, but for different to_add's. The sparsity pattern of all these to_add's will be the same though. So if I somehow can profit from defining this sparsity pattern once that would be great.

Thanks in advance

Comments

For future reference:

(...)
V = df.FunctionSpace(mesh, 'CG', 1)
u = df.TrialFunction(V)
v = df.TestFunction(V)
mat = lambda m: df.as_backend_type(m).mat()
R = df.FunctionSpace(mesh, 'R', 0)
r = df.TrialFunction(R)
s = df.TestFunction(R)
t1 = df.assemble(s * u * f * ds(i))
t2 = df.assemble(r * v * f * ds(i))
to_add = mat(t2).matMult(mat(t1))
A.axpy(1.0, to_add)
A = df.PETScMatrix(A)

r and s here merely make sure that a matrix is assembled rather than a vector. petsc4py is then used to do the matrix multiplication and add it to another petsc4py matrix A. I may have made a mistake in de lines above since it is a simplification of my own code, but the idea should be clear from this.