Displaying finite-element basis functions

[tucker]

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.

[nate]

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.

[tucker]

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.

[dokken]

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.

[nate]

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.

[tucker]

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

[tucker]

Very good! Thanks so much! @nate

[dohertyw1]

Hi Nate.

Would you be able to share the code for these plots?

Best,
Will