 # Displaying finite-element basis functions

I was hoping to display the basis elements for a degree 2 Lagrange space with a mesh element on the unit interval on got some unexpected behavior. I am hoping that someone can help provide clarity.

For a function u\in\mathcal{V}_{h}, I expect

u(x) =\sum_{i=0}^{\text{dim}(\mathcal{V}_{h})-1}u_{i}\phi_{i}(x)

where \lbrace\phi_{i}\rbrace_{i} is a basis for \mathcal{V}_{h} and u_{i} are the finite-element expansion coefficients. I assume that the expansion coefficients are represented in the vector data of a FEniCS Function. With this in mind I then attempted to construct \lbrace \phi_{i}(x)\rbrace as a Function, where only the $i$th entry of the associated vector data is nonzero.

import numpy as np
import dolfin as dl
import matplotlib.pyplot as plt
#alter the plot style
plt.style.use(‘classic’)
plt.rcParams.update({‘font.size’: 16})
#define a unit interval mesh with nElem elements
nElem = 1
mesh = dl.UnitIntervalMesh(nElem)
dl.plot(mesh)
plt.show()
#define degree deg continuous elements over the mesh
deg = 2
Ph = dl.FiniteElement(“CG”, mesh.ufl_cell(), deg)
Vh = dl.FunctionSpace(mesh, Ph)

#I now desire to visualize the basis functions
#I expect that the vector entry data of a dl.Function(), u
#corresponds to the finite element expansion of u
#that is u(x) = sum_i phi_i(x) u_i
#where u_i are the finite-element expansion coefficients
#and phi_i(x) are the finite-element basis functions

n = Vh.dim() # dim(Vh) = number of basis functions of Vh
basis_functions = [dl.Function(Vh) for i in range(n)]
ei = np.zeros(n) # numpy unit vector
for i in range(n):
ei[i] = 1.0
# assign the unit vector data to the dl.Function(Vh) function
# so that it now corresponds to a basis function
basis_functions[i].vector().set_local(ei[:])
ei[i] = 0.0
# plot the result
dl.plot(basis_functions[i])
plt.title(r’$k$th finite-element basis function k = {0:1d} '.format(i))
plt.show()

In the above example I specified degree deg 2 elements and so did not expect to see the following linear functions.

Please let me know if anything I have described is unclear.

I can’t understand your code, but is this what you mean?

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh, "CG", 2)
u = Function(V)
x = np.linspace(0, 1, 50)

for j in range(V.dim()):
u.vector()[:] = 0.0
u.vector()[j] = 1.0

y = np.fromiter(map(u, x), dtype=np.double)

plt.plot(x, y, label=f"$j = {j}$")

plt.xlim(0.0, 1.0)
plt.xlabel("$x$")
plt.ylabel(r"$\phi_j(x)$")
plt.legend()
plt.show()


E.g.    1 Like

Thanks so much @nate, it appears that I was just not properly visualizing the basis functions. I needed to sample the basis functions at more points, as you have via x=np.linspace(0, 1, 50). I wonder why the FEniCS plot function did not take care of that for me.

The fenics plot function only visualizes CG 1 functions (it interpolates higher order functions to CG 1). This is because the plot function Also supports plotting on triangles and tetrahedra.

1 Like

FEniCS’s built in plot interpolates at the mesh vertices because it’s cheap. The operation you see above is very expensive in comparison.

It’s often encouraged to output solutions using XDMFFile.write_checkpoint for use external visualisation software for this very purpose (e.g. paraview).

E.g. the following which are both visualisations of the same basis function; however, one the cheap interpolant at the mesh verticies, and the other with more expensive interpolant sampling throughout the element.  2 Likes

Thanks @dokken for explaining the reason why things occurred as they did. This is useful information.

Very good! Thanks so much! @nate