Co-normal to embedded open surface boundary

General
1

Hi,

I have a problem where I need to use the co-normal to a surface boundary in order to find a weak version of the (vector) mean curvature. An example code is given below

from ngsolve import *
from netgen.csg import *
from netgen.meshing import MeshingStep
from ngsolve.webgui import Draw

geo          = CSGeometry()
sphere       = Sphere(Pnt(0,0,0), 1)
bot          = Plane(Pnt(0,0,0), Vec(0,0,-1))
finitesphere = sphere * bot

geo.AddSurface(sphere, finitesphere.bc("surface"))
geo.NameEdge(sphere,bot, "bottom")

mesh = Mesh(geo.GenerateMesh(maxh=0.1))

V1 = VectorH1(mesh, order=1, definedon=mesh.Boundaries("surface"))

nu = specialcf.normal(mesh.dim)
Ps = Id(mesh.dim) - OuterProduct(nu, nu)

tang = specialcf.tangential(mesh.dim)
xsi = CF((0, 0, -1))
# xsi = CF((0, 0, -sqrt(tang[0]**2+tang[1]**2))) also works

(kappa, eta) = V1.TnT()
M = BilinearForm(V1)
l = LinearForm(V1)
M += kappa*eta*ds
M.Assemble()
l += -InnerProduct(Ps, Grad(eta).Trace())*ds
l += InnerProduct(xsi, eta)*ds(definedon=mesh.BBoundaries("bottom"))
l.Assemble()

kappa_h = GridFunction(V1)
kappa_h.vec.data = M.mat.Inverse(V1.FreeDofs())*l.vec

Draw(kappa_h)

I tried to use

xsi = Cross(tang, nu)

but the result is a zero Coefficient Function. What should I do in order to get that co-normal working no matter the shape of the boundary? Is there some specific CF I am not aware? Or should I use other CFs like specialcf.JacobianMatrix(3,) and special.xref() to compute it manually by rotating tang?

Thanks a lot for any help.
Best,
Alessandro

2

Hi Alessandro,

The surface’s normal vector is not well-defined on the co-dimension two boundary. It is therefore zero.

You can use an indicator function being 1 only on the specific edges

gfF = GridFunction(FacetSurface(mesh, order=0))
gfF.Set(1, definedon=mesh.BBoundaries("bottom"))

and use it in combination with element_boundary=True. It loops over all surface elements, and for each element, it loops over its boundaries

xsi = Cross(nu, tang)

l += InnerProduct(xsi, eta) * gfF * ds(element_boundary=True)

Please note, that the edges will be evaluated from both sides, computing a ‘‘jump’’. If you want it only from one boundary you can use an additional indicator function via the L2Surface space.

Best,
Michael

1 Like
3

Thanks a lot for the quick answer, works perfectly!

Best,
Alessandro