Nitsche’s Method for boundary and interface conditions

40. Nitsche’s Method for boundary and interface conditions#

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

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

A penalty approximation to the Dirichlet boundary condition \(u = u_D\) is:

\[ \int_\Omega \nabla u \nabla v + \int_{\Gamma_D} \alpha u v = \int_\Omega f v + \int_{\Gamma_D} \alpha u_D v \qquad \forall \, v \]
fes = H1(mesh, order=2)
u,v = fes.TnT()

pen = 100
a = BilinearForm(grad(u)*grad(v)*dx + pen*u*v*ds).Assemble()
f = LinearForm(10*v*dx).Assemble()

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

print ("error bc:", sqrt(Integrate((gfu-0)**2, mesh.Boundaries(".*"))))
Draw (gfu, deformation=True);

error bc: 0.05275907536407711

Exercises:

  • How does the error depend on the penalty parameter \(\alpha\) ?

  • Does the error get reduced when the approximation space gets enrichred ?

  • Modify the right hand side to set \(u = x+y\) on \(\partial \Omega\)

40.1. Nitsche’s method:#

We start from the strong form, perform integration by parts, keep the boundary term. Then, we add some more consistent terms to make the bilinear-form symmetric, and discrete coercive:

\[ \int_\Omega \nabla u \nabla v - \int_{\Gamma_D} \frac{\partial u}{\partial n} v - \int_{\Gamma_D} \frac{\partial v}{\partial n} u + \int_{\Gamma_D} \frac{\alpha p^2}{h} u v = \int_\Omega f v - \int_{\Gamma_D} \frac{\partial v}{\partial n} u_D + \int_{\Gamma_D} \frac{\alpha p^2}{h} u_D v \qquad \forall \, v \]
order = 4
alpha = 5

fes = H1(mesh, order=order)
u,v = fes.TnT()

h = specialcf.mesh_size
n = specialcf.normal(mesh.dim)

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx + alpha*order**2/h*u*v*ds
a += (-n*grad(u)*v-n*grad(v)*u) * ds(skeleton=True)
a.Assemble()

f = LinearForm(10*v*dx).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky") * f.vec

print ("error bc:", sqrt(Integrate((gfu-0)**2, mesh.Boundaries(".*"))))
Draw (gfu);

error bc: 1.9357393093205987e-05

The bilinear-form is coercive w.r.t. the norm induced by \(\int_\Omega \nabla u \nabla v dx + \int_{\partial \Omega} \frac{p^2}{h} u v ds \) if \(\alpha\) is sufficiently large.

We check positive definite by computing the few smallest eigenvalues of

\[ A x = \lambda N x, \]

where the matrix \(N\) is defined by the norm.

from ngsolve.solvers import PINVIT

bfnorm = BilinearForm(grad(u)*grad(v)*dx + order/h*u*v*ds).Assemble()
eval,evec = PINVIT(a.mat, bfnorm.mat, pre=bfnorm.mat.Inverse(), \
                   num=5, printrates=False)
print (eval)

 0.967068
 0.975969
 0.97722
 0.97906
 0.990897

Exercise:

  • Extend to non-homogeneous Dirichlet boundary conditions

  • How does the error depend on \(p\) ?

  • How does the error depend on \(\alpha\) ?

40.2. Interfaces#

Consider an electric motor with a rotating rotor. The fixed and the rotating part are meshed independently, and continuity at the interface is achieved by a Nitsche method.

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')

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)

without gluing together we compute independent solutions for both regions:

mesh.SetDeformation(MeshRotation(0.8))

fes = H1(mesh, order=3, dirichlet="outer")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx).Assemble()
f = LinearForm(1e3*(x-0.5)*v*dx("bar")).Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
Draw (gfu);

The ContactBoundary class computes integrals between two different boundaries. It integrates over the primary boundary, finds the closest point on the secondary boundary and evaluates the other function there.

deform = MeshRotation(0.8)
mesh.UnsetDeformation()
a = BilinearForm(grad(u)*grad(v)*dx + u*v*dx)

contact = ContactBoundary(mesh.Boundaries("gammai"), mesh.Boundaries("gammao"), volume=True)
h = specialcf.mesh_size

contact.AddIntegrator (3/h*(u-u.Other()) * (v-v.Other()))
# consisteny term not yet implemented for contact boundary
contact.AddIntegrator (n*grad(u)*(v.Other()-v)+n*grad(v)*(u.Other()-u))

contact.Update (deform, bf=a, intorder=20)
a.Assemble()

f = LinearForm(1e3*(x-0.5)*v*dx("bar")).Assemble()

gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs())*f.vec
mesh.SetDeformation(deform)

Draw (gfu)
Draw (grad(gfu), mesh);

The current implementation is very simple and not highly accurate. It uses Gauss-rules on the primary boundary which are accurate for finite element functions on the primary boundary, but not for the other boundary. One can observe oscillations for the gradient near to the boundary.

40.3. Hybrid Interfaces#

In a hybrid interface method one introduces another field \(\hat u\) only at the interface. The functions \(u\) from both sides are glued to this common interface field by a Nitsche method:

\[ \int_{\Omega_i} \nabla u \nabla v - \int_\gamma \partial_n u (v - \hat v) - \int_\gamma \partial_n u (u - \hat u) + \int_\gamma \frac{\alpha p^2}{h} (u - \hat u)(v - \hat v) \]

Now one has to perform integrals between boundary values of finite element functions, and the interface functions. Often the interface is geometrically simple (a circle or a cylinder), and one can choose global functions there. On the circle we choose trigonometric functions.

In NGSolve there is the GlobalInterfaceSpace. It allows to provide the coordinate (the angle) as a function of the global coordinates. This mapping can be different on both sides of the interface, which allows the gluing of shifted functions.

angle = Parameter(0.8) 
mapping = atan2(y-0.5, x-0.5) + \
    mesh.MaterialCF({"inner" : angle}, default=0)

interface = mesh.Boundaries("gammao|gammai")
print (mesh.GetMaterials())
print (mesh.GetBoundaries())
print (interface.Mask())

('outer', 'inner', 'bar')
('outer', 'outer', 'outer', 'outer', 'gammao', 'gammai', 'default', 'default', 'default', 'default')
0: 0000110000

order=3
V = H1(mesh, order=order, dirichlet="outer")
from ngsolve.comp import GlobalInterfaceSpace
Vhat = GlobalInterfaceSpace(mesh, mapping=mapping, \
                            order=10, periodic=True, definedon=mesh.Boundaries("gamma.*"))
X = V *Vhat

mesh.UnsetDeformation()

u,uhat = X.TrialFunction()
v,vhat = X.TestFunction()

h = specialcf.mesh_size
n = specialcf.normal(mesh.dim)

a = BilinearForm(X) 
a += grad(u)*grad(v)*dx 
a += (3*order**2/h*(u-uhat)*(v-vhat) + \
    n*grad(u)*(vhat-v)+n*grad(v)*(uhat-u))* \
    ds("gammao|gammai", skeleton=True)
a.Assemble();

f = LinearForm(1e3*(x-0.5)*v*dx("bar")).Assemble()
gf = GridFunction(X)

gf.vec.data = a.mat.Inverse(freedofs=X.FreeDofs())*f.vec
gfu, gfuhat = gf.components

mesh.SetDeformation(MeshRotation(angle))
Draw (gfu)
Draw (grad(gfu), mesh);
mesh.UnsetDeformation()
Contents