Gradient calculation on the boundary with dirichlet bc

from ngsolve import * from netgen.geom2d import * from netgen.occ import * U = 10          # velocity of surface rho = 890       # density of the oil eta = 0.0078792 # viscosity of the oil filename_geometry = 'wedge.step' # maximum element size of mesh mesh_max_h = 1e-1 # import geometry from given path geo = OCCGeometry(filename_geometry) # generate and refine mesh mesh = Mesh(geo.GenerateMesh(maxh=mesh_max_h)) # generate finite elemente spaces V = VectorH1(mesh,order=3, dirichlet='bottom|top|NONE')   # velocity #V = VectorH1(mesh,order=3, dirichlet='inlet')   # velocity Q = H1(mesh,order=2, dirichlet='inlet|outlet|right|left')    # pressure X = V * Q # obtain test and trial functions u,p = X.TrialFunction() v,q = X.TestFunction() # space for solution data gfu = GridFunction(X) # boundary condition uin = CoefficientFunction((U,0,0)) gfu.components[0].Set(uin, definedon=mesh.Boundaries("bottom")) # draw the boundary conditions Draw (Norm(gfu.components[0]), mesh, "velocity", sd=3) nu = (eta/rho) SetNumThreads(8) with TaskManager():     # right hand side of the equation     f = LinearForm(X)     f.Assemble()     # left hand side of the equation     stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q     #stokes = (eta/rho)*InnerProduct(grad(u), grad(v)) + div(u)*q - div(v)*p + InnerProduct(v,div(u)*u) - 1e-12*p*q     a = BilinearForm(X, symmetric=True)     a += stokes*dx     a.Assemble()     # calculate solution data     inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")     res = f.vec.CreateVector()     res.data = f.vec - a.mat*gfu.vec     gfu.vec.data += inv_stokes * res # density free pressure --> pressure gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho velocity = gfu.components[0] # friction calculation c_u = GridFunction(Q) c_u.Set(velocity[0]) c_w = GridFunction(Q) c_w.Set(velocity[2]) Draw(Norm(Grad(velocity)),mesh,'velocity gradient') dudz = Grad(c_u)[2] dwdx = Grad(c_w)[0] tau = eta * (dudz + dwdx) #Draw(Norm(Grad(velocity)),mesh,'text') #Draw(tau,mesh,'dudz') F_fric = Integrate(cf = tau,mesh=mesh,VOL_or_BND=BND,definedon=mesh.Boundaries('bottom')) print('Friction force:',F_fric)