import netgen.occ as occ
import netgen.csg as csg
import netgen as ng
from netgen.meshing import Mesh

import numpy as np
import ngsolve as ngs




def generate_geometry():
    domain = occ.Box(occ.Pnt(-2, -2, -2), occ.Pnt(2,2,2))
    domain.faces.Min(occ.Z).name = "bc_low"
    domain.faces.Max(occ.Z).name = "bc_high"
    dielectric_sphere = occ.Sphere( occ.Pnt(0,0,0), 1).mat("dielectric")
    domain = occ.Glue([ domain, dielectric_sphere])
    domain.WriteStep("testcase.step")
    return domain


def solve(domain_mesh):

    fes = ngs.H1(domain_mesh, order=2, dirichlet="bc_low|bc_high")
    phi = ngs.GridFunction(fes)
    bcf = domain_mesh.BoundaryCF({"bc_low": -10, "bc_high": 10})
    phi.Set(bcf,ngs.BND)

    mat_dict = {"dielectric": 1.0, "air": 1.0}
    eps_r = domain_mesh.MaterialCF(mat_dict, default=1.0)

    u, v = fes.TnT()
    a = ngs.BilinearForm(fes, symmetric=True)
    a += ngs.grad(u)*ngs.grad(v)*ngs.dx
    a.Assemble()

    r = - a.mat * phi.vec
    phi.vec.data += a.mat.Inverse(freedofs=fes.FreeDofs()) * r
    vtk = ngs.VTKOutput(domain_mesh, coefs=[phi], names=["phi"], subdivision=1, filename="testcase")
    vtk.Do()
    return phi