Accessing basis functions

variational formulation
1

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!!!

2

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
3

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.

4

How do I get a numpy matrix of A?

5

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.