import numpy
import netgen.meshing as ngm
from netgen.webgui import Draw

dim = 3
#backend = "spline"
#backend = "csg"
backend = "occ"

isoccgeom = (backend == "occ")

def mesh(dim, geo_type):
    maxh = 0.75
    if dim == 2:
        if geo_type == "occ":
            from netgen.occ import Circle, OCCGeometry
            circle = Circle((0, 0), 1.0).Face()
            circle.edges.name = "surface"
            geo = OCCGeometry(circle, dim=2)
        elif geo_type == "spline":
            from netgen.geom2d import SplineGeometry
            geo = SplineGeometry()
            geo.AddCircle(c=(0, 0), r=1.0, bc="surface")
        elif geo_type == "csg":
            from netgen.geom2d import CSG2d, Circle
            geo = CSG2d()
            geo.Add(Circle(center=(0, 0), radius=1, bc="surface"))
        else:
            raise ValueError(f"Unexpected geometry backend {geo_type}")
    elif dim == 3:
        if geo_type == "occ":
            from netgen.occ import Sphere, OCCGeometry
            sphere = Sphere((0, 0, 0), 1.0)
            sphere.faces.name = "surface"
            geo = OCCGeometry(sphere, dim=3)
        elif geo_type == "csg":
            from netgen.csg import CSGeometry, Sphere, Pnt
            geo = CSGeometry()
            sphere = Sphere(Pnt(0, 0, 0), 1)
            sphere.bc("surface")
            geo.Add(sphere)
            maxh = 0.5
        else:
            raise ValueError(f"Unexpected geometry backend {geo_type}")
    else:
        raise ValueError(f"Unexpected dimension {dim}")
    ngmesh = geo.GenerateMesh(maxh=maxh)
    return ngmesh
    

ngmesh = mesh(dim, backend)
dim = ngmesh.dim
geo = ngmesh.GetGeometry()

# This refines the mesh without projecting the points to the boundary
ngmesh.SetGeometry(None)
ngmesh.Refine()
coordinates = ngmesh.Coordinates()

# Construct a new mesh and project the points on the boundary
rngmesh = ngm.Mesh(dim=dim)
rngmesh.SetGeometry(geo)
rngmesh.AddPoints(coordinates)
                
def addSimplices(ngMesh, dim, index, data, project_geometry, isoccgeom):
    if len(data) == 0:
        return
    if dim == 1:
        d = ngm.EdgeDescriptor()
        d.index = index
        d.edgenr = index
        index = ngMesh.Add(d)
    elif dim == 2:
        surfnr = index if isoccgeom else index-1
        index = ngMesh.Add(ngm.FaceDescriptor(bc=index, surfnr=surfnr))
    base = 1
    ngMesh.AddElements(dim=dim, index=index, data=data, base=base,
                       project_geometry=project_geometry)


els = [ngmesh.Elements0D, ngmesh.Elements1D, ngmesh.Elements2D, ngmesh.Elements3D]

start = 1       # fails with occ 3D
#start = dim-1

for d in range(start, dim+1):
    els_np = els[d]().NumPy()
    index = els_np["index"]
    
    for i in set(index):
        nodes = els_np["nodes"][index == i, :d+1]
        addSimplices(rngmesh, d, i, nodes, True, isoccgeom)
        

rngmesh.Save("disk.vol")
Draw(rngmesh, filename="mesh.html")

rngmesh = ngmesh
rels = [rngmesh.Elements0D, rngmesh.Elements1D, rngmesh.Elements2D, rngmesh.Elements3D]

rels_np = rels[dim-1]().NumPy()
index = rels_np["index"]
nodes = rels_np["nodes"]
#print(numpy.linalg.norm(ngmesh.Coordinates()[nodes[index != 0, :dim]], axis=1))