Direct access to the finite-element basis functions of a pre-specified FunctionSpace

Question

I need to get direct access to the finite-element basis functions (of Continuous-Galerkin type, degree 1) on a structured mesh over the unit square $\Omega := [0,1] \times [0,1]$.

I have made the following initializations:

from dolfin import *

nx = 3
ny = 3

mesh_obj = UnitSquareMesh(nx, ny)

V = FunctionSpace(mesh_obj, "CG", 1)

explicit_func = Expression(’1 + x[0]*x[0] + 2*x[1]*x[1]’)
explicit_func_proj = project(explicit_func, V)

and I would like to:

a) ask for confirmation if the mathematical description of the V (defined above) is "the set of $v \in C(\overline{\Omega})$ such that $v$ is affine on each triangle defined by mesh_obj with $v|_{\partial \Omega}$ is arbitrary"; if not, what is the correct mathematical description of the V which was defined above?

b) obtain a Python list with all finite-element basis functions $(\phi_j)_{j=1}^N$ for the vector space V (in order to be able to compute the quantities $\int_\Omega \phi_j \phi_k \, dx$ for all $j, k = 1, ..., N$)

c) compute the numbers $(c_j)_{j=1}^N$ such that explicit_func_proj equals $\sum_{j=1}^N c_j \phi_j$

I would appreciate it very much it if someone could advise me on these three issues.

Answer 1

You should be aware that while FEniCS follows an elegant symbolic approach, the underlying data structures for the matrices and vectors are in most cases still accessible. When you work through the tutorials this is explained in some detail. That being said, you never really need to have access to the basis functions for basic tasks. This is something which you may only need if you consider writing non-standard assemblers, e.g., for stabilizations which cannot be formulated using the UFL of for implementing things like limiters, postprocessings, etc..

Regarding your questions:

a) yes, the CG1 is simply the space of continuous piecewise linear functions. Boundary conditions are treated by adding the additional constraints in the assembly stage.

b) if you just want to obtain the entries of the mass matrix, why don't you just compute it using the usual workflow, i.e., assemble the form u*v*dx (u and v are trial and test functions, respectively) into a sparse matrix M and extract the components from there? You don't really need to have access to the local function spaces for that, neither do you have to assemble these integrals individually. Of course, you can create hat functions and compute the global integrals using assemble, but this is extremely inefficient compared to assembling everything in one sweep.

c) looks like you just want to access the components of an interpolation. This can again be done by accessing the underlying vector of a function.

I don't have FEniCS on the machine where I am writing this, but given these pointers you should be able to search for help. Everything I said above was covered in this Q&A sometime and some basic things like matrix assembly etc. are treated in the respective lessons of the tutorial.

Comments

Thank you for your suggestions.

In order to make sure that I've taken the correct message away, I would like to ask some clarifying questions about the following block of code, if I may:

from dolfin import *

nx = 3
ny = 3

mesh_obj = UnitSquareMesh(nx, ny)

V = FunctionSpace(mesh_obj, "CG", 1)

u = TrialFunction(V)
v = TrialFunction(V)

a = u*v*dx
# compute the stiffness/mass matrix associated with 
# the bilinear form a
A = assemble(a)
A_array = A.array()

explicit_func = Expression(’1 + x[0]*x[0] + 2*x[1]*x[1]’)
explicit_func_proj = project(explicit_func, V)
explicit_func_proj_arr = explicit_func_proj.vector().array()

d) If the dimension of the finite-element vector space V is $N$, is it true that $N$ equals len(explicit_func_proj_arr)?

e) If $(\varphi_j)_{j=1}^N$ is the finite-element basis for the vector space V, is it true that the value of $\int_{\Omega} \varphi_j \varphi_k \, dx$ is stored in A_array[j][k]?

f) Does explicit_func_proj_arr hold the numbers $(c_j)_{j=1}^N$, which were introduced in question c) from my initial post?

---

Three times yes, but be aware that the numbering of degrees of freedom must not necessarily coincide with the numbering of vertices in the mesh. There is however a mapping from vertices to DOFs available (search this site for details). I am also not sure if A.array (numpy) is dense or sparse... You should check this and maybe also have a look on how to access the linear algebra backend directly.

---

@Christian Waluga: (Overly delayed) Thanks for the feedback!

Regarding your comment that "the numbering of degrees of freedom must not necessarily coincide with the numbering of vertices in the mesh": does this refer to explicit_func_proj_arr?

If yes, is MiroK's answer from http://fenicsproject.org/qa/4045/dofmap-and-nodal-solutions?show=4045#q4045 what you had in mind (in order to create a two-dimensional numpy array X of shape $N$-by-$2$, where $N$ is the same as in d) from my post on August 17, so that the value of explicit_func at the point X[k] is stored in explicit_func_proj_arr[k] for $k = 1, ..., N$.)?