%matplotlib notebook
import netgen.gui
from ngsolve import *
import matplotlib.pyplot as plt
from netgen.geom2d import CSG2d, Circle, Rectangle, Solid2d

geo = CSG2d()

# define some primitives
inner_circle = Circle(center=(0,0), radius=1, mat="Omega", bc="omega_bc")
tri = Solid2d( [(0.5,0.4), (0.1,-0.2), (0.9,-0.2)], mat="mat2", bc="bc_tri" )
rect1 = Rectangle( pmin=(-0.5,-0.5), pmax=(-0.1,-0.1), mat="mat2", bc="bc_rect1")
rect2 = Rectangle( pmin=(-0.1,0.1), pmax=(-0.5,0.5), mat="mat2", bc="bc_rect2")
pml_circle = Circle(center=(0,0), radius=2.0, mat="pml", bc="pml_bc")


# use operators +, - and * for union, difference and intersection operations
domain1 = rect1 + rect2 + tri
domain1.Mat("solid").Maxh(0.03)
#domain1.Maxh(0.01) # change domain name and maxh
domain2 = inner_circle - (rect1 + rect2 + tri)
domain2.Mat("Omega1")
domain3 = pml_circle - inner_circle
domain3.Mat("pml")
#domain4 = inner_circle
#domain4.Mat("Omega")

# add top level objects to geometry
geo.Add(domain1)
geo.Add(domain2)
geo.Add(domain3)
#geo.Add(domain4)

mesh = Mesh(geo.GenerateMesh (maxh=0.1))

mesh.SetPML(pml.Radial(rad=1.75,alpha=1j,origin=(0,0)), "pml") #Alpha is the strenth of the PML.

mesh.Curve(3)

import numpy as np

omega_0 = 20
omega_tilde = 19.5

domain_values = {'solid': omega_tilde,'Omega': omega_0, 'Omega1': omega_0, 'pml': omega_0}
values_list = [domain_values[mat] for mat in mesh.GetMaterials()]
omega = CoefficientFunction(values_list)

fes = H1(mesh, complex=True, order=5)

s = 3/2*np.pi

u_in =exp(1j*omega_0*(cos(s)*x + sin(s)*y)) #Can use any vector as long as it is on the unit circle.

#Defining our test and solution functions.
u = fes.TrialFunction()
v = fes.TestFunction()

#Defining our LHS of the problem.
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx - omega**2*u*v*dx
a += -omega*1j*u*v * ds("omega_bc")
a.Assemble()

#Defining the RHS of our problem.
f = LinearForm(fes)
f += u_in * (omega**2 - omega_0**2) * v * dx
f.Assemble()

#Solving our problem.
u_s = GridFunction(fes, name="u")
u_s.vec.data = a.mat.Inverse() * f.vec

u_tot = u_in + u_s

Draw(u_in, mesh, "u_in")
Draw(u_s, mesh, "u_s")
Draw(u_tot, mesh, "u_tot")

def f1_field(r,s):
   
    temp = Integrate(exp(-1j*omega*((cos(r)-cos(s))*x+(sin(r)-sin(s))*y))*(omega**2 - omega_0**2),
              mesh,definedon=mesh.Materials("Omega1"))
   
    return (((1/(8*omega_0*np.pi))**0.5)*exp(1j*np.pi/4))*temp


r = 1/4 * np.pi   
   
point = mesh(0.8,0.6)

print(n(point))

print("First Order:", f1_field(r,s))