This page was generated from unit-2.9-fourthorder/fourthorder.ipynb.
2.9 Fourth order equations¶
We consider the Kirchhoff plate equation: Find
A conforming method requires
We have the following two alternatives:
Hybridized -continuous interior penalty method:¶
A simple way out is to use continuous elements, and treat the missing
[Baker 77, Brenner Gudi Sung, 2010]
We consider its hybrid DG version, where the normal derivative is a new, facet-based variable:
The facet variable is the normal derivative
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh (unit_square.GenerateMesh(maxh=0.1))
[2]:
order = 3
V1 = H1(mesh, order=order, dirichlet="left|bottom|right|top")
V2 = NormalFacetFESpace(mesh, order=order-1, dirichlet="left|bottom|right|top")
V = V1*V2
w,what = V.TrialFunction()
v,vhat = V.TestFunction()
Some proxy-functions and gridfunctions provide additional differential operators. We can get them via the Operator function. w.Operator("hesse") provides the Hessian, a matrix-valued function. Note that we can use the InnerProduct(.,.) for
[3]:
n = specialcf.normal(2)
h = specialcf.mesh_size
def jumpdn(v,vhat):
return n*(grad(v)-vhat)
def hesse(v):
return v.Operator("hesse")
def hessenn(v):
return InnerProduct(n, hesse(v)*n)
dS = dx(element_boundary=True)
a = BilinearForm(V)
a += InnerProduct (hesse(w), hesse(v)) * dx \
- hessenn(w) * jumpdn(v,vhat) * dS \
- hessenn(v) * jumpdn(w,what) * dS \
+ 3*order*order/h * jumpdn(w,what) * jumpdn(v,vhat) * dS
a.Assemble()
f = LinearForm(1*v*dx).Assemble()
[4]:
u = GridFunction(V)
u.vec.data = a.mat.Inverse(V.FreeDofs()) * f.vec
Draw (u.components[0], mesh, "disp_DG")
Draw (grad (u.components[0]), mesh, "grad")
Draw (hesse (u.components[0]), mesh, "hesse");
The Hellan-Herrmann-Johnson Mixed Method:¶
[Hellan 67, Herrmann 67, Johnson 73, Arnold+Brezzi 85, Comodi 89]
We introduce the bending moments as a new, matrix valued variable
Integration by parts in upper-right and lower-left terms:
The Hellan-Herrmann-Johnson method uses finite elements
Since
The Hellan-Herrmann-Johnson method does not require access to the neighbour element as in DG methods. In [Pechstein-Schöberl 11] the function space
[5]:
order = 3
V = HDivDiv(mesh, order=order-1)
Q = H1(mesh, order=order, dirichlet="left|bottom|top|right")
X = V*Q
print ("ndof-V:", V.ndof, ", ndof-Q:", Q.ndof)
sigma, w = X.TrialFunction()
tau, v = X.TestFunction()
n = specialcf.normal(2)
def tang(u): return u-(u*n)*n
a = BilinearForm(X, symmetric=True)
a += (InnerProduct (sigma, tau) + div(sigma)*grad(v) \
+ div(tau)*grad(w) - 1e-10*w*v )*dx \
+ (-(sigma*n) * tang(grad(v)) - (tau*n)*tang(grad(w)))*dx(element_boundary=True)
a.Assemble()
f = LinearForm(-1 * v * dx).Assemble()
gfu = GridFunction(X)
gfu.vec.data = a.mat.Inverse(X.FreeDofs()) * f.vec
Draw (gfu.components[0], mesh, name="sigma")
Draw (gfu.components[1], mesh, name="disp");
ndof-V: 3192 , ndof-Q: 1105
[ ]: