This page was generated from unit-2.13-interfaces/interfaceresistivity.ipynb.

2.13 Interface resisitivity

We model a problem with resistivity interface conditions \begin{eqnarray*} \lambda_l \frac{\partial u_l}{\partial n_l} + \lambda_r \frac{\partial u_r}{\partial n_r} & = & 0 \\ \alpha (u_l - u_r) & = & \lambda_l \frac{\partial u_l}{\partial n_l} \end{eqnarray*}

This is easily modeled by an additional Robin-like term on the interface:

\[\int_\Omega \lambda \nabla u \nabla v \, dx + \int_\gamma \alpha [u][v] \, ds,\]

where \([u] = u_l - u_r\) denotes the jump across the interface.

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

inner = Rectangle(1,1).Face()
inner.edges.name="interface"
inner.edges.Max(X).name="interface_right"

outer = MoveTo(-1,-1).Rectangle(3,3).Face()
outer.edges.name="dir"
outer = outer-inner

inner.faces.name="inner"
outer.faces.name="outer"
geo = Glue ([inner,outer])
mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.2))

Draw (mesh);

print (mesh.GetMaterials(), mesh.GetBoundaries())

('inner', 'outer') ('interface', 'interface_right', 'interface', 'interface', 'dir', 'dir', 'dir', 'dir')

Interfaces without boundaries

In the first case, we assume the interface is the boundary of the inner domain. For this case, we define spaces on both sides separately, and link them via the Robin - term.

In our problem we choose zero conductivity across the right edge of the interface.

fesi = H1(mesh, order=3, definedon=mesh.Materials("inner"))
feso = H1(mesh, order=3, definedon=mesh.Materials("outer"), dirichlet="dir")
fes = fesi*feso

(ui,uo), (vi,vo) = fes.TnT()

alpha = 10

a = BilinearForm(fes)
a += grad(ui)*grad(vi)*dx("inner")
a += grad(uo)*grad(vo)*dx("outer")
a += alpha * (ui-uo) * (vi-vo) * ds("interface")
a.Assemble()

f = LinearForm(1*vi*dx("inner")).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec

sol = mesh.MaterialCF({"inner":gfu.components[0], "outer":gfu.components[1]})
Draw (sol, mesh);

Interfaces with boundaries

Now we model an interface resistivity only at the right edge. For this we need a space discontinuous only at that edge. We obtain that by enriching an H1-space with functions defined by values on the interface, and supported only on one side. The jump is exactly the enrichment function.

We define the enrichment space first on one domain, and then restrict degrees of freedom to the interior of the right interface edge.

fes1 = H1(mesh, order=3, dirichlet="dir")

fesext = H1(mesh, order=3, definedon="outer")
actdofs = fesext.GetDofs(mesh.Boundaries("interface_right")) & \
    ~fesext.GetDofs(mesh.Boundaries("interface"))
fesext = Compress (fesext, active_dofs=actdofs)
print (actdofs)

0: 00000000000011110000000000000000000000000000000000
50: 00000000000000000000000000000000000000000000000000
100: 00000000000000000000000000000000000000000000000000
150: 00000000000000000000000000000000000000000000000000
200: 00000000000000000000000000000000000000000000000000
250: 00000000000000000000000000000000000000000000000000
300: 00110000001100000000000000000000000000000000000000
350: 00000000000000000000110000110000110000000000000000
400: 00000000000000000000000000000000000000000000000000
450: 00000000000000000000000000000000000000000000000000
500: 00000000000000000000000000000000000000000000000000
550: 00000000000000000000000000000000000000000000000000
600: 00000000000000000000000000000000000000000000000000
650: 00000000000000000000000000000000000000000000000000
700: 00000000000000000000000000000000000000000000000000
750: 00000000000000000000000000000000000000000000000000
800: 00000000000000000000000000000000000000000000000000
850: 00000000000000000000000000000000000000000000000000
900: 00000000000000000000000000000000000000000000000000
950: 00000000000000000000000000000000000000000000000000
1000: 00000000000000000000000000000000000000000000000000
1050: 00000000000000000000000000000000000000000000000000
1100: 00000000000000000000000000000000000000000000000000
1150: 00000000000000000000000000000000000000000000000000
1200: 00000000000000000000000000000000000000000000000000
1250: 00000000000000000000000000000000000000000000000000
1300: 00000000000000000000000000000000000000000000000000
1350: 00000000000000000000000000000000000000000000000000
1400: 00000000000000000000000000000000000000000000000000
1450: 00000000000000000000000000000000000000000000000000
1500: 00000000000000000000000000000000000000000000000000
1550: 00000000000000000000000000000000000000000000000000
1600: 00000000000000000000000000000000000000000000000000
1650: 00000000000000000000000000000000000000000000000000
1700: 00000000000000000000000000000000000000000000000000
1750: 00000000000000000000000000000000000000000000000000
1800: 00000000000000000000000000000000000000000000000000
1850: 00000000000000000000000000000000000000000000000000
1900: 00000000000000000000000000000000000000000000000000
1950: 00000000000000000000000000000000000000000000000000
2000: 00000000000000000000000000000000000000000000000000
2050: 00000000000000000000000000000000000000000000000000
2100: 00000000000000000000000000000000000000000000000000
2150: 00000000000000000000000000000000000000000000000000
2200: 0000000000000000000000000000000000

A typical enrichment function looks like

gfext = GridFunction(fesext)
gfext.Set(sin(10*y))
Draw (gfext);

fes = fes1 * fesext
(u,uext), (v,vext) = fes.TnT()

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx("inner")
a += (grad(u)+grad(uext))* (grad(v)+grad(vext))*dx("outer")
a += 3 * uext*vext * ds("interface_right")
a.Assemble()

f = LinearForm(1*v*dx("inner")).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec

sol = mesh.MaterialCF({"inner":gfu.components[0], "outer":gfu.components[0]+gfu.components[1]})
Draw (sol, mesh);