Hey guys, I am currently trying to evaluate the integral of the test functions over a line domain as you can see below

```
import numpy as np
import matplotlib
from ufl import FiniteElementBase
from fenics import *
from mshr import *
timesh = IntervalMesh(1,0.0,0.5)
Th = FunctionSpace(timesh, "Lagrange", 1)
vh = TestFunction(Th)
I1 = np.array(assemble(vh*dx))
print(I1)
```

It turns out that I1 returns a vector with the value of each Lagrange base function over the domain and then I simply sum the elements of the vector to get the integral of the sum. Now I am trying to get the value of the integral of the square of the sum, however when I do this

```
I1 = np.array(assemble(vh*vh*dx))
```

I get the error

```
Traceback (most recent call last):
(...)
File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 390, in __len__
raise NotImplementedError("Cannot take length of non-vector expression.")
NotImplementedError: Cannot take length of non-vector expression.
```

From what I could understand, the assemble function cannot read the size of the product of two forms. Is there a way to overcome this issue and get the integral of the products of all base functions? Is there a way to extract or access the entries of a TestFunction object?

Thanks in advance!!!

First, the wording of the text is a bit unclear, but, in mathematical notation, your first vector `I1`

is a vector whose i^\text{th} entry is

\int_\Omega\phi_i\,d\Omega\text{ ,}

where \phi_i is the i^\text{th} basis function of the space `Th`

and \Omega is the interval (0,0.5). It looks like what youâ€™re trying to compute next is a vector whose i^\text{th} entry is

\int_\Omega\phi_i^2\,d\Omega\text{ .}

This second vector is the diagonal of the matrix whose i,j entry is

\int_\Omega\phi_i\phi_j\,d\Omega\text{ ,}

which you can get from

```
import numpy as np
from fenics import *
timesh = IntervalMesh(1,0.0,0.5)
Th = FunctionSpace(timesh, "Lagrange", 1)
vh = TestFunction(Th)
uh = TrialFunction(Th)
A = assemble(uh*vh*dx)
diag_A = PETScVector(as_backend_type(A).mat().getDiagonal())
I1 = np.array(diag_A)
print(I1)
```

(There is some extra â€śwastedâ€ť computation of the off-diagonal entries, but it still scales linearly with the number of elements.)

1 Like

Sorry if I was not very clear. I am trying to compute \int_\Omega \phi_i\phi_j dx where the \phi_i's are the basis functions of the space Th. I think I got from your example an idea of what I needed. By what you did, The matrix A has entries a_{i,j} = \int_\Omega \phi_i\phi_j dx and thatâ€™s exactly what I need to compute the integrals. I donâ€™t need only the diagonal of this matrix, but all of its entries.

How do I get a numpy matrix of A?

I got the answer. When the object is of type

`dolfin.cpp.la.` `Matrix` ( **args* )`

we just simply use the attribute .array() and the matrix object becomes a numpy type matrix.