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

4

I have been trying to extend the above to 1D open line in 2D without success, any tip?
I tried:

  • using the FacetSurface, but the code crashes.
  • using Nodal spaces but it says there are no evaluators for bbnd
  • Using :
jac = specialcf.JacobianMatrix(2,1)[:,0]
xref = specialcf.xref(1)
conormal = (2*xref[0]-1)*jac

as discussed in previous questions together with different integrations on the bboundary or on the boundary done elementwise.

Thanks a lot,
Alessandro

5

Hi Alessandro,

in 2D most spaces become either H1 or L2 and so some spaces are not implemented in 1D. The specialcf.tangential takes the role of the co-normal vector for 1D domains.

The indicator function can be constructed using an H1-GridFunction, setting the corresponding vertex dofs to 1.

Best,
Michael

1 Like
6

Yes, I see, I was eventually able to figure it out but forgot to update this. Thanks a lot for the help!