Best approximations in 2D

2. Best approximations in 2D#

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 best approximations from these spaces. Best approximations tell us how good a finite element space approximates functions outside it.

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

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 later.

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

BaseWebGuiScene

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, typically denoted by ⟦ $ w $ ⟧, but here we will denote it simply by $\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, 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$.

Questions for discussion:

How to transform the condition on $\mathrm{jump}(w)$ to a sum over element boundaries?
Is it possible to compute the $L_2$ best approximation of $f$ from a generalized CR-like space $U_{hpk}$?

Here is a code that does the distance minimization on $U_{hpk}$ via constrained optimization. We use a Lagrange multiplier in a finite element space on the facets. To implement $\mathrm{jump}(w)$ on the facets, we use an indirect method that produces the jump while summing over elements and integrating over boundaries of each element. Note that the optimization need not succeed for all combinations of $p$ and $k$. I leave you to investigate this on your own.

# Compute best approximation from U_hpk
p = 3    
k = 1

W = ng.L2(mesh, order=p)                  # Approximation lives here 
L = ng.NormalFacetFESpace(mesh, order=1)  # 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.939612696576264e-17

0.0002470504716182514

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

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.

Best approximations in 2D

Contents

2. Best approximations in 2D#

2.1. Shape functions in the lowest order case#

2.2. Higher order cases and coupling dofs#

2.3. Best approximation#

2.4. Constrained minimization#

2.5. Summary#