33. Boundary Conditions

Usually, Dirichlet boundary conditions (essential bc) are built into the space: find \(u \in H^1\) such that \(u = u_D\) on \(\Gamma_D\). Now, we want to pose the Dirichlet bc by an extra equation. For this, we start from the strong form

\[ -\Delta u = f, \]

multiply the equation by a test-function in \(H^1\), integate by parts, and keep the boundary term:

\[ \int \nabla u \nabla v - \int_{\partial \Omega} \frac{\partial u}{\partial n} v = \int f v \]

We introduce a new variable \(\lambda\) for \(-\frac{\partial u}{\partial n}\) on the Dirichlet boundary. Natural bc are treated as usual, we assume we only have homogeneous Neumann bc for ease of notation:

\[ \int \nabla u \nabla v + \int_{\Gamma_D} \lambda v = \int fv \]

The Dirichlet bc \(u = u_D\) is now enforced also by a test function \(\mu\) living on the Dirichlet boundary. Thus, the whole equation is now:

Find \(u \in V := H^1(\Omega)\) and \(\lambda \in Q := H^{-1/2}(\Gamma_D)\) such that

\[\begin{split} \begin{array}{ccccll} \int_\Omega \nabla u \nabla v & + & \int_{\Gamma_D} \lambda v & = & \int_\Omega f v & \forall \, v \in V \\ \int_{\Gamma_D} \mu u & & & = & \int_{\Gamma_D} \mu u_D & \forall \, \mu \in Q \end{array} \end{split}\]

Why do we have the space \(H^{-1/2}(\Gamma_D)\) ? Functions from \(H^1(\Omega)\) have boundary values (the so called trace) exactly in the space \(H^{1/2}(\Gamma_D)\). We can pair these functions with elements from its dual space, called \(H^{-1/2}(\Gamma_D)\). To be precise, the integral is a convenient notation for the duality pairing:

\[ \int_{\Gamma_D} u \mu \quad \text{in sense of} \quad \left< u|_{\Gamma_D}, \mu \right>_{H^{1/2}(\Gamma_D) \times H^{-1/2}(\Gamma_D)} \]

The space \(H^{-1/2}\) is weaker (i.e. larger) than \(L_2\). Thus, we can use discontinuous \(L_2\) finite elements for its discretization.

from ngsolve import * 
from ngsolve.webgui import Draw
# mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesh = Mesh(unit_cube.GenerateMesh(maxh=0.2))

order=3
V = H1(mesh, order=order, orderface=order+1)
# Q = SurfaceL2(mesh, order=order-2)
Q = H1(mesh, order=order, definedon=mesh.Boundaries(".*"))
X = V*Q
print ("V.ndof =", V.ndof, "Q.ndof =", Q.ndof)

V.ndof = 6542 Q.ndof = 1652

WARNING: kwarg 'orderface' is an undocumented flags option for class <class 'ngsolve.comp.H1'>, maybe there is a typo?

u,lam = X.TrialFunction()
v,mu = X.TestFunction()
a = BilinearForm(X)
a += grad(u)*grad(v)*dx + (u*mu+v*lam)*ds

f = LinearForm(X)
f += 10*x*v*dx + y*mu*ds

a.Assemble()
f.Assemble()

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

sol_u, sol_lam = gfu.components
Draw (sol_u, mesh, clipping={ "function" : True })
Draw (sol_lam, mesh);

Choosing the test-function \(v = 1\) in the first equation $\( \int_\Omega \nabla u \nabla 1 + \int_{\Gamma_D} \lambda 1 = \int_\Omega f 1 \)$ we observe that the total flux is exactly in balance with the total source

\[ \int_{\Gamma_D} \lambda = \int_{\Omega} f \]

We compute the integral over the whole boundary as

Integrate(sol_lam, mesh, BND)

4.999999999999967

bndparts = Integrate(sol_lam, mesh, BND, region_wise=True)
print ("boundary parts:", bndparts)
print ("sum: ", sum(bndparts))

boundary parts:  0.31624
 1.34946
 1.81158
 -0.145481
  0.8347
 0.833505

sum:  4.999999999999967

Can we use continuous finite elements for the normal derivative \(\lambda\) ?

33.1. Interface conditions

Instead of enforcing Dirichlet boundary conditions, one can use the Lagrange multiplier also to enforce continuity across interfaces. Let the domain be \(\Omega_1 \cup \Omega_2\), and \(\gamma\) the common interface.

The interface conditions are continuity of the function, and continuity of the normal derivative:

\[\begin{align*} u_1 &= u_2 \\ \lambda \frac{\partial u_1}{\partial n_1} & = - \lambda_2 \frac{\partial u_2}{\partial n_2} \end{align*}\]

Define \(V = H^1 (\Omega_1 \cup \Omega_2)\), and \(Q = H^{-1/2}(\gamma)\), and search for \(u \in V\) and \(\lambda \in Q\) such that:

\[\begin{split} \begin{array}{ccccll} \int_{\Omega \cup \Omega_2} \nabla u \nabla v & + & \int_{\gamma} \lambda (v_1-v_2) & = & \int_\Omega f v & \forall \, v \in V \\ \int_{\gamma} \mu (u_1-u_2) & & & = & 0 & \forall \, \mu \in Q \end{array} \end{split}\]

The continuity of \(u\) is enforced by the second equation. The continuity of the normal derivative is obtained weakly from the first equation: Choosing a arbitrary test function \(v_1\), and setting \(v_2 = 0\) we obtain

\[ \int_{\Omega_1} (-\Delta u) v + \int_{\partial \Omega_1} \frac{\partial u}{\partial n} v + \int_\gamma \lambda v = \int_{\Omega_1} f v \]

Which implies \(\lambda = - \frac{\partial u}{\partial n_1}\). The same we get from \(v_2\), i.e. \(\lambda = -\frac{\partial u_2}{\partial n_2}\). The opposite sign comes from opposite outer normal vectors \(n_1 = -n_2\).

Imagine a rotating part of machine:

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

square = MoveTo(0,0).Rectangle(1,1).Face()
circo = Circle((0.5,0.5), 0.3).Face()
circ = Circle((0.5,0.5), 0.3).Face()
bar = MoveTo(0.3,0.45).Rectangle(0.4,0.1).Face()

square.edges.name="outer"
circ.edges.name="gammai"
circo.edges.name="gammao"
outer = square-circo
outer.faces.name = "outer"

circ.faces.name = "inner"
bar.faces.name = "bar"
inner = circ-bar

both = Compound([outer, inner, bar])
mesh = Mesh(OCCGeometry(both, dim=2).GenerateMesh(maxh=0.05)).Curve(3)
print (mesh.GetMaterials(), mesh.GetBoundaries())
Draw (mesh);

('outer', 'inner', 'bar') ('outer', 'outer', 'outer', 'outer', 'gammao', 'gammai', 'default', 'default', 'default', 'default')

We can virtually rotate the inner domain by setting a deformation function:

def MeshRotation(angle):
    mesh.UnsetDeformation()
    deform = GridFunction(VectorH1(mesh, order=3))

    rotmat = CF( (cos(angle), -sin(angle), sin(angle), cos(angle))).Reshape( (2,2))
    center = CF( (0.5, 0.5) )
    pos = CF( (x,y) )

    deform.Set( (rotmat-Id(2))*(pos-center), definedon=mesh.Materials("inner|bar"))
    return deform

from time import sleep
scene = Draw (mesh)

for i in range(30):
    mesh.SetDeformation(MeshRotation(i/30))
    scene.Redraw()
    sleep(0.03)

Forming the integrals of \(\lambda\) defined on one side, against an \(v\) from the other side is tricky. This is a very recent feature in NGSolve called ContactBoundary:

mesh.SetDeformation(MeshRotation(0.8))

fesu = H1(mesh, order=3, dirichlet="outer")
feslam = H1(mesh, order=3, definedon=mesh.Boundaries("gammao"))
fes = fesu*feslam

(u,lam),(v,mu) = fes.TnT()

a = BilinearForm(grad(u)*grad(v)*dx)

contact = ContactBoundary(mesh.Boundaries("gammao"), mesh.Boundaries("gammai"), volume=False)
contact.AddIntegrator (mu*(u-u.Other()) + lam*(v-v.Other()))
contact.Update (bf=a, intorder=20)

a.Assemble()
mesh.SetDeformation(MeshRotation(0))
f = LinearForm(1e3*(x-0.5)*v*dx("bar")).Assemble()
mesh.SetDeformation(MeshRotation(0.8))

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
Draw (gfu.components[0]);

Draw (grad(gfu.components[0]), mesh);