Hi Michael

Thanks for your reply. I did what you mentioned above, but it returns 0 when I print instead of -5.

Here is my code:

from netgen.occ import *

import ngsolve

from ngsolve import *

from ngsolve.webgui import Draw

from ngsolve.krylovspace import CGSolver

box = Box((0,0,0), (3,0.6,1))

geo = box

geo.faces.Min(X).name = “fix”

geo.faces.Max(X).name = “force”

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)

E, nu = 210, 0.2

mu = E / 2 / (1+nu)

lam = E * nu / ((1+nu)*(1-2*nu))

def Stress(strain):

return 2*mu*strain + lam*Trace(strain)*Id(3)

fes = VectorH1(mesh, order=3, dirichlet=“fix”)

u,v = fes.TnT()

gfu = GridFunction(fes)

with TaskManager():

a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)*

pre = Preconditioner(a, “bddc”)

a.Assemble()

force = CF((-5,0,0))

f = LinearForm(forcev*ds(“force”)).Assemble()

inv = CGSolver(a.mat, pre, printrates=True, tol=1e-6)

gfu.vec.data = inv * f.vec

with TaskManager():

fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)

gfstress = GridFunction(fesstress)

gfstress.Interpolate (Stress(Sym(Grad(gfu))))

signed_fesstress = H1(mesh,order=3)

gf_signed_stress = GridFunction(signed_fesstress)

n = specialcf.normal(3)

gf_signed_stress.Set(Stress(Sym(Grad(gfu)))*n *n, BND)

gf_none_signed_stress = Norm(gfstress)

print(gf_signed_stress(mesh(3,0.3,0.5)))

print(gf_none_signed_stress(mesh(3,0.3,0.5)))

Can you please tell me what I am doing wrong?

Regards

Hamed