from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
from xfem import CutInfo, IF, InterpolateToP1

import numpy as np

# Background mesh
L = 1.5
max_h = 0.375
mesh_refinements = 1

background = CSGeometry()
background.Add(OrthoBrick(Pnt(-L, -L, -L), Pnt(L, L, L)))
mesh = Mesh(background.GenerateMesh(maxh=max_h, quad_dominated=False))

# Level set function
radius = 1.0
order_norm = 2
phi = (x**order_norm + y**order_norm + z**order_norm)**(1/order_norm) - radius
        
# calculating mesh refinements
refine_H1_space = H1(mesh, order=1, autoupdate=True)
refine_levelset = GridFunction(refine_H1_space, autoupdate=True)

for i in range(mesh_refinements):
    InterpolateToP1(phi, refine_levelset)
    cut_info = CutInfo(mesh, refine_levelset).GetElementsOfType(IF)
    mesh.SetRefinementFlags(cut_info)
    mesh.Refine()

global_coordinates = np.array([-0.63795623, -0.47848099, -0.61152876], dtype=np.float64, ndmin=2)
local_coordinates = mesh(global_coordinates[:,0], global_coordinates[:,1], global_coordinates[:,2])
print(local_coordinates)