Extended documentation for InteriorFacetBasis

[gatling-nrl]

import gmsh
import matplotlib.pyplot as plt
import meshio
import numpy as np
import skfem
from skfem.helpers import dot, grad
import skfem.io
import skfem.visuals.matplotlib

Let's project a goofey piecewise constant function into P0 on our familiar toy mesh. This function is contrived to index the elements for us, so we can see what happens when we build an interior facet basis.

mesh = skfem.MeshTri().refined(1)
cell_basis_p0 = skfem.Basis(mesh, skfem.ElementTriP0())
f = np.arange(cell_basis_p0.zeros().shape[0])

fig, ax = plt.subplots(figsize=(9,8))
skfem.visuals.matplotlib.plot(cell_basis_p0, f, ax=ax, colorbar=True, cmap='Set2')
skfem.visuals.matplotlib.draw(mesh, ax=ax)

skfem can construct a special kind of basis from the interior facets of this mesh. "Interior" facets are those that are an edge for two triangles. This can be contrasted with boundary facets which are an edge for only one triangle.

interior_facet_basis_p0 = skfem.InteriorFacetBasis(mesh, cell_basis_p0.elem)

This basis spans the same space as the CellBasis we constructed from the same mesh (since we used the same element), but has a different set of quadrature points.

fig, ax = plt.subplots(1, 2, figsize=(12,6))
ax[0].set_title('CellBasis quadrature points')
basis_p0_pts = cell_basis_p0.global_coordinates().value
skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
ax[0].plot(basis_p0_pts[0], basis_p0_pts[1], 'or')
ax[1].set_title('InteriorFacetBasis quadrature points')
interior_facet_basis_p0_pts = interior_facet_basis_p0.global_coordinates().value
skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
ax[1].plot(interior_facet_basis_p0_pts[0], interior_facet_basis_p0_pts[1], 'or')

Computing the integral of our function on the CellBasis is straight forward. The function has a unique value at every quadrature point. We can use the Functional decorator and assembly to compute this integral.

@skfem.Functional
def test1(w):
    return w['any_discrete_field']  # return the thing we want the integral of

skfem.asm(test1, cell_basis_p0, any_discrete_field=f)

3.5

This should be the area of each triangle times the value in that triangle, all summed up, since this is just a piecewise constant projection.

area_of_each_cell = 0.5*0.5*0.5 # 1/2 * base * height
(f * area_of_each_cell).sum()

3.5

Now suppose we want the integral along the interior edges only. In the present consideration, that's a single integral over length (as opposed to the double integral over area we just calculated for the CellBasis). This is were the InteriorFacetBasis comes in handy. We should be able to do the exact same assembly... except which value should be used for each edge? Comparing the annotations in the two plots below reveals the problem.

fbasis_s0 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=0)
fbasis_s1 = skfem.InteriorFacetBasis(mesh, skfem.ElementTriP0(), side=1)


def plot_annotations(ax):
    def _functional(w):
        pts = w.x.reshape(2,-1).T
        vals = w['any_discrete_field'].value.reshape(-1)
        for p, v in zip(pts, vals):
            ax.plot(p[0], p[1], 'or')
            ax.annotate(f'{v:.1f}', p)
        return w['any_discrete_field']
    return skfem.Functional(_functional)

fig, ax = plt.subplots(1, 2, figsize=(12,6))
skfem.visuals.matplotlib.plot(cell_basis_p0, f, ax=ax[0], cmap='Set2')
skfem.asm(plot_annotations(ax[0]), fbasis_s0, any_discrete_field=f)
skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
skfem.visuals.matplotlib.plot(cell_basis_p0, f, ax=ax[1], cmap='Set2')
skfem.asm(plot_annotations(ax[1]), fbasis_s1, any_discrete_field=f)
skfem.visuals.matplotlib.draw(mesh, ax=ax[1])

This snippet also shows how the side argument is used when constructing an InteriorFacetBasis. It is important to note that this class doesn't directly allow specifying which side to use in a facet-by-facet manner. Nevertheless, changing side from 0 to 1 will change every the element association for every facet in the set. So, while in theory these integrals over interior facets form a combintorial problem, skfem makes two (arbitrary) ones readily available:

print('using side==0:', skfem.asm(test1, fbasis_s0, any_discrete_field=f))
print('using side==1:', skfem.asm(test1, fbasis_s1, any_discrete_field=f))

using side==0: 18.606601717798213
using side==1: 19.19238815542512

We can easily check these values using some visual inspection of the above plots and integrating over the line segments, keeping in mind for this toy problem everything is piecewise constant.

print('using side==0:', 0.5*np.sqrt(2)*np.array([6,3,5,1]).sum() + 0.5*np.array([4,3,7,2]).sum())
print('using side==1:', 0.5*np.sqrt(2)*np.array([0,4,7,2]).sum() + 0.5*np.array([6,7,1,6]).sum())

using side==0: 18.606601717798213
using side==1: 19.19238815542512

[gatling-nrl]

Another very important feature of the InteriorFacetBasis class is its concept of a normal vector. This is available when defining forms and functionals as w['n'] and w.n. Note, however, that interior facets do not have a consistent way to define the normal. There is no notion of "outward", so in that sense the normal vectors are arbitrarily chosen from the two directions associated with each facet. We can plot the ones associated with our example above.

def plot_normal_arrows(ax):
    arrow_style = dict(width=.0, length_includes_head=True, head_width=.025, head_length=.025, facecolor='k', fill=True)
    def _form(w):
        print('w.x.shape, w.n.shape =', w.x.shape, w.n.shape)
        pp = w.x.reshape(2,-1).T  # array is shaped to group by facet, but we just want all the points by x,y
        nn = w.n.reshape(2,-1).T  # array is shaped to group by facet, but we just want all the normals by dx,dy
        for p, n in zip(pp, nn):
            # when calling functions, (...,*arg,...) will expand the list `arg` as if we had used `...,arg[0],arg[1],etc,...`
            ax.arrow(*p, *(0.075*n), **arrow_style)  # 0.075 so arrows look nice on this plot
            ax.plot(p[0], p[1], 'or')
        return w.x[0]
    return skfem.Functional(_form)

fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].set_title('normal vectors with side==0')
skfem.visuals.matplotlib.plot(cell_basis_p0, f, ax=ax[0], cmap='Set2')
skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
skfem.asm(plot_normal_arrows(ax[0]), fbasis_s0)
ax[1].set_title('normal vectors with side==1')
skfem.visuals.matplotlib.plot(cell_basis_p0, f, ax=ax[1], cmap='Set2')
skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
skfem.asm(plot_normal_arrows(ax[1]), fbasis_s1)

w.x.shape, w.n.shape = (2, 8, 2) (2, 8, 2)
w.x.shape, w.n.shape = (2, 8, 2) (2, 8, 2)

This snippet also shows how w.x and w.n are related. They will always be the same shape, with the first dimension of w.x corresponding to x and y and the first dimension of w.n corresponding to dx and dy. (Or likewise in higher or lower dimensional mesh spaces.)

[gdmcbain]

I like these pictures; might be good for rendering ElementTriRT0 #824. But why do the two look the same? I thought the arrows would be flipped.

[kinnala]

Normal vector is always outward from side=0.