So i'm running some code to solve a quadratic eigenvalue problem but encountering an error when I run it.
Code:
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
geo = SplineGeometry()
geo.AddCircle( (0.0, 0.0), r=1.0, leftdomain=0, rightdomain=1, bc = "scatterer")
geo.AddCircle ( (0.0, 0.0), r=1.4, leftdomain=1, rightdomain=2)
geo.AddCircle ( (0.0, 0.0), r=1.9, leftdomain=2, rightdomain=0)
ngmesh = geo.GenerateMesh(maxh=0.02)
mesh = Mesh(ngmesh)
mesh.SetPML(pml.Radial(origin=(0.0,0.0), rad=1.4, alpha=1j), definedon=2)
nn=2;
fes = H1(mesh, order=4, complex=True, dirichlet="dir")
u = fes.TrialFunction()
v = fes.TestFunction()
k = BilinearForm (fes, symmetric=True)
k += -1*grad(u)*grad(v)*dx
c = BilinearForm (fes, symmetric=True)
c += 1j*u*v*ds("outerbnd")
m = BilinearForm (fes, symmetric=True)
m += (1+nn)*u*v*dx
k.Assemble()
c.Assemble()
m.Assemble()
I = IdentityMatrix(fes.ndof, complex=False)
II = -1*I
B1 = BlockMatrix([[m.mat, None], [None, I]])
B2 = BlockMatrix([[c.mat, k.mat], [II, None]])
u = GridFunction(fes, multidim=20, name='resonances')
with TaskManager():
lam = ArnoldiSolver(B2, B1, fes.FreeDofs(), u.vecs, shift=2)
Draw(u)
Code:
---------------------------------------------------------------------------
NgException Traceback (most recent call last)
<ipython-input-16-af6cb1fe1cdb> in <module>
42 u = GridFunction(fes, multidim=20, name='resonances')
43 with TaskManager():
---> 44 lam = ArnoldiSolver(B2, B1, fes.FreeDofs(), u.vecs, shift=2)
45 Draw(u)
46
NgException: VHeight does not make sense for BlockMatrix