from ngsolve import *

# # If you want to use the mesh from the ngspy meshing exercise, use:
#
# from mygeom import MakeMesh
# mesh = MakeMesh()
#
# # Else, import prebuilt mesh by:

mesh = Mesh("circleplate.vol.gz")


V = H1(mesh, order=2, dirichlet=[1,2])  # Approximation space 

epsl = DomainConstantCF([1, 10])        # Dielectric coefficient

a = BilinearForm(V, symmetric=True)     # Bilinear form
a+= Laplace(epsl)


extendone = DomainConstantCF([1,0])     
u1 = GridFunction(V, name="extension")
u1.Set(extendone, BND)                  # Boundary sources
Draw(u1)                                # Put extension in visual menu

u = GridFunction(V, name="solution")
u.vec.data = u1.vec

a.Assemble()

f = u.vec.CreateVector()
f.data = a.mat * u1.vec
u.vec.data -= a.mat.Inverse(V.FreeDofs(), inverse="sparsecholesky") * f

Draw(u)
Draw(-u.Deriv(), mesh, "Flux")