Code:
from netgen.occ import *
from ngsolve import *
air = Sphere((0,0,0), 1)
air.faces.name = "air"
ground = HalfSpace((0,0,0), (0,0,-1))
ground.faces.name = "ground"
box = Box((-0.1,-0.1,-0.1), (0.1,0.1,0.1))
box.faces.name = "volt"
box.faces.maxh = 0.02
domain = (air * ground) - box
geo = OCCGeometry(domain)
mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(3)
Draw(mesh)
fes = H1(mesh, order=3, dirichlet="ground|volt")
u,v = fes.TnT()
eps = 1
a = BilinearForm(fes)
a += eps * grad(u) * grad(v) * dx
u = GridFunction(fes)
with TaskManager():
r = u.vec.CreateVector()
a.Assemble()
u.Set(1, definedon=mesh.Boundaries("volt"))
r.data = - a.mat * u.vec
u.vec.data += a.mat.Inverse(fes.FreeDofs()) * r
Draw(u, mesh, "potential")
Draw(-grad(u), mesh, "E")