Best approximations in 2D

2. Best approximations in 2D#

The ability of a finite element space to approximate functions outside it is foundational to its utility. Best approximations provide a quantitative measure of this ability.

Having seen finite elements on one-dimensional domains previously, we now consider two-dimensional domains and three finite element spaces on them, denoted by \(W_{hp}, V_{up}\) and \(U_{h}\), defined as follows:

  • Perhaps the simplest finite element space is \(W_{hp} = \{ w: \; w|_K\) is a polynomial of degree \( \le p\) on each mesh element \(K\},\) often called the DG space (or the “discontinuous Galerkin” space).

  • With more inter-element continuity restrictions, one obtains the smaller often-used space, the Lagrange space (whose lowest order case is also often called the Courant space), namely \(V_{hp} = \{ v: \; v \) is continuous and \(v|_K\) is a polynomial of degree \(\le p\) on each mesh element \( K\}.\)

  • The CR space (also known as the scalar Crouzeix-Raviart space or the \(P_1\)-nonconforming space) is \(U_h = \{ u: u|_K \) is linear on each mesh element \(K\) and is continuous at each mesh edge midpoint\(\}.\)

In this notebook, we examine these spaces and minimize a distance functional to compute the closest approximation of a function from these spaces. Such best approximations tell us how good each space approximates a given external function.

import ngsolve as ng
from ngsolve import Mesh, H1, L2, FESpace, GridFunction, x, y, dx, exp 
from netgen.occ import unit_square
from ngsolve.webgui import Draw

2.1. Shape functions in the lowest order case#

Let us begin by examining the lowest order cases of the three above-mentioned spaces on a coarsely meshed square domain.

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

The Draw facility above is used to produce an interactive plot of the mesh that you can rotate, pan, and zoom. It is also used to visualize GridFunctions (as seen below) and CoefficientFunction objects.

The DG space is obtained in ngsolve using L2 as seen below. For linear or higher degrees, the vector space given by L2 and Discontinuous(H1(...)) (seen previously) are the same. The lowest order DG space consists of piecewise constants:

W = ng.L2(mesh, order=0)
w = GridFunction(W)
w.vec[:] = 0; w.vec[W.ndof//2] = 0.2 
Draw(w, deformation=True, euler_angles=[-70,0,0]);

Clearly there is one shape function per mesh triangle that forms a basis for this space.

Next, consider the lowest order Lagrange space, whose functions, unlike the DG space, are continuous. The lowest order that gives a nontrivial space now is the case order=1. Note how that there is one shape function per vertex.

V = ng.H1(mesh, order=1)
v = GridFunction(V)
v.vec[:] = 0; v.vec[mesh.nv-1] = 0.2 
Draw(v, deformation=True, euler_angles=[-50,0,0]);

This is the familiar “hat function” associated to a vertex.

The third space is the CR space made below. It has one shape function associated to each mesh edge. The continuity at the edge midpoint is evident from the plot of the shape functions:

U = ng.FESpace('nonconforming', mesh) # only degree=1 available
u = GridFunction(U)
u.vec[:] = 0; u.vec[mesh.nedge-1] = 0.2 
Draw(u, deformation=True, euler_angles=[-60,0,-50]);

Question for discussion: For the Lagrange and CR spaces, how do your answers for the following differ?

  • What are the element spaces?

  • What are the element degrees of freedom?

  • How are element degrees of freedom fused together to form global degrees of freedom?

2.2. Higher order cases and coupling dofs#

Higher degree versions of the Lagrange and DG spaces are routinely used. Higher degree CR-type spaces are typically not used and are not implemented in ngsolve, but one can constrain the DG space to effectively work in such spaces (as we will see in a later section below).

The quadratic Lagrange space is made next and one of its shape functions is plotted below:

V2 = ng.H1(mesh, order=2)
v = GridFunction(V2)
v.vec[:] = 0;
v.vec[40] = -1
Draw(v, deformation=True, euler_angles=[-60,0,0],
     settings={"Misc": {"subdivision": 15, "line_thickness": 5, "fast_draw": True},
               "Colormap":{"autoscale": True, "ncolors": 16}});

The cubic DG space is made next and one of its shape functions is shown:

W3 = ng.L2(mesh, order=3)
w = GridFunction(W3)
w.vec[:] = 0;
w.vec[W3.ndof-2] = 0.3
Draw(w,deformation=True, euler_angles=[-80,0,-30],
     settings={"Misc": {"subdivision": 15, "line_thickness": 5, "fast_draw": True},
               "Colormap":{"autoscale": True, "ncolors": 16}});

In a previous notebook, we saw how to distinguish between coupling and local degrees of freedom. As you can imagine, the DG space, not having any continuity requirements, has no coupling degrees of freedom.

all( [W3.CouplingType(i) == ng.COUPLING_TYPE.LOCAL_DOF  # true if this holds for all dofs
      for i in range(W3.ndof)] )

True

In contrast, the CR space has no local degrees of freedom.

any( [U.CouplingType(i) == ng.COUPLING_TYPE.LOCAL_DOF 
      for i in range(U.ndof)] )

False

Questions for discussion:

  • For what degrees does the Lagrange space have local degrees of freedom in two space dimensions?

  • What would be your answer in higher dimensions?

  • How do we make finite element spaces on three-dimensional domains in ngsolve?

2.3. Best approximation#

There are multiple ways to approximate a given function f by finite element spaces.

  • One is by the Oswald approximation, implemented in ngsolve by the GridFunction.Set(f) method, as we have already seen.

  • Another is the canonical finite element interpolant, obtainable in ngsolve using the same Set method, but with an additional keyword argument: GridFunction.Set(f, dual=True).

  • While the above two are local constructs, there is a global approximation method which gives the best possible approximation to f from the finite element space, which we now discuss.

Given a function \(\newcommand{\om}{\varOmega}\) \(f \in L_2(\om)\), the best approximation to \(f\) from a finite element subspace \(S_h\) of \(L_2(\om)\) is obtained as the unique minimizer \(f_h \in S_h\) satisfying

\[ \| f - f_h \|^2_{L_2(\om)} = \min_{v \in S_h } \| f -v \|^2_{L_2(\om)}. \]

To find this minimizer, we can use ngsolve’s automatic Newton-based minimization facilities as shown in the next bit of code: behind the scenes, it automatically differentiates the quadratic functional \(Q(v) = \| f -v \|^2_{L_2(\om)}\) and uses the Newton iteration on the linearized differential to solve for the minimum. Note how the square of the \(L_2\) norm is expressed as an integral over the domain in the next code cell.

def MinDistL2(fespace, function):
    """
    Return the best approximation of "function" from the "fespace"
    where optimality is measured in the L2 norm.
    """

    # Set up the distance minimization 
    approx = fespace.TrialFunction()
    a = ng.BilinearForm(fespace)
    a += ng.Variation((function - approx)**2 * dx)

    # Solve for the minimizer
    best_approx = GridFunction(fespace)
    ng.solvers.Newton(a=a, u=best_approx, printing=False)
    
    return best_approx

Here approx is an ngsolve TrialFunction that takes the role of \(v\) in the minimization above, function is an expression representing the given function \(f\), dx tells ngsolve to integrate over all elements of the mesh, Variation asks ngsolve to compute the linearization automatically, and BilinearForm allows for storage and representation of the linearization. We will become much more familiar with some of these objects as the course proceeds. If you think that just calling Newton is too opaque, rest assured that we will see other more transparent ways to compute the minimizer in a later activity.

Let us now apply the above function to compute the best approximation to a given function from three finite element spaces.

This is also a good time to hint at the ease of specifying complex domains through the tools at your disposal now. The renowned OCCT kernel for geometric modeling from the Open Cascade company is made available in ngsolve via the netgen.occ module. Here is an example domain that goes a bit beyond the unit square.

from netgen.occ import WorkPlane, Glue, OCCGeometry

n = 8
wp = WorkPlane()

# make a polygon of n sides 
for i in range(n): 
    wp.Line(1).Rotate(360/n)
face = wp.Face()

# make three holes to cut out later
wp.MoveTo(0, 1.5)
eye1 = wp.Rectangle(l=0.3, w=0.2).Face()
wp.MoveTo(0.8, 1.5)
eye2 = wp.Rectangle(l=0.3, w=0.2).Face()
wp.MoveTo(0, 0.7)
smile = wp.Rectangle(l=1, w=0.15).Face()
wp.MoveTo(0.45, 1)

# make a subdomain (not to be cut out)
nose = wp.Rectangle(l=0.2, w=0.5).Face()
nose.name = 'nose'
nose.maxh = 0.05

# final structure, using union and set subtraction
domain = Glue([nose, face - eye1 - eye2 - smile])
mesh = Mesh(OCCGeometry(domain, dim=2).GenerateMesh(maxh=0.1))
# visualize: you can either Draw(mesh) or Draw(domain)

Using the names assigned while creating the domain, we can create a piecewise function which takes values of distinct expressions on distinct subdomains. Here is a simple example for this domain, which defines the function \(f\) we want to approximate using finite elements.

mesh.GetMaterials()  # recall names assigned during domain creation

('nose', 'default')

f = mesh.MaterialCF({'nose': exp(-100*((x-0.55)**2 + (y-1)**2/10))}, default=0)
Draw(f, mesh, deformation=True);

Here is the best approximation to this \(f\) from the DG space of degree one:

wh = MinDistL2(ng.L2(mesh, order=1),f)
print('Error in best approximation by DG = ', ng.Integrate((wh - f)**2, mesh))
Draw(wh, deformation=True,  euler_angles=[-70, 0, -70]);

Error in best approximation by DG =  6.717392133705651e-06

What if we use the Lagrange space of the same degree?

vh = MinDistL2(ng.H1(mesh, order=1), f)
print('Error in best approximation by H1 = ', ng.Integrate((vh - f)**2, mesh))
# Draw(vh, deformation=True) 

Error in best approximation by H1 =  0.0014162711410865595

And how does the CR space of the same degree perform?

uh = MinDistL2(ng.FESpace('nonconforming', mesh), f)
print('Error in best approximation by CR = ', ng.Integrate((uh - f)**2, mesh))
# Draw(uh, deformation=True)

Error in best approximation by CR =  0.0007872914948498292

Clearly, minimization over smaller subspaces result in larger errors. For this particular \(f\), the DG space is very well suited for approximation because \(f\) is discontinuous along a curve that is a union of some element interfaces (where DG functions are also allowed to be discontinuous).

Questions for discussion:

  • Is the best approximation mapping \(f \mapsto f_h\) an orthogonal projection?

  • Could we use a linear solver rather than a nonlinear solver to compute this map?

2.4. Constrained minimization#

Recall how local degrees of freedom were fused together to form the global CR space and that note that edge midpoint continuity is equivalent to the zero mean of the jump of the function across the edge. We thus recognize that the CR space can also be written as the following subspace the linear DG space:

\[ U_h = \{ w \in W_{h1}: \text{ jump of $w$ across every interior mesh edge has zero mean} \}. \]

The jump of \(w\) is a function on element interfaces. On an interface \(F\) shared by two elements \(K_\pm\), with an orienting normal \(n_F\), the jump ⟦ \( w \) ⟧ at a point \(x\) in \(F\) equals the difference of limits of \(w(x \pm \epsilon n_F)\) as \(\epsilon \to 0\). (Later we give an orientation-independent definition.) The jump of \(w\) is typically extended even to facets on the boundary by considering \(w\) as extended by zero outside the domain. Although the jump is typically denoted by ⟦ \( w \) ⟧, to ease notation, here we will simply refer to it as \(\mathrm{jump}(w)\). With this notation, the best approximation \(f_h\) to \(f\) from \(U_h\) is also the solution of a constrained optimization problem on \(W_{h1}\):

\[\begin{split} \begin{aligned} \| f - f_h\|_{L_2(\om)} & = \min_{u_h \in U_h} \| f - u_h\|_{L_2(\om)} \\ & = \min_{\overset{\scriptstyle{w_h \in W_{h1}}}{\int_F \mathrm{jump}(w_h)=0 \text{ on facets } F}} \| f - w_h\|_{L_2(\om)} \end{aligned} \end{split}\]

Thus, the unconstrained minimization over \(U_h\) is equivalent to a constrained minimization over \(W_{h1}\) with one constraint per edge \(F\).

From this perspective, it is natural to ask what happens if we constrain \(\mathrm{jump}(w_h)\) further, e.g., instead of just the zero-order moment (or the mean value), one could require further higher-order moments of the jump to vanish, i.e., what if we generalize \(U_h\) to

\[ U_{hpk} = \Big\{ w \in W_{hp}: \int_F \mathrm{jump}(w) q =0 \text{ for all } q \in P_k(F) \text{ on all interior edges } F\Big\}. \]

The \(p=1, k=0\) case gives \(U_h\).

Is it possible to compute the \(L_2\) best approximation of \(f\) from a generalized CR-like space \(U_{hpk}\)?

In the next code cell you will see a code that does the distance minimization on \(U_{hpk}\) via constrained optimization. The idea behind it is the same as what you have seen in calculus: to minimize a function \(F(x)\) over \(x \in \mathbb R^d\) subject to a constraint \(G(x)=0\), one introduces a Lagrange multiplier \(\lambda\) and find the stationary points of the augmented functional \(L(x, \lambda) = F(x) + \lambda G(x)\), often called the Lagrangian. We apply this idea with \(F\) equal to the distance functional and \(G\) set to the above jump constraint; the only difference here is that both the variable being optimized over (\(w\)) and the Lagrange multiplier (\(\lambda\)) are now functions in finite element spaces:

\[ L(w, \lambda) = \| f - w\|^2_{L_2(\om)} + \sum_F \int_F \lambda \; \mathrm{jump}(w)\, ds, \]

where \(w \in W_{hp}\) and \(\lambda|_F \in P_k(F)\) on each facet \(F\).

We see that the Lagrange multiplier \(\lambda\) should be set in a finite element space on facets. To implement the last term with \(\mathrm{jump}(w)\) on the facets, we use the NormalFacetFESpace provided by ngsolve. This space consists of vector functions \(\varLambda\) such that \(\varLambda|_F = \lambda|_F n_F\) on each facet \(F\), where \(n_F\) is a unit normal vector on \(F\) (of arbitrarily fixed orientation) and \(\lambda|_F\) is a scalar function in \(P_k(F)\). Using \(n\) (without subscript) to denote the outward unit normal on element boundaries, we can rewrite the last term above as

\[ \sum_F \int_F \lambda \; \mathrm{jump}(w)\, ds = \sum_K \int_{\partial K} \Lambda \cdot w n \, ds. \]

This is just a rearrangement of the sum over facets into a sum over element boundaries \(\partial K\), and can be immediately understood from the following illustration of two triangles \(K_\pm\) sharing a facet \(F\). Indeed, for any \(\lambda\) supported only on \(F\), the sum over integrals over \(\partial K_\pm\) of \(\Lambda \cdot w n\) gives the integral of \(\lambda \,\mathrm{jump}(w)\) on \(F\).

Facet shared by 2 triangles

You will see this rearranged form of the constraint implemented in the next code cell, which computes the best approximation from \(U_{hpk}\).

# Given p, k, compute best approximation from U_hpk
p = 3    
k = 1

W = ng.L2(mesh, order=p)                  # Approximation lives here 
L = ng.NormalFacetFESpace(mesh, order=k)  # Lagrange multiplier lives here
n = ng.specialcf.normal(mesh.dim)
WL = W * L   # Cartesian product of two finite element spaces

approx, multiplier = WL.TrialFunction()
distance2 = (f - approx)**2 
constraint = ng.InnerProduct(multiplier, approx*n)

a = ng.BilinearForm(WL)
a += ng.Variation(distance2 * dx + constraint * dx(element_boundary=True))
                  
best_approx = ng.GridFunction(WL) # beware of nonconvergence/singularity
ng.solvers.Newton(a=a, u=best_approx, printing=True)

uh2 = best_approx.components[0]
ng.Integrate((uh2 - f)**2, mesh)

Newton iteration  0
err =  0.21643588155795607
Newton iteration  1
err =  3.925567747942326e-17

0.00024705047161825136

Draw(uh2, mesh, deformation=True, euler_angles=[-70, 0, -70], 
     settings={"Misc": {"subdivision": 15, "line_thickness": 5, "fast_draw": False}});

Note that the constrained optimzation code need not succeed for all combinations of \(p\) and \(k\). I leave you to investigate this on your own.

2.5. Summary#

We have seen

  • the DG, Lagrange, and Crouzeix-Raviart finite element spaces on 2D domains,

  • how to optimize a nonlinear functional using a built-in Newton solver in ngsolve,

  • how to compute best approximations from spaces by minimizing distance,

  • some OpenCascade geometry modeling constructs,

  • unit outward normal vector on element boundaries,

  • use of a facet finite element space,

  • jump functions through element boundary sums,

  • how to integrate inside an element,

  • how to integrate over element boundaries,

  • product of finite element spaces,

  • constrained optimization by Lagrange multiplier finite elements.

Contents