Unclear functionality of specialcf.normal

Dear NGSolve developers and community,

Although there are already quite some questions regarding the specialcf.normal coefficient function, it is still unclear to me why it gives me only zero vectors in the attached example (checking e.g. “n(mesh(0,0,0))” ) and how I obtain the outer surface normal vectors (on the entire surface or only on subdomains/faces like “top” in the example) in a proper way.

Any help is greatly appreciated, thank you!

Best,
Chris

from netgen.occ import *
from ngsolve import *

cyl1 = Cylinder(Pnt(0, 0, 0), Z, r=1, h=5)
cyl1.faces.Max(Z).name = "top"
geo = OCCGeometry(cyl1)
ngmesh = geo.GenerateMesh(maxh=0.4)
mesh = Mesh(ngmesh)

n = specialcf.normal(3)

#normal = GridFunction(VectorH1(mesh, order=1))
#normal.Set(n)
#vtk = VTKOutput(ma=mesh, coefs=[normal], names=["normal"], filename="test_normalcf", subdivision=0)
#vtk.Do()

Dear Chris,

The normal vector is well-defined only on boundaries and surfaces. If you use this in the volume (as in Set) it evaluates to zero. However, if you evaluate it on the boundary (of the domain or an element) it has the expected meaningful entry. Try for example normal.Set(n,BND) and you should get a finite element function normal that has the normal field on the boundaries (averaged on edges with different normals).

Best,
Christoph

Thanks for the explanation!

And besides this handling/visualisation via GridFunctions, I can use the Coefficient Function (specialcf.normal) in bilinear forms in surface-related terms (those related to some “ds”)?

Yes, you can also point evaluate, but you need to pass the BND flag there too:

n(mesh(0,0,0, BND))