JKM mixed finite elements

7. JKM mixed finite elements#

J. Gopalakrishnan, J. Guzman, J.J. Lee: “The Johnson-Krizek-Mercier elasticity element in any dimensions”

A stress finite element on Alfeld split. Stress space is symmetric and row-wise in \(H(\operatorname{div})\)

NGSolve implementation:

  • JKM space is \(nn\)-continuous

  • \(nt\) continuity by variational formulation

from ngsolve import *
from ngsolve.webgui import Draw

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

7.1. JKM - TDNNS#

Find \(\sigma \in JKM\), \(v \in {\mathcal N} \subset H(\operatorname{curl})\)

\[\begin{split} \begin{array}{ccccll} (\sigma, \tau) & + & \left< \operatorname{div} \tau, u \right> & = & 0 & \forall \, \tau \\ \left< \operatorname{div} \sigma, v \right> & & & = & f(v) & \forall \, v \end{array} \end{split}\]
fesu = HCurl(mesh, order=2, orderinner=3, type1=True, dirichlet=".*")
fesS = JKM_FESpace(mesh, order=1)
fes = fesu*fesS

(u,sigma), (v, tau) = fes.TnT()

n = specialcf.normal(2)
def tang(v): return v - (n*v)*n
irs = fesS.GetIntegrationRules(bonus_intorder=4)

def Div(sigma,v):
    return div(sigma)*v*dx(intrules=irs) - tang(v)*(sigma*n)*dx(element_boundary=True)

def Strain(stress): return Deviator(stress) + 1e-4/2 * Trace(stress) * Id(2)
def Stress(strain): return Deviator(strain) + 1e4/2 * Trace(strain) * Id(2)

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

rhs = CF( (y-1/2,0) )

term = InnerProduct(Strain(sigma),tau)*dx(intrules=irs) + Div(tau,u) + Div(sigma,v)
a = BilinearForm(term).Assemble()
f = LinearForm(-v*rhs*dx).Assemble()

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

Draw (gfu.components[0], mesh)
Draw (gfu.components[1][0,0], mesh)
gfsigma1 = gfu.components[1]

7.2. The original method#

fesu = VectorL2(mesh, order=1)
fesuf = TangentialFacetFESpace(mesh, order=1, dirichlet=".*")

fesS = JKM_FESpace(mesh, order=1)
fes = fesu*fesuf*fesS

(u,uf, sigma), (v, vf, tau) = fes.TnT()

n = specialcf.normal(2)
irs = fesS.GetIntegrationRules()

def Div(sigma,v,vf):
    return div(sigma)*v*dx(intrules=irs) - tang(vf)*(sigma*n)*dx(element_boundary=True)

term = InnerProduct(Strain(sigma),tau)*dx(intrules=irs) + Div(tau,u,uf) + Div(sigma,v,vf)
a = BilinearForm(term).Assemble()
f = LinearForm(-v*rhs*dx).Assemble()

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

Draw (gfu.components[0], mesh)
Draw (gfu.components[2][0,0], mesh)
gfsigma2 = gfu.components[2]

diff = gfsigma1-gfsigma2
Draw (diff[0,0], mesh)

WebGLScene

7.3. TDNNS#

fesu = HCurl(mesh, order=2, type1=True, dirichlet=".*")
fesS = HDivDiv(mesh, order=1, plus=True)
fes = fesu*fesS

(u,sigma), (v, tau) = fes.TnT()

n = specialcf.normal(2)
def tang(v): return v - (n*v)*n

def Div(sigma,v):
    return div(sigma)*v*dx - tang(v)*(sigma*n)*dx(element_boundary=True)

term = InnerProduct(Strain(sigma),tau)*dx(intrules=irs) + Div(tau,u) + Div(sigma,v)
a = BilinearForm(term).Assemble()
f = LinearForm(-v*rhs*dx).Assemble()

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

Draw (gfu.components[0], mesh)
Draw (gfu.components[1][0,0], mesh)
gfsigma3 = gfu.components[1]

7.4. The primal method:#

fesprimal = VectorH1(mesh, order=8, dirichlet=".*")
u,v = fesprimal.TnT()
aprimal = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v)))*dx).Assemble()
fprimal = LinearForm(v*rhs*dx).Assemble()
gfprimal = GridFunction(fesprimal)
gfprimal.vec.data = aprimal.mat.Inverse(freedofs=fesprimal.FreeDofs()) * fprimal.vec
Draw (gfprimal, mesh)
Draw (Stress(Sym(Grad(gfprimal)))[0,0], mesh)

WebGLScene

print ("err JKM-tdnns = ", Integrate ( ( Norm(gfsigma1-Stress(Sym(Grad(gfprimal))))**2), mesh))
print ("err mixed = ", Integrate ( ( Norm(gfsigma2-Stress(Sym(Grad(gfprimal))))**2), mesh))
print ("err tdnns = ", Integrate ( ( Norm(gfsigma3-Stress(Sym(Grad(gfprimal))))**2), mesh))

err JKM-tdnns =  2.2243517407269708e-06
err mixed =  2.3585620579637012e-06
err tdnns =  1.1725698768853879e-06
Contents