Tutorial 4.2: Heavy Sampling of an Intricate Edge¶

We may sometimes have an edge that has fine details that need to be resolved by increasing the sampling parameter $n$, with the edge being sampled at $2n+1$ points, including the end points.

When Things Go Right¶

For example, consider a unit square with one of the edges being sinusoidal.

In [1]:

import puncturedfem as pf

# define vertices
verts: list[pf.Vert] = []
verts.append(pf.Vert(x=0.0, y=0.0))
verts.append(pf.Vert(x=1.0, y=0.0))
verts.append(pf.Vert(x=1.0, y=1.0))
verts.append(pf.Vert(x=0.0, y=1.0))

# define edges
edges: list[pf.Edge] = []
edges.append(
    pf.Edge(
        verts[0],
        verts[1],
        pos_cell_idx=0,
        curve_type="sine_wave",
        amp=0.1,
        freq=4,
    )
)
edges.append(pf.Edge(verts[1], verts[2], pos_cell_idx=0))
edges.append(pf.Edge(verts[2], verts[3], pos_cell_idx=0))
edges.append(pf.Edge(verts[3], verts[0], pos_cell_idx=0))

# define mesh cell
K_simple = pf.MeshCell(idx=0, edges=edges)

# parameterize edges
K_simple.parameterize(quad_dict=pf.get_quad_dict(n=64))

# set up Nystrom solver
nyst = pf.NystromSolver(K_simple, debug=True)

# plot boundary
pf.plot.MeshPlot(K_simple.get_edges()).draw()

import puncturedfem as pf # define vertices verts: list[pf.Vert] = [] verts.append(pf.Vert(x=0.0, y=0.0)) verts.append(pf.Vert(x=1.0, y=0.0)) verts.append(pf.Vert(x=1.0, y=1.0)) verts.append(pf.Vert(x=0.0, y=1.0)) # define edges edges: list[pf.Edge] = [] edges.append( pf.Edge( verts[0], verts[1], pos_cell_idx=0, curve_type="sine_wave", amp=0.1, freq=4, ) ) edges.append(pf.Edge(verts[1], verts[2], pos_cell_idx=0)) edges.append(pf.Edge(verts[2], verts[3], pos_cell_idx=0)) edges.append(pf.Edge(verts[3], verts[0], pos_cell_idx=0)) # define mesh cell K_simple = pf.MeshCell(idx=0, edges=edges) # parameterize edges K_simple.parameterize(quad_dict=pf.get_quad_dict(n=64)) # set up Nystrom solver nyst = pf.NystromSolver(K_simple, debug=True) # plot boundary pf.plot.MeshPlot(K_simple.get_edges()).draw()

Setting up Nyström solver... 512 sampled points on 4 edges
	debug-NystromSolver: Computing condition number...
	debug-NystromSolver: Condition number = 4.44e+02

It is simple to verify that $v\in V_1(K)$ given by $v(x_1,x_2) = x_2$ has a square $L^2$ norm of \begin{align*} \int_K v^2 ~dx = \frac13~. \end{align*} Let's verify this:

In [2]:

# define v to have a Dirichlet trace of x_2 on each edge
x2 = pf.Polynomial([(1.0, 0, 1)])
v_trace = pf.DirichletTrace(edges=K_simple.get_edges(), funcs=x2)

# the local function v = x_2 is harmonic
v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False)

# compute area and error
L2_exact = 1 / 3
L2_computed = v.get_l2_inner_prod(v)
print(f"Error = {abs(L2_exact - L2_computed):.4e}")

# define v to have a Dirichlet trace of x_2 on each edge x2 = pf.Polynomial([(1.0, 0, 1)]) v_trace = pf.DirichletTrace(edges=K_simple.get_edges(), funcs=x2) # the local function v = x_2 is harmonic v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False) # compute area and error L2_exact = 1 / 3 L2_computed = v.get_l2_inner_prod(v) print(f"Error = {abs(L2_exact - L2_computed):.4e}")

Error = 9.1862e-12

When Things Go Wrong¶

Let's make this example more interesting by increasing the frequency of the sinusoid on the bottom of the square.

In [3]:

# crazy edge
edges[0] = pf.Edge(
    verts[0],
    verts[1],
    pos_cell_idx=0,
    curve_type="sine_wave",
    amp=0.1,
    freq=32,  # increase frequency
)

# define and parameterize a new mesh cell
K = pf.MeshCell(idx=0, edges=edges)
K.parameterize(quad_dict=pf.get_quad_dict(n=64))

# and look at it
pf.plot.MeshPlot(K.get_edges()).draw()

# crazy edge edges[0] = pf.Edge( verts[0], verts[1], pos_cell_idx=0, curve_type="sine_wave", amp=0.1, freq=32, # increase frequency ) # define and parameterize a new mesh cell K = pf.MeshCell(idx=0, edges=edges) K.parameterize(quad_dict=pf.get_quad_dict(n=64)) # and look at it pf.plot.MeshPlot(K.get_edges()).draw()

That doesn't look right... We can change the sampling parameter $n$ when initializing a MeshPlot instance to get more resolution. We also need to set the reparameterize flag to True.

In [4]:

pf.plot.MeshPlot(K.get_edges(), reparameterize=True, n=512).draw()

pf.plot.MeshPlot(K.get_edges(), reparameterize=True, n=512).draw()

/home/docs/checkouts/readthedocs.org/user_builds/puncturedfem/envs/stable/lib/python3.11/site-packages/puncturedfem/mesh/quad.py:135: UserWarning: Quad: n > 128 may cause numerical instability
  warn("Quad: n > 128 may cause numerical instability")

That looks pretty good, but note that MeshPlot didn't overwrite the sampled points we got above with n=64:

In [5]:

print(f"n = {K.num_pts // K.num_edges // 2}")

print(f"n = {K.num_pts // K.num_edges // 2}")

n = 64

Since this is not a high enough sampling rate to capture the high frequency of the bottom edge, we might expect our computation of the area to not be very accurate. Let's confirm this suspicion:

In [6]:

# set up Nystrom solver
nyst = pf.NystromSolver(K, debug=True)

# the harmonic function v = x_2
v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False)

# compute square L^2 norm and error
L2_exact = 1 / 3
L2_computed = v.get_l2_inner_prod(v)
print(f"Error = {abs(L2_exact - L2_computed):.4e}")

# set up Nystrom solver nyst = pf.NystromSolver(K, debug=True) # the harmonic function v = x_2 v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False) # compute square L^2 norm and error L2_exact = 1 / 3 L2_computed = v.get_l2_inner_prod(v) print(f"Error = {abs(L2_exact - L2_computed):.4e}")

Setting up Nyström solver... 512 sampled points on 4 edges
	debug-NystromSolver: Computing condition number...
	debug-NystromSolver: Condition number = 1.05e+05
warn-NystromSolver: GMRES failed to converge after 5120 iterations, residual norm = 1.50e-01
Error = 1.2192e-01

One might expect that if we increase the sampling parameter, this error will get smaller. However, we soon discover that this crashes the NystromSolver initialization.

In [7]:

# get 1024 sampled points on each edge
K.parameterize(quad_dict=pf.get_quad_dict(n=512))

try:
    # (WARNING!) this line will result in an exception being thrown
    nyst = pf.NystromSolver(K, debug=True)
except ZeroDivisionError as e:
    print("Indeed, an exception was thrown!\n", e)

# get 1024 sampled points on each edge K.parameterize(quad_dict=pf.get_quad_dict(n=512)) try: # (WARNING!) this line will result in an exception being thrown nyst = pf.NystromSolver(K, debug=True) except ZeroDivisionError as e: print("Indeed, an exception was thrown!\n", e)

Setting up Nyström solver... 4096 sampled points on 4 edges
Indeed, an exception was thrown!
 Nystrom system could not be constructed

Changing the Kress parameter (optional)¶

As we saw in Tutorial 1.1, we can change the Kress parameter $p$ to adjust how much the sampled points are "clustered" near the endpoints. The default value is $p=7$, but changing this to its lowest value $p=2$ results in sampled points that are more spread out, perhaps enough so that we can avoid division by machine zero.

NOTE: The condition number of the Nyström matrix is very high and GMRES will not converge quickly, if at all. Uncomment the following cell to see this.

In [8]:

# # get 1024 sampled points on each edge with lower Kress parameter
# K.parameterize(quad_dict=pf.get_quad_dict(n=512, p=2))
# nyst = pf.NystromSolver(K, debug=True)

# # the harmonic function v = x_2
# v = pf.LocalFunction(nyst=nyst, trace=v_trace)

# # (WARNING!) this line will take a long time to run
# v.compute_all()

# # compute square L^2 norm and error
# L2_exact = 1 / 3
# L2_computed = v.get_l2_inner_prod(v)
# print(f"Error = {abs(L2_exact - L2_computed):.4e}")

# # get 1024 sampled points on each edge with lower Kress parameter # K.parameterize(quad_dict=pf.get_quad_dict(n=512, p=2)) # nyst = pf.NystromSolver(K, debug=True) # # the harmonic function v = x_2 # v = pf.LocalFunction(nyst=nyst, trace=v_trace) # # (WARNING!) this line will take a long time to run # v.compute_all() # # compute square L^2 norm and error # L2_exact = 1 / 3 # L2_computed = v.get_l2_inner_prod(v) # print(f"Error = {abs(L2_exact - L2_computed):.4e}")

Splitting Edges¶

As we saw in Example 0.1, we can split edges in two using the split_edge() function. Let's try splitting the 'bad' edge into smaller edges.

In [9]:

# replace edge 0 with eight new edges
edges += pf.split_edge(edges[0], num_edges=8)
del edges[0]

# define mesh cell
K = pf.MeshCell(idx=0, edges=edges)

# replace edge 0 with eight new edges edges += pf.split_edge(edges[0], num_edges=8) del edges[0] # define mesh cell K = pf.MeshCell(idx=0, edges=edges)

In the previous section, we tried sampling each edge with $2n = 1024$ points. Notice, though, that only the bottom edge is problematic, and we might get away with sampling the straight edges at a lower rate. To keep the number of sampled points on the bottom edge the same, which has now been split into 8 edges, we need to set the sampling parameter to $n=64=512/8$.

In [10]:

# bottom edge sampled at 1024 points
K.parameterize(quad_dict=pf.get_quad_dict(n=64))

# set up Nystrom solver
nyst = pf.NystromSolver(K, debug=True)

# bottom edge sampled at 1024 points K.parameterize(quad_dict=pf.get_quad_dict(n=64)) # set up Nystrom solver nyst = pf.NystromSolver(K, debug=True)

Setting up Nyström solver... 1408 sampled points on 11 edges
	debug-NystromSolver: Computing condition number...
	debug-NystromSolver: Condition number = 1.69e+04

The NystromSolver didn't crash this time. Let's define the local function $v = x_2$ and take a peek at its trace:

In [11]:

# Dirichlet trace of the harmonic function v = x_2
x2 = pf.Polynomial([(1.0, 0, 1)])
v_trace = pf.DirichletTrace(K, funcs=x2)

# the harmonic function v = x_2
v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False)

# plot the trace of v
pf.plot.TracePlot(v_trace, K, quad_dict=pf.get_quad_dict(n=64)).draw()

# Dirichlet trace of the harmonic function v = x_2 x2 = pf.Polynomial([(1.0, 0, 1)]) v_trace = pf.DirichletTrace(K, funcs=x2) # the harmonic function v = x_2 v = pf.LocalPoissonFunction(nyst=nyst, trace=v_trace, evaluate_interior=False) # plot the trace of v pf.plot.TracePlot(v_trace, K, quad_dict=pf.get_quad_dict(n=64)).draw()

Finally, let's see if we can accurately compute our quantity of interest:

In [12]:

L2_exact = 1 / 3
L2_computed = v.get_l2_inner_prod(v)
print(f"Error = {abs(L2_exact - L2_computed):.4e}")

L2_exact = 1 / 3 L2_computed = v.get_l2_inner_prod(v) print(f"Error = {abs(L2_exact - L2_computed):.4e}")

Error = 4.5874e-06