Lagrange finite element space

In this activity, we become familiar with a computational object that can represent elements of the simplest finite element space, called the Lagrange finite element space. These are vector spaces of functions on a spatial domain. We start from the case of a one-dimensional domain and proceed to higher dimensions.

1D: Lowest order one-dimensional case

Partition the interval $[0, 1]$ using a uniform mesh $\{ x_i = i h: i=0, \ldots, N \}$ of $N+1$ points, where the mesh size is $h = 1/N$ and set
$$ V_{h, 1} = \{ v: \; v \text{ is continuous and } v|_{[x_i, x_{i+1}]} \text { is linear }\}. $$ This vector space is called the lowest order Lagrange finite element space in the one-dimensional domain $[0,1]$. To work with functions in this space, we first make the mesh of points $x_i$.

import ngsolve as ng
from ngsolve.meshes import Make1DMesh

N = 5   # the number of mesh elements
mesh = Make1DMesh(5)

The resulting ngsolve mesh object mesh can be queried in many ways - type help(mesh) for the documentation. For example, you can double check it has the points you intended to set as the grid points:

for x in mesh.vertices:
    print(x.point)

(0.0,)
(0.2,)
(0.4,)
(0.6,)
(0.8,)
(1.0,)

In higher dimensions, the printed tuple x.point will show more coordinates.

The next steps use functionalities that are common for any finite element space (be it in one, two, or three dimensions) implemented in NGSolve. Let us make an ngsolve object representing the space $V_{h1}$ we introduced above.

Vh1 = ng.H1(mesh)

The syntax H1 for the Lagrange space comes from the fact that Lagrange finite elements are used to approximate weak formulations in the Sobolev space $H^1$. (More about that later.) Functions in this space are represented in ngsolve by GridFunction objects. Let us see how a function in this space looks like.

v = ng.GridFunction(Vh1, 'myfun')

This is now an uninitialized function in $V_{h1}$. We can set the values of v in various ways. One way is to set the values of an underlying array of numbers in v. Another way is to interpolate a known function, which we see first.

Interpolation

Interpolation into a finite element space is accomplished using the Set method of a GridFunction object. To help declare functions in terms of coordinates (just $x$ in 1D, $x, y, z$ in 3D etc), these coordinates are available as NGSolve CoefficientFunction objects. Let us import (the first and only coordinate in 1D) x to make up a function of $x$ and interpolate it into $V_{h1}$:

from ngsolve import x     

v.Set(x * x)

The GridFunction object v now has the interpolant of the function $f(x) = x^2$ into the piecewise linear space $V_{h1}$. To visualize v, we may use the common matplotlib python module as follows.

import matplotlib.pyplot as plt
import numpy as np

def plotv(v=v, sty='.-'):
    pts = [vtx.point[0] for vtx in mesh.vertices]
    plt.plot(pts, np.array(v.vec), sty)
    plt.xlabel('x'); plt.ylabel('$v$');

plotv()

Clearly, this is the expected continuous piecewise linear approximation of $x^2$ from the FE space $V_{h1}$.

Basis expansion

The functions $\psi_i \in V_{h1}$ with the property $$ \psi_i(x_j) = \delta_{ij} $$ at every mesh point $x_j$ are often called hat functions. Here $\delta_{ij}$ denotes the Kronecker delta (equals 1 if $i=j$ and 0 otherwise). Note that although the statement $\psi_i(x_j) = \delta_{ij}$ only gives the values of $\psi_i$ at the mesh vertices $x_j$, that is enough to determine the function $\psi_i$ everywhere, since it is linear in between the mesh points. Here is a visualization of one of these hat functions.

i = N//2
v.vec[:] = 0
v.vec[i] = 1
plotv()

And similarly, here are all of the hat functions associated to every mesh vertex:

for i in range(N+1):
    v.vec[:] = 0
    v.vec[i] = 1
    plotv()
plt.title('$\psi_i(x)$ for each $i$');

Notice how we made these hat functions using v.vec. To understand this, it is essential to understand how functions in $V_{h1}$ are represented using a basis expansion. A critical observation is the following statement.

Exercise: Prove that $\{\psi_i: i=0, \ldots, N\}$ defined above is a basis for the space $V_{h1}$.

In fact, when any $v \in V_{h1}$ is expressed in term of finite element basis of hat functions $\{\psi_i\}$, $$ v(x) = \sum_{i=0}^N v_i \psi_i(x), $$ the vector of coefficients $v_i$ in this basis expansion gives all information contained in v. This is precisely the vector we have been accessing and setting using the array v.vec. This is what allows functions to be represented using arrays of numbers.

2D: Lagrange spaces of higher degrees

Consider the degree $p$ Lagrange finite element space $$ V_{hp} = \{ v: \; v \text{ is continuous and $v|_K$ is a polynomial of degree } \le p \text{ on each mesh element}\}. $$ We examine the basis functions in the degree $p=2$ case. As in the 1D case, we begin by making a mesh.

from netgen.occ import unit_square
from ngsolve.webgui import Draw

mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.4))
Draw(mesh)

BaseWebGuiScene

Vh2 = ng.H1(mesh, order=2)  # order = p, the polynomial degree 
v = ng.GridFunction(Vh2)
Vh2.ndof

49

The number of basis functions of this space is given by the ndof attribute (which stands for number of degrees of freedom).

v.vec[:] = 0 
v.vec[15] = 1
Draw(v);

Open the controls and turn on deformation to see the appearance of a "hat" when this function is plotted as a three-dimensional graph. This is a "2D hat function" analogous to the 1D hat function we saw previously.

Because we set $p=2$ using the order=2 argument, this space also has quadratic functions within elements. Here is a continuous piecewise quadratic basis function of a different "shape". Indeed such basis functions, including the hat functions, are all collectively called shape functions in the finite element literature.

v.vec[:] = 0 
v.vec[30] = 1
Draw(v,
     settings={"Misc": {"subdivision": 15, "line_thickness": 5, "fast_draw": True},
               "Colormap":{"autoscale": True, "ncolors": 16}});

For discussion:

• How would you prove that the two shape functions shown above are linearly independent? What if their support overlaps?
• How would you make a 3D Lagrange space? (Look up documentation in ngsolve.org if needed.)

Before next class:

• Please review NGSolve i-Tutorial 1.2 to learn more about CoefficientFunction objects and ask questions if any.

---

$\ll$ Table Of Contents
$\ll$ Jay Gopalakrishnan