Vector elements in 2D

4. Vector elements in 2D#

The finite element spaces we have seen so far, such as the Lagrange, the DG, and the nonconforming finite element spaces, are all spaces of scalar-valued functions. How can we approximate vector fields? We could certainly use a Cartesian product of scalar finite element spaces to approximate a vector field component by component. However, better properties can often be obtained using tailored vector finite element spaces, such as the Raviart-Thomas (RT) space or the Nedelec space. In this notebook we see how to access these spaces in ngsolve and compute with them.

4.1. The Raviart-Thomas space#

\(\newcommand{\om}{\varOmega} \newcommand{\oh}{\varOmega_h} \newcommand{\dive}{\mathop{\mathrm{div}}} \newcommand{\grad}{\mathop{\mathrm{grad}}} \renewcommand{\div}{\mathop{\mathrm{div}}} \newcommand{\R}{\mathbb{R}} \newcommand{\jump}{\mathop{\mathrm{jump}}} \) Consider a facet \(F\) shared by two adjacent elements \(K^\pm\) in two-dimensional mesh \(\oh\) of a domain \(\om\). In a previous notebook, we examined jumps of a scalar piecewise function \(w\) on the mesh, written as ⟦ \( w\) ⟧ or simply \(\jump(w)\). Here, we consider jumps of the normal component of a piecewise smooth vector field \(q: \om \to \R^2\) across the facet \(F\). Define ⟦ \(q\cdot n\) ⟧, also denoted by \(\jump(q\cdot n)\), by

\[\tag{1} \jump(q \cdot n)|_F := \Big(q|_{K^+} \cdot n_+ + q|_{K^-} \cdot n_-\Big) \Big|_F, \]

where \(n^\pm\) denote the outward unit normal on the boundary of the element \(K^\pm\). Since the sign of \(n_+\) and \(n_-\) are opposite on \(F\), this indeed measures the discontinuity in the normal component of \(q\). This definition, facet by facet, defines the normal jump function \(\jump(q \cdot n)\) of the vector field \(q\) on the union of all interior facets. (Note that this jump definition if independent of the facet orientation.)

The Raviart-Thomas (RT) finite element space on a two-dimensional mesh \(\newcommand{\oh}{\varOmega_h}\) \(\oh\) (of a domain \(\newcommand{\om}{\varOmega} \om\)) is defined by

\[\begin{split} \begin{aligned} R_{hp} = \Big\{ q : \; \newcommand{\jump}{\mathop{\mathrm{jump}}} & q|_K \in P_p(K)^2 + \begin{pmatrix} x \\ y \end{pmatrix} P_p(K) \; \text{ on all mesh elements } K, \\ &\text{ and } \jump(q\cdot n) =0 \text{ on all interior mesh facets} \Big\} \end{aligned} \end{split}\]

for any \(p \ge 0\). Here, \(P_p(K)\) denote the set of polynomials of degree at most \(p\). From the definition, it is clear that plots of any function in the RT space should show a possibly discontinuous vector field with 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. Fluxes, such as the flux of a fluid flow, current flow, or magnetic flux, are examples of vector fields whose normal components remain continuous across a material interface, even when its other components jump across the interface.

4.2. Shape functions#

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

Let us examine the (local) element shape functions of the lowest order RT space on a single triangle. We first create a mesh of a single triangle.

wp = WorkPlane() # Make a single triangle
wp.MoveTo(0,0).Line(1).Rotate(45+90).Line(ng.sqrt(2)).Close()
singletrg = wp.Face()
trg = ng.Mesh(OCCGeometry(singletrg, dim=2).GenerateMesh(maxh=100))

The RT space on this triangle or any mesh is created using ng.HDiv as follows. The lowest order case is obtained by setting order=0 and RT=True.

R = ng.HDiv(trg, order=0, RT=True)
R.ndof # Number of degrees of freedom

3

This shows that there are three degrees of freedom in this cas. Recall that basis functions dual to the set of degrees of freedom are called shape functions. Since we cannot visualize the degrees of freedom functionals, visualizing shape functions is one way to understand what degrees of freedom are implemented in a code.

To visualize the shape functions, as usual, we create a GridFunction in the RT space and set its vector to be zero except for one entry corresponding to the desired shape function. Plotting the resulting vector field then shows the desired shape function. Here are the shape functions for all three degrees of freedom.

for i in range(R.ndof):
    shapenumber = i
    shape = ng.GridFunction(R, name='shape')
    shape.vec[:] = 0
    shape.vec[shapenumber] = 1
    Draw(shape, vectors={'grid_size': 20});

It is clear from these plots that the normal component of each shape function is zero on the boundary of the triangle, except for one edge.

The shape functions of the global RT space are obtained by fusing together such local RT shape functions on each mesh element in such a manner as to obtain continuity of the normal component across the facet. Plotting the global shape functions (by the same procedure as above) shows this normal component continuity. Here is how a global shape function which has zero normal component on all mesh facets except one facet looks like.

mesh1 = ng.Mesh(ng.unit_square.GenerateMesh(maxh=0.5))
R = ng.HDiv(mesh1, order=0, RT=True)

shapenumber = R.ndof-1
shape = ng.GridFunction(R, name='shape')
shape.vec[:] = 0
shape.vec[shapenumber] = 1
Draw(shape, 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. Every edge is thus associated to a global degree of freedom of the RT space.

Note the number of mesh facets and the dimension of the lowest order RT space by querying both numbers:

mesh1.nfacet   # or nedge

13

R.ndof

13

As expected, these two numbers are equal.

4.2.1. Approximation#

Consider approximating the following vector field using the RT space:

\[\begin{split} q = \left\{ \begin{aligned} & (\sin(x), y), && x \le 1/2, \\ & (\sin(x), 1-y), && x > 1/2. \end{aligned} \right. \end{split}\]

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.

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))

We make the needed piecewise constant coefficient function using the MaterialCF method in ngsolve.

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. (We have already seen a similar phenomena around discontinuities of scalar function approximation.) Instead, knowing that the normal component is continuous, we can try to approximate the function into the RT space. Such an approximation liberates its tangential component from continuity constraints, while preserving the normal continuity.

The Set method, as in other finite element spaces in NGSolve, can be used to perform the Oswald approximation in the RT space. The procedure for Oswald approximation, already described previously, is the same for all finite element spaces: one first \(L_2\)-projects (component by component) the (vector) function to be approximated into each local finite element space, and then assigns to every coupling degree of freedom the average of those degrees of freedom of the projection that must be fused to form that coupling degree of freedom. In particular, in the RT case, the Oswald interpolation does not average the tangential components across the facets (since those are not coupling degrees of freedom in the RT space) but only the normal components.

The \(p\) in the next code cell indicates the polynomial degree in the definition of \(R_{hp}\) given earlier.

p = 5   # This is the p in the definition of R_hp 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.4319667352337606e-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.11081020823065106

Clearly the error in the product space approach is much 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, as seen from the plot of the error below.

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

4.2.2. Superconvergence of divergence#

Ordinarily, when a piecewise polynomial approximation of a function \(q\) converges in \(L_2\) at some rate, we expect its derivatives to converge at a lower rate. We already saw this rate reduction in a prior notebook.

However, somewhat miraculously, the Raviart-Thomas finite element interpolant \(I_h^{RT} q\) has the property that \(\newcommand{\dive}{\mathrm{div}} \dive(I_h^{RT} q - q)\) converges at the same rate as \(I_h^{RT} q - q\). All phenomena where convergence occurs at a rate higher than expected are called superconvergence.

To see the superconvergence of the divergence of RT interpolation error, we compute the interpolant using the Set method with dual=True option. We then compute the error measures

\[ \left( \sum_{K} \int_K |q - I_h^{RT} q|^2 \, dx \right)^{1/2} \quad \text{ and } \quad \left(\sum_{K} \int_K \dive(q - I_h^{RT} q)^2 \, dx \right)^{1/2} \]

on successively refined meshes and tabulate the convergence rates much the same way as in a previous notebook.

from prettytable import PrettyTable

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)

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\). We can compute the piecewise divergence of our simple q “by hand” (or using Diff on each piece). It is set below.

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.079646309574 |  *   |
|  0.25/2  |  0.039825655461 | 1.00 |
|  0.25/4  |  0.019913141104 | 1.00 |
|  0.25/8  |  0.009956609748 | 1.00 |
| 0.25/16  |  0.004978309774 | 1.00 |
| 0.25/32  |  0.002489155500 | 1.00 |
| 0.25/64  |  0.001244577826 | 1.00 |
| 0.25/128 |  0.000622288923 | 1.00 |
+----------+-----------------+------+
+----------+----------------------+------+
|    h     | ||div(q - IhRT(q))|| | rate |
+----------+----------------------+------+
|   0.25   |    0.030651751057    |  *   |
|  0.25/2  |    0.015370264479    | 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 |
+----------+----------------------+------+

Questions for discussion:

  • Do all the other derivatives (in the full gradient) superconverge as well?

  • Does the divergence superconverge if you use the Oswald approximation instead of the canonical RT interpolant?

  • What happens for higher degrees \(p\)?

4.3. The Nédélec finite element space#

Consider the curl operator in two space dimensions, acting on a vector field \(u \equiv \begin{bmatrix}u_0 \\ u_1 \end{bmatrix}: \varOmega \to \mathbb R^2\),

\[ \newcommand{\curl}{\mathop{\mathrm{curl}}} \curl u = \frac{\partial u_1}{\partial x} - \frac{\partial u_0}{\partial y}. \]

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 function \(\curl u\).

This 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.,

\[\begin{split} J_{\pi/2} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} b \\ -a \end{pmatrix}. \end{split}\]

Now, if two vector fields \(u\) and \(q\) are related by

\[ u = J_{\pi/2} q \]

then obviously

\[ \newcommand{\dive}{\mathop{\mathrm{div}}} \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.\) The definition of the tangential jump ⟦ \(u\cdot t\) ⟧ on an interface \(F = \partial K_+ \cap \partial K_-\) is completely analogous to how we defined the normal jump and just involves replacing \(n_\pm\) in that definition (1) with \(t_\pm\).

We define the Nédélec space in two dimensions \(N_{hp}\) as the rotated RT space, i.e.,

\[ N_{hp} := J_{\pi/2} R_{hp}. \]

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_{hp}\), we also find that \(N_{hp}\) can be equivalently described by

\[\begin{split} \newcommand{\curl}{\mathop{\mathrm{curl}}}\newcommand{\om}{\varOmega} \begin{aligned} N_{hp} = \Big\{ u : \; \newcommand{\jump}{\mathop{\mathrm{jump}}} & u|_K \in P_p(K)^2 + \begin{pmatrix} \phantom{-}y \\ -x \end{pmatrix} P_p(K) \text{ on all elements } K, \\ &\text{ and } \jump(u\cdot t) =0 \text{ on all interior mesh facets} \Big\}. \end{aligned} \end{split}\]

Here, the vector \((y, -x)\) arises as the rotation of the vector \((x, y)\) used in the RT space definition.

Note that although the 2D Nédélec space is isomorphic to the 2D Raviart-Thomas space (via \(J_{\pi/2}\)), the 3D Nédélec space is truly different from the 3D Raviart-Thomas space. (You can be convinced of this by counting their dimensions.)

NGSolve provides an implementation of the Nedelec space in 2D and 3D accessible through HCurl as follows. We provide the argument type1 because there are two spaces named Nédélec and the one we denoted by \(N_{hp}\) above is often referred to as the Nédélec space of type 1.

N = ng.HCurl(mesh, order=0, type1=True)

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

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

Questions for discussion:

  • Is there an RT shape function that looks like the above but rotated 90 degrees?

  • To approximate a piecewise vector field with tangential continuity, would you use the RT space or the Nédélec space?

  • Is there a derivative of the Nédélec interpolant that superconverges?

The Nédélec finite element is sometimes called the edge element in some engineering literature. It is an essential prerequisite for advanced electromagnetic simulations, as we shall see later. 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: what we have plotted above are indeed some scalar multiples of Whitney forms. Their expressions in terms of barycentric coordinates are left as exercises.

4.4. Summary#

We have seen

  • the 2D RT space with normal continuity,

  • the 2D Nédélec space with tangential continuity,

  • that in 2D, one is obtained by a rotation of the other,

  • the scalar 2D curl operator and its relation to divergence,

  • how to interpolate vector fields with normal or tangential continuity using these spaces,

  • the shape functions of RT and Nédélec spaces,

  • the superconvergence of divergence of RT interpolant, and

  • when not to expect such superconvergence for RT and Nédélec approximations.

Contents