Generation of a higher dimensional integral

Question

Hi everyone, I am trying to generate a linear problem $Ax = b$, where

$$ A = \int f_1 \otimes f_1 $$
$$ b = \int f_2 $$

where $f_1, f_2$ are some vector fields. I am currently assembling each coordinate

A[i,j] = assemble( f_1[i] * f_1[j] * dx )
b[i] = assemble ( f_2[i] * dx ) 

which is of course very slow when I use fine meshes. I would like to formulate this as a variational problem so that the underlying matrix handler (Petsc in my case) does everything faster, but my problem is that it assembly only works on scalar integrals. Is there a trick for this?

Thanks!

Comments

what are the FE spaces you're employing for your problem?

---

Well, only the FunctionSpaces for generating the functions $f_1, f_2$, which are (CG, 2) per coordinate. Each $f$ field has 6 entries, so the general problem consists of a 6x6 matrix A and a vector, 6x1, b. Right now I am using fenics for the integration that happens in each matrix and vector component. After assembling, I simply invert the $A$ matrix with numpy.

** I need fenics because the $f$ functions are a bit complicated, as they consist in the composition of a function which changes in each iteration. This system arrives from a linearized version of an inverse problem.

---

Do you need strict block structure?

---

Nope, only to find x given $f_1, f_2$.

---

Finally I solved the problem as follows: Define the space of constants (6 constants in this case)

el = VectorElement('R', triangle, 0, 6)
V = FunctionSpace(mesh, el)
b = Function( V )
u = TrialFunction(V)
v = TestFunction(V)

with this it is enough to write the mentioned forms as

A = dot(f1, u)*dot(f1,v)*dx
L = dot(f2, v)*dx

and finally solve as

solve( A==L, b)

To update a current solution (this problem comes from a Newton-like formulation), do

a += b.vector()

It was as simple as getting the weak formulation of the problem $Ax = b$, which was just seeing the equations as per coordinate. Best regards!