This page was generated from unit-1.4-staticcond/staticcond.ipynb.

1.4 Static Condensation¶

Static condensation refers to the process of eliminating unknowns that are internal to elements from the global linear system. They are useful in standard methods and critical in methods like the HDG method. NGSolve automates static condensation across a variety of methods via a classification of degrees of freedom.

In this tutorial, you will learn

how to turn on static condensation (eliminate_internal or condense),
how to get the full solution from the condensed system using data members
- harmonic_extension_trans
- harmonic_extension
- inner_solve,
their relationship to a Schur complement factorization,
types of degrees of freedom such as COUPLING_TYPE.LOCAL_DOF, and
an automatic utility for solving via static condensation.

[1]:

import netgen.gui
from ngsolve import *
from netgen.geom2d import unit_square

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

fes = H1(mesh, order=4, dirichlet='bottom|right')
u, v = fes.TnT()

Asking `BilinearForm` to `condense`¶

Keyword arguments condense and eliminate_internal are synonymous.

[2]:

a = BilinearForm(fes, condense=True)
a += grad(u) * grad(v) * dx
a.Assemble()

f = LinearForm(fes)
f += 1 * v * dx
f.Assemble()

u = GridFunction(fes)

The assembled matrix \(A=\) a.mat can be block partitioned into

\[\begin{split}A = \left( \begin{array}{cc} A_{LL} & A_{LE} \\ A_{EL} & A_{EE} \end{array} \right)\end{split}\]

where \(L\) denotes the set of local (or internal) degrees of freedom and \(E\) denotes the set of interface degrees of freedom.

In our current example \(E\) consists of edge and vertex dofs while \(L\) consists of triangle dofs. (Note that in practice, \(L\) and \(E\) may not be ordered contiguously and \(L\) need not appear before \(E\), but such details are immaterial for our discussion here.)
The condensed system is known as the Schur complement:

\[S = A_{EE} - A_{EL} A_{LL}^{-1} A_{LE}.\]

When eliminate_internal (or condense) is set to True in a, the statement a.Assemble actually assembles \(S\).

A factorization of the inverse¶

NGSolve provides
- a.harmonic_extension_trans \(= \left(\begin{array}{cc} 0 & 0 \\ -A_{EL} A_{LL}^{-1} & 0 \end{array}\right)\)
- a.harmonic_extension \(= \left(\begin{array}{cc} 0 & -A_{LL}^{-1} A_{LE} \\ 0 & 0 \end{array}\right)\)
- a.inner_solve \(=\left(\begin{array}{cc} A_{LL}^{-1} & 0 \\ 0 & 0 \end{array}\right)\).
To solve

\[\begin{split} \left(\begin{array}{cc} A_{LL} & A_{LE}\\ A_{EL} & A_{EE} \end{array}\right) \left(\begin{array}{c} u_L \\ u_E \end{array}\right)= \left(\begin{array}{c} f_L \\ f_E \end{array}\right)\end{split}\]

we use a factorization of \(A^{-1}\) that uses \(S^{-1}\). Namely, we use the following identity:

\[\begin{split}\left(\begin{array}{c} u_L \\ u_E \end{array}\right) = \left(\begin{array}{cc} I & -A_{LL}^{-1} A_{LE} \\ 0 & I \end{array}\right) \left(\begin{array}{cc} A_{LL}^{-1} & 0 \\ 0 & S^{-1} \end{array}\right) \underbrace{ \left(\begin{array}{cc} I & 0 \\ -A_{EL} A_{LL}^{-1} & I \end{array}\right) \left(\begin{array}{c} f_L \\ f_E \end{array}\right)}_{ \left(\begin{array}{c} f'_L \\ f'_E \end{array}\right) }\end{split}\]

We implement this formula step by step, starting with the computation of \(f_L'\) and \(f_E'\).

Full solution from the condensed system¶

The following step implements

\[\begin{split}\left(\begin{array}{c} f'_L \\ f'_E \end{array}\right) = \left(\begin{array}{cc} I & 0 \\ -A_{EL} A_{LL}^{-1} & I \end{array}\right) \left(\begin{array}{c} f_L \\ f_E \end{array}\right).\end{split}\]

[3]:

f.vec.data += a.harmonic_extension_trans * f.vec

The next step implements part of the next matrix application in the formula.

\[\begin{split}\left(\begin{array}{c} 0 \\ u_E \end{array}\right) = \left(\begin{array}{cc} 0 & 0 \\ 0 & S^{-1} \end{array}\right) \left(\begin{array}{c} f'_L \\ f'_E \end{array}\right).\end{split}\]

[4]:

u.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs(coupling=True)) * f.vec

Note:
- Because a has condense set, the inverse a.mat.Inverse actually computes \(S^{-1}\).
- Note that instead of the usual fes.FreeDofs(), we have used fes.FreeDofs(coupling=True) or simply fes.FreeDofs(True) to specify that only the degrees of freedom that are not local and not Dirichlet should participate in the inverse computations. (The underlying assumption is that Dirichlet dofs cannot be local dofs.)
Next, we compute

\[\begin{split}\left(\begin{array}{c} u'_L \\ u_E \end{array}\right) = \left(\begin{array}{c} 0 \\ u_E \end{array}\right) + \left(\begin{array}{cc} A_{LL}^{-1} & 0 \\ 0 & 0 \end{array}\right) \left(\begin{array}{c} f'_L \\ f'_E \end{array}\right).\end{split}\]

[5]:

u.vec.data += a.inner_solve * f.vec

Finally:

\[\begin{split}\left(\begin{array}{c} u_L \\ u_E \end{array}\right) = \left(\begin{array}{cc} I & -A_{LL}^{-1} A_{LE} \\ 0 & I \end{array}\right) \left(\begin{array}{c} u_L' \\ u_E \end{array}\right)\end{split}\]

[6]:

u.vec.data += a.harmonic_extension * u.vec
Draw(u)

Behind the scenes: `CouplingType`¶

How does NGSolve know what is in the index sets \(L\) and \(E\)?

Look at fes.CouplingType to see a classification of degrees of freedom (dofs).

[7]:

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

COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.INTERFACE_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF
COUPLING_TYPE.LOCAL_DOF

[8]:

dof_types = {}
for i in range(fes.ndof):
    ctype = fes.CouplingType(i)
    if ctype in dof_types.keys():
        dof_types[ctype] += 1
    else:
        dof_types[ctype] = 1
dof_types

[8]:

{COUPLING_TYPE.WIREBASKET_DOF: 12,
 COUPLING_TYPE.INTERFACE_DOF: 75,
 COUPLING_TYPE.LOCAL_DOF: 42}

The LOCAL_DOF forms the set \(L\) and the remainder forms the set \(E\). All finite element spaces in NGSolve have such dof classification.

Through this classification a bilinear form is able to automatically compute the Schur complement and the accompanying extension operators. Users need only specify the eliminate_internal option. (Of course users should also make sure their method has an invertible \(A_{LL}\)!)

Inhomogeneous Dirichlet boundary conditions¶

In case of inhomogeneous Dirichlet boundary conditions, we combine the technique of Dirichlet data extension in 1.3 with the above static condensation principle in the following code.

[9]:

U = x*x*(1-y)*(1-y)          # U = manufactured solution
DeltaU = 2*((1-y)*(1-y)+x*x) # Source: DeltaU = ∆U

f = LinearForm(fes)
f += -DeltaU * v * dx
f.Assemble()

u = GridFunction(fes)
u.Set(U, BND)               # Dirichlet b.c: u = U on boundary

r = f.vec.CreateVector()
r.data = f.vec - a.mat * u.vec
r.data += a.harmonic_extension_trans * r

u.vec.data += a.mat.Inverse(fes.FreeDofs(True)) * r
u.vec.data += a.harmonic_extension * u.vec
u.vec.data += a.inner_solve * r

sqrt(Integrate((U-u)*(U-u),mesh))  # Compute error

[9]:

7.560675945217909e-16

If you find the above steps onerous, you can again use the automatic utility solvers.BVP, which will perform static condensation if condense is turned on in its input bilinear form.

Automatic utility `BVP`¶

[10]:

u = GridFunction(fes)
u.Set(U, BND)

c = Preconditioner(a,"direct")
c.Update()
solvers.BVP(bf=a, lf=f, gf=u, pre=c)
sqrt(Integrate((U-u)*(U-u),mesh))

iteration 0 error = 0.4613481630629836
iteration 1 error = 9.396751530627324e-16

[10]:

7.725365972187299e-16

1.4 Static Condensation¶

Asking BilinearForm to condense¶

A factorization of the inverse¶

Full solution from the condensed system¶

Behind the scenes: CouplingType¶

Inhomogeneous Dirichlet boundary conditions¶

Automatic utility BVP¶

Asking `BilinearForm` to `condense`¶

Behind the scenes: `CouplingType`¶

Automatic utility `BVP`¶