Computation of Curvature

11. Computation of Curvature#

work with Michael Neunteufel, Max Wardetzky and Jay Gopalakrishnan

  • 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

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

Define Gauss curvature from Weingarten tensor:

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

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

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

11.1.2. 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.999839532080625

11.1.3. 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.0002596199782285

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

11.1.4. 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 = 1.0000000000007234

11.2. Curvatur computation from metric tensor on parameter domain#

  • Evan S. Gawlik: High-order approximation of Gaussian Curvature with Regge finite elements, SINUM ‘20

  • J. Gopalakrishnan, M. Neunteufel, JS, M. Wardetzky: Analysis of curvature approximation via covariant curl and incompatibility for Regge metrics (WIP)

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import CSG2d, Circle
from math import pi

order    = 4
addorder = 2

geo = CSG2d()
circle = Circle(center=(0,0), radius=1.0)
geo.Add(circle)

mesh = Mesh(geo.GenerateMesh(maxh=0.1))
mesh.Curve(10)
Draw(mesh);

gfset = GridFunction(H1(mesh,order=order)**3)

r = sqrt(x**2+y**2)
phi = atan2(y,x)

gfset.Interpolate( (cos(phi)*sin(r*pi/2)-x,sin(phi)*sin(r*pi/2)-y,cos(r*pi/2)))

gfsetxy = GridFunction(H1(mesh,order=order)**2)
gfsetxy.Set ( gfset[0:2] )
mesh.SetDeformation(gfsetxy)
Draw (gfset[2], mesh, deformation=True)
mesh.UnsetDeformation()

Interpolating metric

\[ C = \nabla \varphi^T \nabla \varphi \]

into a Regge finite element space:

F = CF( (1,0, 0,1, 0,0),dims=(3,2)) + Grad(gfset)
C = F.trans*F

fesCC = HCurlCurl(mesh, order=3) 
gfC = GridFunction(fesCC)
gfC.Set(C, dual=True, bonus_intorder=4)   # half-sphere
# gfC.Set(Id(2), dual=True)    # disk

Draw (gfC[0,0], mesh);

fesR = H1(mesh, order=3)
fR = LinearForm(fesR)
u,v = fesR.TnT()

charVertex = GridFunction(H1(mesh,order=1))
charVertex.Set(1)
for el in charVertex.space.Elements(BND):
    for d in el.dofs:
        charVertex.vec[d] = 1     # set to 0 to skip bnd vert

charEdge = GridFunction(FacetFESpace(mesh,order=0))  
charEdge.vec[:] = 1
for el in charEdge.space.Elements(BND):
    for d in el.dofs:
        charEdge.vec[d] = 1     # set to 0 to skip bnd edge
       
# vertex - contributions = angle deficit
vertextang  = specialcf.VertexTangentialVectors(2)
vt1 = vertextang[:,0]
vt2 = vertextang[:,1]

fR += charVertex*v*acos(gfC[vt1,vt2]/sqrt(gfC[vt1,vt1]*gfC[vt2,vt2]))*dx(element_vb=BBND)
fR += -charVertex*v*acos((vt1*vt2)/sqrt((vt1*vt1)*(vt2*vt2)))*dx(element_vb=BBND)

t  = specialcf.tangential(2,True)
n = specialcf.normal(2)  
edgecurve = specialcf.EdgeCurvature(2)
Gamma_ttn = gfC.Operator("christoffel2")[t,t,:]
fR += charEdge*v*sqrt(Det(gfC))/gfC[t,t]*n*(Gamma_ttn+edgecurve)*dx(element_vb=BND, bonus_intorder=4)


Riemann = gfC.Operator("Riemann")
fR += 1/sqrt(Det(gfC))*Riemann[0,1,0,1] * v * dx(bonus_intorder=4)

fR.Assemble()

gfR = GridFunction(fesR)
mass = BilinearForm(sqrt(Det(gfC))*u*v*dx(bonus_intorder=4)).Assemble().mat
gfR.vec.data = mass.Inverse() * fR.vec
Draw (gfR, min=0.99, max=1.01);

print ("total curvature on half-sphere:", Integrate (gfR, mesh) / pi)

total curvature on half-sphere: -1.0000000000155922

11.3. Current work: numerics for Einstein equation#

  • solve for an evolving metric in 3+1

  • can do now linearized Einstein wave equation

Contents