Hi,

I have the bilinear form `q`

that can be represented as `q(f, g) = int2d(Th)( b * ( dx(uf) * dx(ug) + dy(uf) * dy(ug) ) )`

where `b`

is a fixed function, `uf`

and `ug`

solve a given BVP (details of the specific problem are not important) with Neumann boundary conditions `f`

and `g`

, respectively.

I also have a collection of functions `f1,...,fn`

and I want to compute the matrix `B`

given by `B_jk=q[fj,fk]`

for all pairs `j,k`

.

The question is, *what is the fastest way to do it?*

Using the expression for `q`

each time is extremely slow, so what I am doing now is to compute the matrix `A = q(Vh, Vh)`

and computing the inner products `<A * f[], g[]>`

whenever I want to compute `q(f, g)`

.

Thanks in advace!