Boundary conditions or other problems

#Import needed packages and programs import netgen.gui from netgen.geom2d import SplineGeometry from ngsolve import * import matplotlib.pyplot as plt import netgen.geom2d as geom2d from netgen.geom2d import unit_square import numpy as np import scipy.sparse as scs from scipy.sparse.linalg import eigs, eigsh #Create geometries, first circle is the scatterer, second begins the PML third truncates the infinite domain #third circle (largest) gets a boundary condition due to the sommerfeld radiation condition. geo = SplineGeometry() geo.AddCircle( (0.0, 0.0), r=1.0, leftdomain=1, rightdomain=2, bc = "scatterer") geo.AddCircle ( (0.0, 0.0), r=1.5, leftdomain=2, rightdomain=3) geo.AddCircle ( (0.0, 0.0), r=2.0, leftdomain=3, rightdomain=0,bc = "outerbnd") #Set material 1 to be inner do that only those nodes inside the scatterer recieve the $n*u*v*dx("inner")$ term geo.SetMaterial(1, "inner") geo.SetMaterial(3, "PML") #Generate the mesh ngmesh = geo.GenerateMesh(maxh=0.04) mesh = Mesh(ngmesh) #Tell the program that we are using PML at the second circle mesh.SetPML(pml.Radial(rad=1.5,alpha=1j,origin=(0,0)), "PML") #Set the extra value of the scatterer nn=2; #I think this tells the program what level of accuracy and where our test/trial functions live fes = H1(mesh, order=4, complex=True, dirichlet="dir") #create our trial function u and our test function v u = fes.TrialFunction() v = fes.TestFunction() #create matrix K that is the non lambda term (sign for quad-eig) k = BilinearForm (fes, symmetric=True) k += grad(u)*grad(v)*dx #create matrix C that is our lambda term (sign for quad-eig) c = BilinearForm (fes, symmetric=True) c += 1.j*u*v*ds('outerbnd') #create matrix M which is our squared term and has the nn element (sign for quad-eig) m = BilinearForm (fes, symmetric=True) m += u*v*dx m += nn*u*v*dx('inner') #convert bilinearform into a matrix form k.Assemble() c.Assemble() m.Assemble() #tells the program to get ready u = GridFunction(fes, multidim=20, name='resonances') #moves the matrices from ngsolve into scipy.sparse because ngsolve hates stacking things mcsr = scs.csr_matrix(m.mat.CSR(), shape=(m.mat.height, m.mat.width),dtype='complex') ccsr = scs.csr_matrix(c.mat.CSR(), shape=(c.mat.height, c.mat.width),dtype='complex') kcsr = scs.csr_matrix(k.mat.CSR(), shape=(k.mat.height, k.mat.width),dtype='complex') #finds the size of M and creates the needed identity matrix of that size leng = m.mat.height I = scs.csr_matrix(scs.eye(leng), dtype='complex') #Builds our first block matrix in csc form B1 = scs.csr_matrix(scs.bmat([[mcsr, None], [None, I]], dtype='complex')) #Builds our second block matrix in csc form B2 = scs.csr_matrix(scs.bmat([[ccsr, kcsr],[I, None]], dtype='complex')) omegaval, omegavec = eigs(B2, k=11, M=B1, sigma=2.2-0.2*1.j) for ome in omegaval: print(ome)