6.7. Cooke’s membrane

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()
Draw (shape)
mesh = shape.GenerateMesh(dim=2, maxh=5)
Draw (mesh);

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

if order>=1:
    fessigma = HDivDiv(mesh, order=order, plus=True, dirichlet="default|right")
else:
    fessigma = HDivDiv(mesh, order=1, orderinner=0, dirichlet="default|right")

fes = fesu*fessigma

(u,sigma), (du,dsigma) = fes.TnT()
n = specialcf.normal(2)

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

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

term = InnerProduct(-Strain(sigma,E,nu),dsigma)*dx + Div(dsigma,u) + Div(sigma,du)
a = BilinearForm(term).Assemble()

traction = CF((0,20))
f = LinearForm(traction*du.Trace()*ds("right"))
f.Assemble()

gf = GridFunction(fes)
gf.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
gfu, gfsigma = gf.components

Draw (gfu[1],mesh, deformation=0.001*gfu, order=5);

Draw (gfsigma[0,0], mesh, min=-50, max=50); # min=-20e6, max=20e6);

gfu(mesh(48, 60)) 

(-6042.552020446337, 8053.028572175645)

fesu.ndof, fessigma.ndof

(1302, 2135)

mesh.GetNE(VOL), mesh.nedge, mesh.nv

(119, 196, 78)

"order 0 dofs", 1*mesh.nedge, 2*mesh.nedge

('order 0 dofs', 196, 392)

"order 1 dofs", 2*mesh.nedge+2*mesh.GetNE(VOL),  2*mesh.nedge+5*mesh.GetNE(VOL)

('order 1 dofs', 630, 987)

"order 2 dofs", 3*mesh.nedge+6*mesh.GetNE(VOL),  3*mesh.nedge+13*mesh.GetNE(VOL)

('order 2 dofs', 1302, 2135)