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.

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).

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”)?