Vector elements & Betti numbers

We have already seen the Lagrange and the DG finite element spaces. They are both spaces of scalar-valued functions. How can we approximate solutions of boundary value problems that are vector fields?

We could certainly use a Cartesian product of scalar finite element spaces to approximate a vector field component by component. However, distinctly better methods can often be constructed using other vector finite element spaces, such as the Raviart-Thomas space or the Nedelec space, both spaces of vector fields we shall see in this notebook.

We will expand on their utility in certain boundary value problems in later lectures. For now, the goal is to introduce you to computing with them and to show you that there might be more than just utilitarian reasons for studying these finite elements. They are connected through rich and beautiful discrete structures. As a prelude to such things, we will see in this notebook how a topological invariant emerges when considering curl as an operator from the Nedelec space in two space dimensions.

The Raviart-Thomas space

$\newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} \newcommand{\dive}{\mathop{\mathrm{div}}} \renewcommand{\div}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \newcommand{\R}{\mathbb{R}} $ The Raviart-Thomas (RT) finite element space on a two-dimensional mesh $\oh$ (of a domain $\om$) is defined by $$ R_{h, p} = \left\{ q\in H(\dive, \om) : \; \text{ for all elements } K \in \oh, \quad q|_K \in P_p(K)^2 + \begin{pmatrix} x \\ y \end{pmatrix} P_p(K) \; \right\} $$ for any $p \ge 0$. Here $P_p(K)$ denote the set of polynomials of degree at most $p$. The distributional divergence of a piecewise smooth vector field was calculated in an earlier lecture, from which it follows that vector fields $q$ in $R_{h, p}$ satisfy ⟦ $q\cdot n$ ⟧ $=0$, i.e., the jump of the normal component of $q$ across element interfaces vanishes. Therefore, plots of any function in the RT space should show a continuously varying normal component across the interior mesh facets.

The continuity of the normal component is a useful property to have when we know that the function being approximated has that same normal continuity property. Often, when $q: \om \to \R^2$ is a vector field representing a physical flux, its normal component remains smooth across a material interface even when its other components jump across the interface.

Approximation

As an example, consider approximating the following vector field using the RT space:
$$ q = \left\{ \begin{aligned} & (\sin(x), y), && x \le 1/2, \\ & (\sin(x), 1-y), && x > 1/2. \end{aligned} \right. $$ Across the interface $x=1/2$, its normal component (the $x$-component in this case) is continuous, while its tangential component is discontinuous. Let's plot $q$ after making two subdomains to the left and right of the $x=1/2$ line.

import ngsolve as ng
import numpy as np
from ngsolve.webgui import Draw
from netgen.occ import OCCGeometry ,WorkPlane, Rectangle, Glue
from ngsolve import x, y, sin, cos, GridFunction, div, curl

wp = WorkPlane()  # Make left and right subdomains
lft = wp.Rectangle(1/2, 1).Face()
rgt = wp.MoveTo(1/2,0).Rectangle(1/2, 1).Face()
lft.faces.name = 'lft'; rgt.faces.name = 'rgt'

geo = Glue([lft,rgt])
Draw(geo)
mesh = ng.Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=1/4))

q = mesh.MaterialCF({'lft': (sin(x), y),   # Piecewise q on the left & right subdomains
                     'rgt': (sin(x), 1-y)})

The normal component (which is the $x$ component in this case) of $q$ is smooth:

Draw(q[0], mesh);

The tangential component of $q$ is discontinuous across the middle interface.

Draw(q[1], mesh);

We can also plot $q$ as a vector field.

Draw(q, mesh, vectors={'grid_size': 15});

If we try to interpolate such a function using a product of Lagrange spaces, then the discontinuity at the interface will generate large interpolation errors. Instead, knowing that the normal component is continuous, we can try to interpolate the function into the RT space. This will liberate the tangential component of the approximation from continuity constraints, while preserving the normal continuity.

The implementation of the RT space $R_{h, p}$ in NGSolve is accessed using ng.HDiv. The Set method, as in other finite element spaces in NGSolve, can be used to perform the interpolation.

p = 5   # This corresponds to the p in the RT definition above.
R = ng.HDiv(mesh, order=p, RT=True)

qh = GridFunction(R)
qh.Set(q)
print('Error in RT interpolation =', 
      ng.sqrt(ng.Integrate((q - qh)**2, mesh)))

Error in RT interpolation = 2.431977823680355e-11

Compare this to what happens when we interpolate the same function in the product of two Lagrange finite element spaces.

V2 = ng.VectorH1(mesh, order=p)  # Cartesian product of Lagrange spaces
qq = GridFunction(V2)
qq.Set(q)
print('Error in product Lagrange space interpolation =', 
      ng.sqrt(ng.Integrate((q - qq)**2, mesh)))

Error in product Lagrange space interpolation = 0.11090298269878275

Clearly the error in the product space approach is several orders higher than that of the error produced by the RT space. The source of the large error while interpolating using ng.VectorH1 is immediately evident if you plot the $y$-component of the interpolant qq or the error. The interpolant, being continuous in both components, has a difficult time approximating the discontinuous component.

Draw(qq - q, mesh, 'Lagrange interpolation error', vectors={'grid_size': 25});

Shape functions

It is illustrative to visualize the basis functions of the RT space. We had gotten used to speaking of "hat functions" but that term is not useful for general finite elements. Basis functions generated from degrees of freedom of any abstract finite element spaces are called shape functions. Here is a visualization of a shape function of the lowest order RT space (the $p=0$ case).

R = ng.HDiv(mesh, order=0, RT=True)

shapenumber = 15
shape = ng.GridFunction(R, name='shape')
shape.vec[:] = 0
shape.vec[shapenumber] = 1
Draw(shape, mesh, vectors={'grid_size': 20});

You can see any shape function you wish by revising the shapenumber in the code above. In all cases, you see that

• the shape functions are vector fields whose normal components are continuous,
• they are supported on a patch of one or two elements,
• their normal component is zero on all mesh facets except one.

Thus in the lowest order case, each shape function is associated to a single mesh facet (an edge, in our 2D mesh). If this edge is on the domain boundary, the shape function is supported on one triangle. If the edge is an interior edge, the shape function is supported on the two triangles which share the interior edge.

One can also directly see the connection between the number of mesh facets and the dimension of the lowest order RT space by querying both numbers:

mesh.nfacet

50

R.ndof

50

Of course, for higher $p$, a number of additional shape functions must be added. You can easily visualize the higher order shape functions by tweaking the above few lines of code.

Convergence of divergence

Ordinarily, when a piecewise polynomial interpolant of a function $q$ converges at a rate $O(h^k)$, we expect its derivatives to converge at a lower rate.

However, somewhat miraculously, the Raviart-Thomas interpolant $I_h^{RT} q$ has the property that $\dive(I_h^{RT} q - q)$ converges at the same rate as $I_h^{RT} q - q$. To obtain this interpolant in NGSolve, we use the Set method with dual=True option. We then compute the errors on successively refined meshes and tabulate the convergence rates much the same way as in a previous notebook.

def InterpolateOnSuccessiveRefinements(q, divq, mesh0, p=0, nrefinements=8):                                   
    """Error in RT interpolant on a sequence of uniformly refined meshes."""
      
    errors = []; diverrors = []; 
    mesh = ng.Mesh(mesh0.ngmesh.Copy())
    
    for ref in range(nrefinements): 
        RT0 = ng.HDiv(mesh, order=p , RT=True)
        qh = GridFunction(RT0)
        qh.Set(q, dual=True)   # Gives the canonical interpolant
      
        err = ng.sqrt(ng.Integrate((q - qh)**2, mesh))
        diverr = ng.sqrt(ng.Integrate((divq - div(qh))**2, mesh))
        errors.append(err); diverrors.append(diverr)      
        mesh.ngmesh.Refine()

    return np.array(errors), np.array(diverrors)

try:   # copied over from a prior notebook
    from prettytable import PrettyTable
except ModuleNotFoundError:
    try:
        import micropip
        await micropip.install("prettytable")
        from prettytable import PrettyTable
    except ModuleNotFoundError:
        !pip3 install prettytable
        from prettytable import PrettyTable
      
def TabulateRate(name, dat, h0=1):
    col = ['h', name, 'rate'];  h0col = ['%g'%h0]
    t = PrettyTable()
    t.add_column(col[0], h0col + [h0col[0] + '/' + str(2**i) for i in range(1, len(dat))])
    t.add_column(col[1], ['%.12f'%e for e in dat])
    t.add_column(col[2], ['*'] +  ['%1.2f' % r for r in np.log(dat[:-1]/dat[1:])/np.log(2)])
    print(t)

To obtain the convergence tables, we call these functions with the above set discontinuous vector field $q$. By the formula for distributional divergence of discontinuous fields from the lectures, we find that $\dive(q)$ has no edge distributions and can simply be given piecewise as follows:

divq = mesh.MaterialCF({'lft': cos(x)+1, 'rgt': cos(x)-1})

Using q and divq, we compute the errors and tabulate the rates:

err, diverr = InterpolateOnSuccessiveRefinements(q, divq, mesh)
TabulateRate('||q - IhRT(q)||', err, h0=1/4)
TabulateRate('||div(q - IhRT(q))||', diverr, h0=1/4)

+----------+-----------------+------+
|    h     | ||q - IhRT(q)|| | rate |
+----------+-----------------+------+
|   0.25   |  0.079574736098 |  *   |
|  0.25/2  |  0.039789854089 | 1.00 |
|  0.25/4  |  0.019895238575 | 1.00 |
|  0.25/8  |  0.009947658253 | 1.00 |
| 0.25/16  |  0.004973833998 | 1.00 |
| 0.25/32  |  0.002486917608 | 1.00 |
| 0.25/64  |  0.001243458880 | 1.00 |
| 0.25/128 |  0.000621729450 | 1.00 |
+----------+-----------------+------+
+----------+----------------------+------+
|    h     | ||div(q - IhRT(q))|| | rate |
+----------+----------------------+------+
|   0.25   |    0.030651751056    |  *   |
|  0.25/2  |    0.015370264478    | 1.00 |
|  0.25/4  |    0.007690652409    | 1.00 |
|  0.25/8  |    0.003846015344    | 1.00 |
| 0.25/16  |    0.001923093787    | 1.00 |
| 0.25/32  |    0.000961557657    | 1.00 |
| 0.25/64  |    0.000480780174    | 1.00 |
| 0.25/128 |    0.000240390255    | 1.00 |
+----------+----------------------+------+

Clearly both $\| q - I_h^{RT} q\|_{L^2(\Omega)}$ and $\| \dive(q - I_h^{RT} q)\|_{L^2(\Omega)}$ converges at the same rate! In a later lecture, we will prove that this superconvergence property of the divergence of the interpolation error when using the RT space arises from a special commuting diagram property involving the RT space.

The 2D Nedelec finite element space

Consider the curl operator in two space dimensions, acting on a vector field $u: \om \to \R^2$, $$ \newcommand{\curl}{\mathop{\mathrm{curl}}} \curl u = \partial_x u_1 - \partial_y u_0. $$ You have worked with $\curl$ as an operator acting on three-dimensional (3D) vector fields. The above two-dimensional (2D) version is the obtained as the $z$-component of the 3D curl when applied to a vector field that has only $x$ and $y$ components with no $z$ dependence. It takes a 2D vector field $u$ and produces a scalar field $\curl u$.

The 2D curl is intimately related to the 2D divergence. To see this, let $J_{\pi/2}$ denote the operator that rotates a vector clockwise by 90 degrees, i.e., $$ J_{\pi/2} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} b \\ -a \end{pmatrix}. $$ Now, if two vector fields $u$ and $q$ are related by
$$ u = J_{\pi/2} q $$ then obviously $$ \dive q = \curl u. $$

Moreover, if $t$ denotes the counterclockwise unit tangent vector on element boundaries, then it is related to the unit outward normal $n$ at the same boundary point by $$ J_{\pi/2} n = t. $$ Hence $$ u\cdot t = q \cdot n. $$ Thus the condition that ⟦ $q\cdot n$ ⟧ $= 0$ is equivalent to the condition that the jump of the tangential component of $u$ vanish, ⟦ $u\cdot t$ ⟧ $= 0.$

We define the Nedelec space in two dimensions as the RT space rotated: $$ N_{h, p} := J_{\pi/2} R_{h,p}. $$ It consists of piecewise polynomial vector fields that are tangentially continuous. Applying the rotation operator $J_{\pi/2}$ to the polynomial functions in our previous definition of $R_{h, p}$, we also find that $N_{h, p}$ can be equivalently described by $$ N_{h, p} = \left\{ u\in H(\curl, \om) : \;\text{ for all elements } K \in \oh, \quad u|_K \in P_p(K)^2 + \begin{pmatrix} -y \\ x \end{pmatrix} P_p(K) \; \right\}. $$

By our prior discussion relating the 2D curl and div, we also conclude that $N_{h, p} \subset H(\curl, \om)$ since $R_{h,p} \subset H(\div, \om)$. Note that this rotational isomorphism between $H(\dive, \om)$ and $H(\curl, \om)$ does not hold in 3D. The Nedelec space in 3D is a truly different space, not isomorphic to the 3D RT space. Postponing the 3D case for later discussions, let us proceed with the 2D case. NGSolve provides an implementation of the Nedelec space accessible through HCurl as follows.

N = ng.HCurl(mesh, order=0, type1=True)  # This is N_{h,0}

shapenumber = 15
shape = GridFunction(N, name='shape')
shape.vec[:] = 0
shape.vec[shapenumber] = 1
Draw(shape, mesh, vectors={'grid_size': 30});

As you can see, this shape function visually appears to be exactly the same shape function plotted earlier for the RT case, except for a rotation by 90 degrees. Note how the tangential component of the shape function varies continuously across the edges. Just as the RT space was useful for interpolating vector fields with normal continuity, the Nedelec space is useful for interpolating vector fields with tangential continuity.

The Nedelec finite element is sometimes called the edge finite element in some engineering literature. Interestingly, the shape functions of the RT and the Nedelec element were known in another mathematical context (they appeared in a 1957 book by Whitney) long before these elements were known. In recognition of this, mathematicians sometimes refer to these shape functions as Whitney forms. Here is a formula for the shape functions (of the lowest order Nedelec space) that you just saw in the above visualization:

Range of the 2D curl

Application of the 2D curl operator to a function in the Nedelec space $N_{h, p}$ results in a function with no continuity restrictions across element interfaces in general. Since the application of curl also reduces the polynomial degree within each element, the result is a piecewise polynomial of degree not more than $p$, i.e., a function in the DG space $P_p(\oh)$. Thus, we may view $\curl$ as an operator from the Nedelec space to the DG space, $$ \curl : N_{h, p} \to P_{p}(\oh). $$ How much of the $P_p(\oh)$ is filled with the range of $\curl$?

To answer this, we will need to understand the range of the above map. Since $\curl$ may act on smooth function spaces or on finite element spaces like $N_{h, p}$, it is a good idea to indicate the specific domain we consider. So we will write $ \text{`` }\curl: N_{h,p} \text{ ''} $ to indicate the $\curl$ when it's viewed as an operator acting (only) on $N_{h, p}$ functions. We now proceed to compute the dimension of its range, i.e., the number $$\newcommand{\rank}{\mathop{\mathrm{rank}}} \rank(\curl: N_{h, p}). $$ To do so, we will need a representation of the operator $\curl: N_{h,p}$ as a matrix obtained using the basis functions of $N_{h,p}$ and $P_p(\oh)$. Of course, we will also need a domain and a mesh. Let us punch some holes into a square to make the domain a bit more interesting than just a square and proceed.

wp = WorkPlane()
sqr = wp.Rectangle(2, 2).Face()
rgthole = wp.MoveTo(1.5,1).Rectangle(1/6, 1/6).Face()                 
lfthole = wp.MoveTo(0.5,1).Rectangle(1/6, 1/6).Face()      
dom = sqr-rgthole -lfthole
Draw(dom);
mesh = ng.Mesh(OCCGeometry(dom, dim=2).GenerateMesh(maxh=1/4))

The Nedelec and DG spaces on this mesh for a given degree $p$ are set as follows.

# To get N_{h, p} for p=0,1,2,3,... use order=0,2,3,4,... in HCurl:

p = 0    
N = ng.HCurl(mesh, order=p+1-max(1-p,0), type1=True)
D = ng.L2(mesh, order=p)

Now let us think about how to make a matrix representation of the operator $\curl: N_{h, p} \to P_p(\oh).$

One way to obtain this matrix is to take each shape function in $N_{h, p}$, apply curl, and then expand the result in terms of the shape functions of the DG space. However, since the application of curl in NGSolve does not output an object in the DG space, we would need to do an additional step of using the Set method of the DG space: if you provide it a function that you know should be contained within the DG space, then the Set method interpolates it with zero error, creating the interpolant as a GridFunction object. Such objects can be queried for their vector of coefficients of basis expansions, which become entries of the curl matrix. An implementation would look like this:

def curlmatrix(N, D):
    """Return the matrix representation of curl: N -> D"""
    curlmat = []
    shapeN = GridFunction(N)
    curlD = GridFunction(D)
    for i in range(N.ndof):
        shapeN.vec[:] = 0;  shapeN.vec[i] = 1   # i-th shape function 
        curlD.Set(curl(shapeN))                 # interpolation error is 0
        curlmat.append(list(curlD.vec)) 
    return np.array(curlmat).T

However, the code above is inefficient. As you know, python loops are slow. Moreover, inside the above loop, the Set does a global computation when what was really needed is just a local computation. Thankfully, NGSolve provides a ConvertOperator which accomplishes the same thing faster.

from ngsolve import curl 
from ngsolve.comp import ConvertOperator
from scipy.sparse import coo_matrix

# Compute the matrix representation of curl: N -> D
u = N.TrialFunction()
Crl = ConvertOperator(N, D, trial_proxy=curl(u))
i, j, val = Crl.COO()
curlmat = coo_matrix((val, (i, j))).toarray()

Now that we have the matrix representation curlmat of the linear operator $\curl: N_{h, p}$, computing its rank can be done by any number of linear algebra techniques. Here is one way, which proceeds by counting the number of nonzero singular values.

from scipy.linalg import svdvals

s = svdvals(curlmat)
print('Rank of curl =', sum(s>1e-14))

Rank of curl = 160

Of course, the rank should be at most the dimension of the codomain, which is immediately retrieved from the ndof attribute of the finite element space:

D.ndof

160

Thus we are computationally led to the interesting discovery, $$\newcommand{\rank}{\mathop{\mathrm{rank}}} \rank(\curl: N_{h, p}) = \dim P_p(\oh), $$ i.e., the 2D operator $\curl : N_{h, p} \to P_p(\oh)$ is surjective!

Please feel free to hit the above code with other meshes. You will end up with the same observation $\ldots$ so it might be a good idea to add this to the things you should rigorously understand, at least for the lowest order case: see Exercise III below.

Having discussed the range, let us now briefly turn to the kernel (or the null space) of the same linear operator $\curl: N_{h, p} \to P_p(\oh)$. Recalling the rank-nullity theorem from linear algebra, we can immediately compute the dimension of the kernel using the above-computed dimension of the range: $$ \begin{aligned} \dim(\ker(\curl:N_{h,p})) & = \dim(N_{h, p}) - \rank(\curl:N_{h, p}) \\ & = \dim(N_{h, p}) - \dim P_p(\oh) \end{aligned} $$

dimkercurl = N.ndof - D.ndof
dimkercurl

100

The emergence of topology

In calculus, we learned that every curl-free vector field can be written as a gradient of a scalar function, provided the domain has the topological property of being simply connected. However, we punched two holes in our domain above, so it is not simply connected. Hence we cannot assert that all curl-free functions on our domain are gradients. The dimension of those that are not turns out to be related to the toplogy of the domain.

The number of holes in a two dimensional domain is a topological invariant called the second Betti number $b_1$. (The first Betti number $b_0$ is the number of connected components of the domain, which is clearly 1 for our example.) A deep theorem of de Rham connects the topological Betti numbers to quantities arising from smooth function spaces (namely the dimensions of certain cohomology spaces of infinitely smooth functions).

We are now going to see that the number $b_1$ also arises naturally from the not-so-smooth (but computable) vector finite elements, as the dimension of curl-free functions that are not gradients within them. We have already computed the dimension of curl-free functions (stored above in dimkercurl variable). So it remains to discuss the dimension of the gradient fields in $N_{h,p}$. All gradients of functions in the Lagrange space $$ V_{h, p+1} = \{ v \in H^1(\om):\; \text{ for all } K \in \oh, \; v|_K \in P_{p+1}(K)\} $$ must have zero curl. Moreover, an elementary computation of the gradient and an investigation of the resulting function's tangential component along mesh interfaces proves that $$ \grad V_{h, p+1} \subset N_{h, p}. $$ (See Exercise IV at the end.) Thus the kernel of $\curl: N_{h, p} \to P_p(\oh)$ must contain at least all the gradient fields generated by differentiating Lagrange functions. The size of this gradient subspace, i.e., the dimension of $\grad V_{h, p+1}$, is easy to compute because we know that on a connected domain, the kernel of $\grad: V_{h, p+1}\to N_{h, p}$ consists of the one-dimensional space of constant functions. Hence, by the rank-nullity theorem, $$ \dim \grad V_{h, p+1} = \dim V_{h, p+1} - 1. $$

V = ng.H1(mesh, order=p+1)
rankgrad = V.ndof -1
rankgrad

98

We computed the exact dimension of the kernel of curl previously. The emergence of topology is evident when we compare the above number with the previously computed $\dim (\ker \curl)$. The number of holes in the domain is obtained as the difference of $\dim(\ker(\curl: N_{h,p}))$ and $\dim(\grad V_{h, p+1})$.

dimkercurl - rankgrad

2

Class Activity: Please return to the code above, remake the domain with as many holes as you like, change the meshsize to whatever you like, change the polynomial degree $p$, etc., and repeat these computations. (Also see exercises at the end.)

After you have experimented enough, your computational experience should lead you to build some confidence in formulating the hypothesis that the second Betti number must be given by $$\newcommand{\grad}{\mathop{\mathrm{grad}}} b_1 = \dim (\ker(\curl: N_{h, p}) - \dim(\grad V_{h, p+1}). $$

We will have more to say regarding this in later lectures. Note that if you replace the Nedelec space $N_{h, p}$ by the Cartesian product of Lagrange spaces in the right hand side above, the resulting number has no topological relevance (see Exercise V below).

To conclude, we have seen that the vector finite elements of the Nedelec and Raviart-Thomas type have special properties not shared by vector fields from simple Cartesian products of Lagrange spaces, including superconvergence of certain derivatives of interpolation errors, and connection with topological invariants. Hopefully these elementary examples, presented with very few prerequisites, serves as an invitation to theoretical studies into finite element exterior calculus. In a later lecture, we will connect simplicial complexes to finite element spaces and provide a theoretical foundation to base some of the above computational observations.

Exercises

Exercise I: Let $N_{h, p}$ and $V_{h, p}$ denote the Nedelec and Lagrange finite element spaces defined above. Let $I_h^N$ denote the canonical interpolant of Nedelec space and let $I_h$ denote the Lagrange interpolant on the Cartesian product of Lagrange spaces $V_{h, p+1} \times V_{h, p+1}$. Consider a smooth vector field, $u = (\sin (xy), y)$ on the unit square. Compute the following errors $$ \|u - I_h u\|_{L^2(\om)}, \quad \|\curl(u - I_h u)\|_{L^2(\om)}, \quad \|u - I_h^N u\|_{L^2(\om)}, \quad \|\curl(u - I_h^N u)\|_{L^2(\om)}, $$ on successive uniform refinements of a mesh. Tabulate rates and compare.

Exercise II: Consider an edge $e_{ij}$ between adjacent mesh vertices $a_i , a_j$ of a triangular mesh. Let $\lambda_i$ denote the barycentric coordinate function associated to the vertex $a_i$. Show that the function $$\newcommand{\grad}{\mathop{\mathrm{grad}}} \phi_{ij} = \lambda_j \grad \lambda_i - \lambda_i \grad \lambda_j $$

has zero tangential components on all mesh edges except $e_{ij}$, where it equals $ \phi_{ij} \cdot t_{ij} = 1/ |e_{ij}|$. Here $|e_{ij}|$ is the length of $e_{ij}$ and $t_{ij} = (a_i - a_j)/ |e_{ij}|$.

Exercise III: Prove that $\curl: N_{h, 0} \to P_0(\oh)$ is surjective on 2D domains. (You may want to start by computing the curl of the shape functions in Exercise 1.)

Exercise IV: Prove that $\grad V_{h, p+1} \subset N_{h, p}$.

Exercise V: We have seen that there is computational evidence that the second Betti number satisfies $$ b_1 = \dim (\ker(\curl: N_{h, p})) - \dim(\grad V_{h, p+1}). $$ Show, through computational examples, that this equality does not hold if we change the Nedelec space $N_{h, p}$ above to the Cartesian product of Lagrange spaces $V_{h, p+1} \times V_{h, p+1}$.

(Suggestion: Use ng.Grad on ng.VectorH1 to get derivatives.)

Exercise VI: Somewhere in this notebook we assumed connectedness of the domain. What happens if we remove this assumption? Understanding this, write a python function that takes as input a mesh and outputs both the first Betti number $b_0$ (the number of connected components of the domain) and the second Betti number $b_1$ (the number of holes in the domain). Display your results from meshes of connected and disconnected domains, with and without holes.

Exercise VII: Consider the Lagrange and Nedelec spaces with boundary conditions, namely $\mathring{V}_{h, p+1} = \{ v \in V_{h, p+1}: v|_{\partial \om } = 0\},$ and $\mathring{N}_{h, p} = \{ u \in N_{h, p}: u\cdot t|_{\partial \om } = 0\}.$ Determine, computationally or theoretically, if $\curl : \mathring{N}_{h,0} \to P_0(\oh)$ is surjective. Is $\grad \mathring{V}_{h, p+1}$ contained in $ \mathring{N}_{h, p}$? Can you compute $b_0$ and $b_1$ like in Exercise VI but using $\mathring{V}_{h, p+1}$ and $ \mathring{N}_{h, p}$ instead?

(Suggestion: Do not forget to use FreeDofs() when there are boundary conditions.)

---

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