This page was generated from unit-2.3-hcurlhdiv/hcurlhdiv.ipynb.

2.3 H(curl) and H(div) function spaces¶

Scalar and vectorial finite elements in NGSolve:

Standard continuous \(H^1\) elements: title

Nedelec’s tangentially-continuous \(H(curl)\)-conforming edge elements:

Raviart-Thomas normally-continuous \(H(div)\)-conforming face elements:

Discontinuous \(L_2\) elements:

These vector-valued spaces allow to represent physical quantities which are either normally or tangentially continuous.

The finite element spaces are related by the de Rham complex:

\[\begin{split}\DeclareMathOperator{\Grad}{grad} \DeclareMathOperator{\Curl}{curl} \DeclareMathOperator{\Div}{div} \begin{array}{ccccccc} H^1 & \stackrel{\Grad}{\longrightarrow} & H(\Curl) & \stackrel{\Curl}{\longrightarrow} & H(\Div) & \stackrel{\Div}{\longrightarrow} & L^2 \\[8pt] \bigcup & & \bigcup & & \bigcup & & \bigcup \\[8pt] W_{h} & \stackrel{\Grad}{\longrightarrow} & V_{h } & \stackrel{\Curl}{\longrightarrow} & Q_{h} & \stackrel{\Div}{\longrightarrow} & S_{h} \: \\[3ex] \end{array}\end{split}\]

NGSolve supports these elements of arbitrary order, on all common element shapes (trigs, quads, tets, prisms, pyramids, hexes). Elements may be curved.

[1]:

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

Generate a higher order \(H^1\)-space. We first explore its different types of basis functions.

[2]:

order=3
fes = H1(mesh, order=order)
gfu = GridFunction(fes)

The first basis functions are hat-functions, one per vertex. By setting the solution vector to a unit-vector, we may look at the individual basis functions:

[3]:

gfu.vec[:] = 0
# vertex nr 17:
gfu.vec[17] = 1
Draw(gfu, min=0, max=1, deformation=True);

The next are edge-bubbles, where we have \((p-1)\) basis functions per edge. A NodeId object refers to a particular vertex, edge, face or cell node in the mesh. We can ask for the degrees of freedom on a node:

[4]:

# basis functions on edge nr:
edge_dofs = fes.GetDofNrs(NodeId(EDGE,10))
print("edge_dofs =", edge_dofs)
gfu.vec[:] = 0
gfu.vec[edge_dofs[0]] = -1
Draw(gfu, order=3, min=-0.05, max=0.05, deformation=True);

edge_dofs = (39, 40)

Finally, we have \((p-1)(p-2)/2\) inner basis functions on every triangle:

[5]:

trig_dofs = fes.GetDofNrs(NodeId(FACE,0))
print("trig_dofs = ", trig_dofs)
gfu.vec[:] = 0
gfu.vec[trig_dofs[0]] = 10
Draw(gfu, order=3, min=0, max=0.3, deformation=True);

trig_dofs =  (103,)

The FESpace also maintains information about local dofs, interface dofs and wire-basket dofs for the BDDC preconditioner:

[6]:

for i in range(fes.ndof):
    print (i,":", fes.CouplingType(i))

0 : COUPLING_TYPE.WIREBASKET_DOF
1 : COUPLING_TYPE.WIREBASKET_DOF
2 : COUPLING_TYPE.WIREBASKET_DOF
3 : COUPLING_TYPE.WIREBASKET_DOF
4 : COUPLING_TYPE.WIREBASKET_DOF
5 : COUPLING_TYPE.WIREBASKET_DOF
6 : COUPLING_TYPE.WIREBASKET_DOF
7 : COUPLING_TYPE.WIREBASKET_DOF
8 : COUPLING_TYPE.WIREBASKET_DOF
9 : COUPLING_TYPE.WIREBASKET_DOF
10 : COUPLING_TYPE.WIREBASKET_DOF
11 : COUPLING_TYPE.WIREBASKET_DOF
12 : COUPLING_TYPE.WIREBASKET_DOF
13 : COUPLING_TYPE.WIREBASKET_DOF
14 : COUPLING_TYPE.WIREBASKET_DOF
15 : COUPLING_TYPE.WIREBASKET_DOF
16 : COUPLING_TYPE.WIREBASKET_DOF
17 : COUPLING_TYPE.WIREBASKET_DOF
18 : COUPLING_TYPE.WIREBASKET_DOF
19 : COUPLING_TYPE.INTERFACE_DOF
20 : COUPLING_TYPE.INTERFACE_DOF
21 : COUPLING_TYPE.INTERFACE_DOF
22 : COUPLING_TYPE.INTERFACE_DOF
23 : COUPLING_TYPE.INTERFACE_DOF
24 : COUPLING_TYPE.INTERFACE_DOF
25 : COUPLING_TYPE.INTERFACE_DOF
26 : COUPLING_TYPE.INTERFACE_DOF
27 : COUPLING_TYPE.INTERFACE_DOF
28 : COUPLING_TYPE.INTERFACE_DOF
29 : COUPLING_TYPE.INTERFACE_DOF
30 : COUPLING_TYPE.INTERFACE_DOF
31 : COUPLING_TYPE.INTERFACE_DOF
32 : COUPLING_TYPE.INTERFACE_DOF
33 : COUPLING_TYPE.INTERFACE_DOF
34 : COUPLING_TYPE.INTERFACE_DOF
35 : COUPLING_TYPE.INTERFACE_DOF
36 : COUPLING_TYPE.INTERFACE_DOF
37 : COUPLING_TYPE.INTERFACE_DOF
38 : COUPLING_TYPE.INTERFACE_DOF
39 : COUPLING_TYPE.INTERFACE_DOF
40 : COUPLING_TYPE.INTERFACE_DOF
41 : COUPLING_TYPE.INTERFACE_DOF
42 : COUPLING_TYPE.INTERFACE_DOF
43 : COUPLING_TYPE.INTERFACE_DOF
44 : COUPLING_TYPE.INTERFACE_DOF
45 : COUPLING_TYPE.INTERFACE_DOF
46 : COUPLING_TYPE.INTERFACE_DOF
47 : COUPLING_TYPE.INTERFACE_DOF
48 : COUPLING_TYPE.INTERFACE_DOF
49 : COUPLING_TYPE.INTERFACE_DOF
50 : COUPLING_TYPE.INTERFACE_DOF
51 : COUPLING_TYPE.INTERFACE_DOF
52 : COUPLING_TYPE.INTERFACE_DOF
53 : COUPLING_TYPE.INTERFACE_DOF
54 : COUPLING_TYPE.INTERFACE_DOF
55 : COUPLING_TYPE.INTERFACE_DOF
56 : COUPLING_TYPE.INTERFACE_DOF
57 : COUPLING_TYPE.INTERFACE_DOF
58 : COUPLING_TYPE.INTERFACE_DOF
59 : COUPLING_TYPE.INTERFACE_DOF
60 : COUPLING_TYPE.INTERFACE_DOF
61 : COUPLING_TYPE.INTERFACE_DOF
62 : COUPLING_TYPE.INTERFACE_DOF
63 : COUPLING_TYPE.INTERFACE_DOF
64 : COUPLING_TYPE.INTERFACE_DOF
65 : COUPLING_TYPE.INTERFACE_DOF
66 : COUPLING_TYPE.INTERFACE_DOF
67 : COUPLING_TYPE.INTERFACE_DOF
68 : COUPLING_TYPE.INTERFACE_DOF
69 : COUPLING_TYPE.INTERFACE_DOF
70 : COUPLING_TYPE.INTERFACE_DOF
71 : COUPLING_TYPE.INTERFACE_DOF
72 : COUPLING_TYPE.INTERFACE_DOF
73 : COUPLING_TYPE.INTERFACE_DOF
74 : COUPLING_TYPE.INTERFACE_DOF
75 : COUPLING_TYPE.INTERFACE_DOF
76 : COUPLING_TYPE.INTERFACE_DOF
77 : COUPLING_TYPE.INTERFACE_DOF
78 : COUPLING_TYPE.INTERFACE_DOF
79 : COUPLING_TYPE.INTERFACE_DOF
80 : COUPLING_TYPE.INTERFACE_DOF
81 : COUPLING_TYPE.INTERFACE_DOF
82 : COUPLING_TYPE.INTERFACE_DOF
83 : COUPLING_TYPE.INTERFACE_DOF
84 : COUPLING_TYPE.INTERFACE_DOF
85 : COUPLING_TYPE.INTERFACE_DOF
86 : COUPLING_TYPE.INTERFACE_DOF
87 : COUPLING_TYPE.INTERFACE_DOF
88 : COUPLING_TYPE.INTERFACE_DOF
89 : COUPLING_TYPE.INTERFACE_DOF
90 : COUPLING_TYPE.INTERFACE_DOF
91 : COUPLING_TYPE.INTERFACE_DOF
92 : COUPLING_TYPE.INTERFACE_DOF
93 : COUPLING_TYPE.INTERFACE_DOF
94 : COUPLING_TYPE.INTERFACE_DOF
95 : COUPLING_TYPE.INTERFACE_DOF
96 : COUPLING_TYPE.INTERFACE_DOF
97 : COUPLING_TYPE.INTERFACE_DOF
98 : COUPLING_TYPE.INTERFACE_DOF
99 : COUPLING_TYPE.INTERFACE_DOF
100 : COUPLING_TYPE.INTERFACE_DOF
101 : COUPLING_TYPE.INTERFACE_DOF
102 : COUPLING_TYPE.INTERFACE_DOF
103 : COUPLING_TYPE.LOCAL_DOF
104 : COUPLING_TYPE.LOCAL_DOF
105 : COUPLING_TYPE.LOCAL_DOF
106 : COUPLING_TYPE.LOCAL_DOF
107 : COUPLING_TYPE.LOCAL_DOF
108 : COUPLING_TYPE.LOCAL_DOF
109 : COUPLING_TYPE.LOCAL_DOF
110 : COUPLING_TYPE.LOCAL_DOF
111 : COUPLING_TYPE.LOCAL_DOF
112 : COUPLING_TYPE.LOCAL_DOF
113 : COUPLING_TYPE.LOCAL_DOF
114 : COUPLING_TYPE.LOCAL_DOF
115 : COUPLING_TYPE.LOCAL_DOF
116 : COUPLING_TYPE.LOCAL_DOF
117 : COUPLING_TYPE.LOCAL_DOF
118 : COUPLING_TYPE.LOCAL_DOF
119 : COUPLING_TYPE.LOCAL_DOF
120 : COUPLING_TYPE.LOCAL_DOF
121 : COUPLING_TYPE.LOCAL_DOF
122 : COUPLING_TYPE.LOCAL_DOF
123 : COUPLING_TYPE.LOCAL_DOF
124 : COUPLING_TYPE.LOCAL_DOF
125 : COUPLING_TYPE.LOCAL_DOF
126 : COUPLING_TYPE.LOCAL_DOF

\(H(curl)\) finite element space¶

In NGSolve we use hierarchical high order finite element basis functions with node-wise exact sequences. The lowest order space \(W_{l.o}\) is the edge-element space:

\[\begin{split}\begin{array}{rcll} W_{hp} & = & W_{p=1} + \sum_E W_E + \sum_F W_F + \sum_C W_C & \subset H^1 \\[0.5em] V_{hp} & = & W_{l.o} + \sum_E V_E + \sum_F V_F + \sum_C V_C & \subset H(curl) \end{array}\end{split}\]

where the edge, face and cell blocks are compatible in the sense that

\[\nabla W_E = V_E, \quad \nabla W_F \subset V_F, \quad \nabla W_C \subset V_C\]

We obtain this by using gradients of \(H^1\) basis functions as \(H(curl)\) basis functions, and some more (see thesis Sabine Zaglmayr):

\[\begin{split}\begin{array}{rcl} V_E & = & \text{span} \{ \nabla \varphi_{E,i}^{H^1} \} \\ V_F & = & \text{span} \{ \nabla \varphi_{F,i}^{H^1} \cup \widetilde \varphi_{F,i}^{H(curl)} \} \\ V_C & = & \text{span} \{ \nabla \varphi_{C,i}^{H^1} \cup \widetilde \varphi_{C,i}^{H(curl)} \} \end{array}\end{split}\]

[7]:

fes = HCurl(mesh, order=2)
uc = GridFunction(fes, name="uc")

[8]:

edge_dofs = fes.GetDofNrs(NodeId(EDGE,10))
print ("edgedofs: ", edge_dofs)
uc.vec[:] = 0
uc.vec[edge_dofs[0]] = 1
Draw (uc, min=0, max=3, vectors = { "grid_size":30})
Draw (curl(uc), mesh, "curl", min=-25, max=25);

edgedofs:  (10, 62, 63)

[9]:

face_dofs = fes.GetDofNrs(NodeId(FACE,10))
print ("facedofs: ", face_dofs)
uc.vec[:] = 0
uc.vec[face_dofs[0]] = 1
Draw (uc, min=0, max=1, vectors = { "grid_size":30})
Draw (curl(uc), mesh, "curl", min=-1, max=1, order=3); # it's a gradient

facedofs:  (135, 136, 137)

\(H(div)\) finite element space¶

NGSolve provides Raviart-Thomas (RT) as well as Brezzi-Douglas-Marini (BDM) finite element spaces for H(div). We obtain the RT-version by setting RT=True, otherwise we get BDM.

[10]:

fes = HDiv(mesh, order=2, RT=True)
ud = GridFunction(fes)
func = x*y*(x,y)
ud.Set (func)
Draw (ud, vectors = { "grid_size":30})
print ("interpolation error:", Integrate ((func-ud)**2, mesh))

interpolation error: 1.3401439357345005e-31

The function spaces know their canonical derivatives. These operations are efficiently implemented by transformation from the reference element.

[11]:

gfu.derivname, ud.derivname, uc.derivname

[11]:

('grad', 'div', 'curl')

But there are additional options, like forming the element-wise gradient of H(div) finite element functions. We can query the available operators via

[12]:

print ("H(div) operators: ", ud.Operators())

H(div) operators:  ['grad', 'dual', 'normalcomponent']

and access them via the Operator() method

[13]:

Draw (grad(ud)[0,1], mesh, "gradud");

[ ]: