Computation of Curvature

34. Computation of Curvature#

work with Michael Neunteufel, Max Wardetzky and Jay Gopalakrishnan reference

  • Given metric tensor \(g \in H(cc)\), compute Riemannian curvature \(R(g)\)

  • Compute curvature of a triangulated manifold

from netgen.occ import *
from netgen import meshing
from ngsolve import *
from ngsolve.webgui import Draw

MESHSURF = meshing.MeshingStep.MESHSURFACE

34.1. Curvature on smooth manifolds#

sp = Sphere(Pnt(0,0,0), r=1)
spmesh = Mesh(OCCGeometry(sp).GenerateMesh(maxh=0.2,perfstepsend=MESHSURF))
spmesh.Curve(3);

Weingarten tensor (a symmetric matrix):

\[ W = \nabla_t n \]

Scalar curvature (Gauss curvature) is product of non-zero eigenvalues of \(W\).

nsurf = specialcf.normal(3)
GaussCurv = Cof(Grad(nsurf))*nsurf*nsurf

Draw (GaussCurv, spmesh, min=0.9, max=1.1);

By Gauss-Bonnet theorem: \(\int_S \kappa = 4 \pi\)

print ("integrated curvature:",
       Integrate (GaussCurv, spmesh.Boundaries(".*")) / (4*pi))

integrated curvature: 1.0121953978547735

34.2. A cylinder has only flat faces:#

cyl = Cylinder( (0,0,0), Z, r=1, h = 2)
cylmesh = Mesh(OCCGeometry(cyl).GenerateMesh(maxh=0.2,perfstepsend=MESHSURF)).Curve(5)
Draw (GaussCurv, cylmesh, min=-0.0001, max=0.0001);

cyl = Cylinder( (0,0,0), Z, r=1, h = 2)
cylf = cyl.MakeFillet(cyl.edges, 0.1)
cylfmesh = Mesh(OCCGeometry(cylf).GenerateMesh(maxh=0.2,perfstepsend=MESHSURF)).Curve(7)
Draw (GaussCurv, cylfmesh);

Integrate (GaussCurv, cylfmesh.Boundaries(".*"), order=10) / (4*pi)

1.000139999505865

34.3. Edge term: jump of geodesic curvature#

\[ \kappa_g = \nu \cdot \nabla_t t \]

with co-normal vector \(\nu = n_F \times t_E\)

nu = Cross(specialcf.normal(3), specialcf.tangential(3)) 
edgecurve = specialcf.EdgeCurvature(3)  # nabla_t t
GeodesicCurv = nu*edgecurve

fes = SurfaceL2(cylmesh, order=5)
u,v = fes.TnT()
f = LinearForm(fes)
f += GaussCurv*v*ds(bonus_intorder=8)
f += -GeodesicCurv*v*ds(element_boundary=True, bonus_intorder=8)
f.Assemble()
gfu = GridFunction(fes)
mass = BilinearForm(u*v*ds).Assemble().mat
gfu.vec.data = mass.Inverse()*f.vec
print ("integrated curv =", Integrate(gfu, cylmesh.Boundaries(".*"))/(4*pi))
Draw(gfu);

integrated curv = 0.9998395320792809

34.4. Vertex term: angle deficit#

cube = Box( (0,0,0), (1,1,1) )
cubef = cube.MakeFillet(cube.edges, 0.1)
cubefmesh = Mesh(OCCGeometry(cubef).GenerateMesh(maxh=0.2,perfstepsend=meshing.MeshingStep.MESHSURFACE)).Curve(7)
print ("integrated curv =", Integrate(GaussCurv, cubefmesh.Boundaries(".*")) / (4*pi))
Draw (GaussCurv, cubefmesh);

integrated curv = 1.0002596199781306

cube = Box( (0,0,0), (1,1,1) )
cubemesh = Mesh(OCCGeometry(cube).GenerateMesh(maxh=0.2,perfstepsend=meshing.MeshingStep.MESHSURFACE)).Curve(7)

fes = H1(cubemesh, order=5)
u,v = fes.TnT()

f = LinearForm(fes)
bbndtang  = specialcf.VertexTangentialVectors(3)
bbndtang1 = bbndtang[:,0] 
bbndtang2 = bbndtang[:,1] 
f += -v*acos(bbndtang1*bbndtang2)*ds(element_vb=BBND)

f.Assemble()
for i in range(f.space.mesh.nv):
    f.vec[i] += 2*pi

gfu = GridFunction(fes)
mass = BilinearForm(u*v*ds).Assemble().mat
gfu.vec.data = mass.Inverse()*f.vec
print ("integrated curv =", Integrate(gfu, cubemesh.Boundaries(".*"))/(4*pi))
Draw(gfu);

integrated curv = 0.9999999999999886

34.5. All contributions together#

sp = Sphere(Pnt(0,0,0), 1)
spmesh = Mesh(OCCGeometry(sp).GenerateMesh(maxh=0.3,perfstepsend=MESHSURF))
spmesh.Curve(2)
Draw (spmesh);

fes = H1(spmesh, order=2)
u,v = fes.TnT()
f = LinearForm(fes)
f += GaussCurv*v*ds(bonus_intorder=12)
f += -GeodesicCurv*v*ds(element_boundary=True, bonus_intorder=12)
f += -v*acos(bbndtang1*bbndtang2)*ds(element_vb=BBND)

f.Assemble()
for i in range(f.space.mesh.nv):
    f.vec[i] += 2*pi

gfu = GridFunction(fes)
mass = BilinearForm(u*v*ds(bonus_intorder=6)).Assemble().mat
gfu.vec.data = mass.Inverse()*f.vec
print ("integrated curv =", Integrate(gfu, spmesh.Boundaries(".*"), order=15)/(4*pi))
Draw(gfu);

integrated curv = 0.9999999999998085
Contents