SubdomainFacetBasis in action

[gatling-nrl]

Some snippets and plots to demonstrate SubdomainFacetBasis.
Edit: This class is in #851

The toy problem here of 2D concentric circles. The inner one is at 1 V and the outer one at 0 V. The conductivity (sigma) is specified between them, 1e6 for the conductor and 1 between the conductor and the outside wall. I want to know the current (in A/m since this is 2D) between the inner and outer ring. Note that as long as the inner boundary is enclosed, it shouldn't matter to the integral what path encloses it. For the stated values, the analytic answer is 21.84072595706403 and COMSOL calculates 21.840744824811342. COMSOL gets this kind of accuracy even with a very coarse mesh, but I don't know the element its using (or any other trickery its up to).

The mesh

The mesh is rather coarse, and I'm only going to use quadratic elements here, giving up on accuracy but making it easier to draw some of the plots for this demonstration.

gmsh.initialize()
try:
    gmsh.option.setNumber("Mesh.MeshSizeFactor", 1)
    r1_tag = gmsh.model.occ.add_disk(xc=0, yc=0, zc=0, rx=r1, ry=r1)
    r2_tag = gmsh.model.occ.add_disk(xc=0, yc=0, zc=0, rx=r2, ry=r2)
    r3_tag = gmsh.model.occ.add_disk(xc=0, yc=0, zc=0, rx=r3, ry=r3)
    r4_tag = gmsh.model.occ.add_disk(xc=0, yc=0, zc=0, rx=r4, ry=r4)
    r5_tag = gmsh.model.occ.add_disk(xc=0, yc=0, zc=0, rx=r5, ry=r5)

    gmsh.model.occ.cut([(2,r5_tag)], [(2,r4_tag)], removeTool=False)
    gmsh.model.occ.cut([(2,r4_tag)], [(2,r3_tag)], removeTool=False)
    gmsh.model.occ.cut([(2,r3_tag)], [(2,r2_tag)], removeTool=False)
    gmsh.model.occ.cut([(2,r2_tag)], [(2,r1_tag)], removeTool=False)

    gmsh.model.occ.synchronize()
    # Label the regions
    gmsh.model.set_physical_name(
        TWO_DIMENSIONAL,
        gmsh.model.add_physical_group(TWO_DIMENSIONAL, tags=[r5_tag]),
        'r5',
    )
    gmsh.model.set_physical_name(
        TWO_DIMENSIONAL,
        gmsh.model.add_physical_group(TWO_DIMENSIONAL, tags=[r4_tag]),
        'r4',
    )
    gmsh.model.set_physical_name(
        TWO_DIMENSIONAL,
        gmsh.model.add_physical_group(TWO_DIMENSIONAL, tags=[r3_tag]),
        'r3',
    )
    gmsh.model.set_physical_name(
        TWO_DIMENSIONAL,
        gmsh.model.add_physical_group(TWO_DIMENSIONAL, tags=[r2_tag]),
        'r2',
    )
    gmsh.model.set_physical_name(
        TWO_DIMENSIONAL,
        gmsh.model.add_physical_group(TWO_DIMENSIONAL, tags=[r1_tag]),
        'r1',
    )
    # Generate and save
    gmsh.model.mesh.generate(2)
    gmsh.write('circle_domain.msh')
finally:
    gmsh.clear()
    gmsh.finalize()

circle_domain.zip

mesh = skfem.io.meshio.from_file('circle_domain.msh', None)
cbasis_p0 = skfem.CellBasis(mesh, skfem.ElementTriP0())
subdomains = cbasis_p0.zeros()
subdomains[cbasis_p0.get_dofs(elements='r1')] = 1
subdomains[cbasis_p0.get_dofs(elements='r2')] = 2
subdomains[cbasis_p0.get_dofs(elements='r3')] = 3
subdomains[cbasis_p0.get_dofs(elements='r4')] = 4
subdomains[cbasis_p0.get_dofs(elements='r5')] = 5

fig,ax = plt.subplots(1,1,figsize=(6,5))
skfem.visuals.matplotlib.plot(cbasis_p0, subdomains, ax=ax, cmap='Set1', colorbar=True, vmin=1, vmax=10)
skfem.visuals.matplotlib.draw(mesh, ax=ax)

The solution to Laplace's equation.

electrode = ['r1', 'r2', 'r3']
medium = ['r4', 'r5']
cbasis = skfem.CellBasis(mesh, skfem.ElementTriP2())
cbasis_p0 = cbasis.with_element(skfem.ElementTriP0())

sigma = cbasis_p0.zeros()
sigma[cbasis_p0.get_dofs(elements=electrode)] = 1000
sigma[cbasis_p0.get_dofs(elements=medium)] = 1

@skfem.BilinearForm
def laplace(u, v, w):
    return dot(w.sigma*grad(u), grad(v))

A = laplace.assemble(cbasis, sigma=cbasis_p0.interpolate(sigma))
u0 = cbasis.zeros()
u0[cbasis.get_dofs(elements=electrode)] = 1
u0_dofs = cbasis.get_dofs() + cbasis.get_dofs(elements='r1')
A,b = skfem.enforce(A, D=u0_dofs, x=u0)
u = skfem.solve(A, b)

fig,ax = plt.subplots(1,2,figsize=(10,5))
ax[0].set_title('sigma')
skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
skfem.visuals.matplotlib.plot(cbasis_p0, sigma, ax=ax[0], cmap='Set3')
ax[1].set_title('solution')
skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
skfem.visuals.matplotlib.plot(cbasis, u, ax=ax[1], cmap='viridis', shading='gouraud')

Calculate Current

The now-familiar asm side effects plotter.

def calculate_current_and_plot_normals(ax):
    def _form(w):
        arrow_style = dict(width=.0, length_includes_head=True, head_width=.025, head_length=.025, facecolor='k', fill=True)
        pp = w.x.reshape(2,-1).T
        nn = w.n.reshape(2,-1).T
        for p, n in zip(pp, nn):
            ax.arrow(*p, *(0.05*n), **arrow_style)  # 0.075 so arrows look nice on this plot
        # calculate current
        return dot(w.n, -w.sigma*grad(w.u))
    return skfem.Functional(_form)

The plotting code

def plot_subdomain_facet_basis(subdomain):
    fig,ax = plt.subplots(1,2,figsize=(12,6))
    skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
    skfem.visuals.matplotlib.plot(cbasis, u, ax=ax[0], cmap='viridis', shading='gouraud')
    skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
    skfem.visuals.matplotlib.plot(cbasis, u, ax=ax[1], cmap='viridis', shading='gouraud')

    sfbasis = skfem.SubdomainFacetBasis(mesh, skfem.ElementTriP2(), elements=subdomain, side=0)
    sfbasis_p0 = sfbasis.with_element(skfem.ElementTriP0())
    calc1 = calculate_current_and_plot_normals(ax[0]).assemble(
        sfbasis, sigma=sfbasis_p0.interpolate(sigma), u=sfbasis.interpolate(u))
    ax[0].set_title(f"side==0\nintegral={calc1}")

    sfbasis = skfem.SubdomainFacetBasis(mesh, skfem.ElementTriP2(), elements=subdomain, side=1)
    sfbasis_p0 = sfbasis.with_element(skfem.ElementTriP0())
    calc2 = calculate_current_and_plot_normals(ax[1]).assemble(
        sfbasis, sigma=sfbasis_p0.interpolate(sigma), u=sfbasis.interpolate(u))
    ax[1].set_title(f"side==1\nintegral={calc2}")

And the results.

plot_subdomain_facet_basis('r1')

Note how the innermost subdomain is used as a Dirichlet boundary condition and should have zero gradient (and thus that integral should be zero, which it is, yay!)

plot_subdomain_facet_basis(['r1', 'r2'])

Notwithstanding the loss in accuracy from the coarse mesh and quadratic element, the current is equal and opposite. It will be this way for any closed contour around the source.

For example, I can also integrate around a really wonky boundary and get the same result.

plot_subdomain_facet_basis(np.arange(86))

Finally, integrating around a contour that doesn't enclose a source gives zero. Hurray Gauss's Law.

plot_subdomain_facet_basis('r3')