6.8. Hybrid TDNNS for Cooke’s membrane benchmark

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

shape = MoveTo(0,0).Line(48,44).Line(0,16,name="right").Line(-48,-16).Close(name="fix").Face()
shape.vertices.hpref=1
mesh = shape.GenerateMesh(dim=2, maxh=10)
# mesh.RefineHP(4, 0.4)
Draw (mesh);

E,nu = 1, 1/3

# plane stress:
def Strain(stress, E,nu):
    return (1+nu)/E * stress - nu/E*Trace(stress)*Id(2)

def Stress(strain, E, nu):
    return E/(1+nu) * strain + (E*nu)/(1-nu**2) * Trace(strain)*Id(2)
    
traction = CF((0,20))    

def SolveTDNNS(order=1):

    fesu = HCurl(mesh, order=order, type1=True, dirichlet="fix")

    if order == 1:
        fesuhat = NormalFacetFESpace(mesh, order=1, dirichlet="fix")
        fessigma = HDivDiv(mesh, order=1, orderinner=0)
    else:
        fesuhat = NormalFacetFESpace(mesh, order=order-1, dirichlet="fix")
        fessigma = HDivDiv(mesh, order=order-1, plus=True)

    fessigma = Discontinuous(fessigma)
    # fessigma = PrivateSpace(Discontinuous(fessigma))

    fes = fesu*fesuhat*fessigma
    (u,uhat, sigma), (du,duhat, dsigma) = fes.TnT()
    
    n = specialcf.normal(2)

    def Div(sigma,v, vhat):
        return InnerProduct(-Grad(v),sigma)*dx + InnerProduct( (v-vhat)*n,sigma[n,n])*dx(element_boundary=True)

    term = InnerProduct(-Strain(sigma,E,nu),dsigma)*dx - Div(dsigma,u,uhat) - Div(sigma,du,duhat)
    a = BilinearForm(term, eliminate_hidden=True, condense=True).Assemble()
    f = LinearForm(traction*du.Trace()*ds("right")).Assemble()

    gf = GridFunction(fes)
    # print (a.mat.nze)
    Solve (a*gf==f)
    return gf

def SolveStd (order=1):
    fes = VectorH1(mesh, order=order, dirichlet="fix")
    u,v = fes.TnT()
    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u)), E,nu), Grad(v))*dx, condense=True).Assemble()
    # print (a.mat.nze)

    f = LinearForm(traction*v*ds("right")).Assemble()
    gf = GridFunction(fes)
    Solve (a*gf==f)

    return gf

gf = SolveTDNNS(order=3)
gfu, gfuhat, gfsigma = gf.components

Draw (gfu[1],mesh, deformation=0.001*gfu, order=5);
Draw (gfsigma[0,0], mesh, min=-50, max=50, order=5); # min=-20e6, max=20e6);

gfu = SolveStd(order=3)

Draw (gfu[1],mesh, deformation=0.001*gfu, order=5);
Draw (Stress(Sym(Grad(gfu)),E,nu)[0,0], mesh, min=-50, max=50, order=5); 

orders = []
uy_std = [] 
uy_TDNNS = []

for order in range(1, 10):
    orders.append(order)
    gfstd = SolveStd(order=order)
    uy_std.append (gfstd[1] (mesh(48, 60)))
    gfu, gfuhat, gfsigma = SolveTDNNS(order=order).components
    uy_TDNNS.append (gfu[1] (mesh(48, 60)))

print (orders, uy_std, uy_TDNNS)

[1, 2, 3, 4, 5, 6, 7, 8, 9] [6432.58487326342, 7890.354822127499, 7987.699680338658, 8018.658168676744, 8032.963804033883, 8041.068827862933, 8046.07988598988, 8049.38623807576, 8051.681683714017] [7591.372436476037, 8072.569671083357, 8047.463513991839, 8047.604961788009, 8052.701496056395, 8056.0015873251905, 8058.046137237404, 8059.376461293968, 8060.26877482604]

import matplotlib.pyplot as plt

plt.figure(figsize=(6.5, 4.0))

plt.plot(orders, uy_std, marker="o", linestyle="-", linewidth=1.8, label="std-FEM")
plt.plot(orders, uy_TDNNS, marker="s", linestyle="--", linewidth=1.8, label="TDNNS")

plt.xlabel("order")
plt.ylabel(r"$u_y$")
plt.title(r"$u_y$ vs order")
plt.grid(True, which="both", linestyle=":", linewidth=0.8)
plt.legend()
plt.tight_layout()
plt.show()

refval = 8060.605428255632   # computed with strong refinement to corners
errstd = [abs(uy-refval) for uy in uy_std]
errTDNNS = [abs(uy-refval) for uy in uy_TDNNS]

import matplotlib.pyplot as plt
plt.figure(figsize=(6.5, 4.0))

plt.plot(orders, errstd, marker="o", linestyle="-", linewidth=1.8, label="std-FEM")
plt.plot(orders, errTDNNS, marker="s", linestyle="--", linewidth=1.8, label="TDNNS")
plt.yscale("log")

plt.xlabel("order")
plt.ylabel(r"$u_y$")
plt.title(r"$u_y$ vs order")
plt.grid(True, which="both", linestyle=":", linewidth=0.8)
plt.legend()
plt.tight_layout()
plt.show()

uy_TDNNS[-1], uy_std[-1]

(8060.26877482604, 8051.681683714017)