from netgen.occ import *

from netgen.meshing import NetgenGeometry, Mesh as NGMesh

import numpy as np

class UnitSphereGeometry(NetgenGeometry):

def midpoint(self, newp, p1, p2, secpoint=0.5):

p1 = np.array([p1[0], p1[1], p1[2]])

p2 = np.array([p2[0], p2[1], p2[2]])

p = p1 + secpoint * (p2 - p1)

self.project(newp, p)

def project(self, newp, p):

pt = np.array([p[0], p[1], p[2]])

pt /= np.linalg.norm(pt)

newp[0] = pt[0]

newp[1] = pt[1]

newp[2] = pt[2]

def PointBetweenEdge(self, p1, p2, secpoint, surfi1, surfi2, ep1, ep2, newp, newgi):

self.midpoint(newp, p1, p2, secpoint)

def PointBetween(self, p1, p2, secpoint, surfi, gi1, gi2, newp, newgi):

self.midpoint(newp, p1, p2, secpoint)

def ProjectPointGI(self, surfind, p, gi):

self.project(p, p)

return True

m = OCCGeometry(Sphere(Pnt(0, 0, 0), 1)).GenerateMesh()

my_geo = UnitSphereGeometry()

m.SetGeometry(my_geo)

m.Refine()

m.Curve(3)

from ngsolve import *

mesh = Mesh(m)

Draw(mesh)