I'm trying to solve this Eigenvalue Problem and to plot some of the modes of the magnetic field (B) belonging to the eigenvalues with the largest real.
But the *.pvd appears to just containing the mesh?
Is it correct, that I have to use "FunctionSpace" if I'm using vector elements like N1curl? Because if I try using VectorSpace, I have to specify a tensor as Dirichlet BC (expecting value rank 2).
If I want to plot the B field, do I have to project or interpolate it onto a nodal VectorSpace first?
I'm still learning, so please excuse if the answers to my question are too obvious.
from dolfin import *
import scipy.special
import Speed as speed
e0 = Constant(1.0)
e1 = Constant(2.0)
e2 = Constant(3.0)
e3 = e2
Re = Constant(1E3)
mesh = Mesh("sphere.xml")
subdomains = MeshFunction("size_t", mesh, "sphere_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "sphere_facet_region.xml")
dx = Measure("dx")[subdomains]
sV = FunctionSpace(mesh, "N1curl", 1)
space = sV #MixedFunctionSpace(sV)
zeroVec = Constant((0.,0.,0.))
bc = DirichletBC(space, zeroVec, boundaries, 1)
trialB=TrialFunction(space)
dB = TestFunction(space)
t01 = speed.TFactor01(element = space.ufl_element())
s02 = speed.SFactor02(element = space.ufl_element())
s2c2 = speed.SFactor2c2(element = space.ufl_element())
s2s2 = speed.SFactor2s2(element = space.ufl_element())
v = (e0*curl(t01) + e1*curl(curl(s02)) + e2*curl(curl(s2c2)) +e3*curl(curl(s2s2)))
form_A_Inner = (-inner(curl(trialB), curl(dB)) + Re * inner(cross(v, trialB), curl(dB)))*dx(0)
form_A_Outer = -inner(curl(trialB), curl(dB))*dx(1)
form_B_Inner = inner(dB, trialB)*dx(0)
form_B_Outer = inner(dB, trialB)*dx(1)
form_A =form_A_Inner + form_A_Outer
form_B =form_B_Inner + form_B_Outer
A = PETScMatrix()
B = PETScMatrix()
assemble(form_A, tensor=A)
assemble(form_B, tensor=B)
bc.apply(A)
bc.apply(B)
eigensolver = SLEPcEigenSolver(A,B)
eigensolver.parameters["spectrum"] = "largest real"
eigensolver.parameters["problem_type"] = "gen_non_hermitian"
N = 5
eigensolver.solve(N)
for n in range(0, N):
r, c, rx, cx = eigensolver.get_eigenpair(n)
print 'Eigenvalue' + str(n) + ' ' + str(r) + ' + i * ' + str(c)
converged = eigensolver.get_number_converged()
print "Converged EV: %r" % converged
r, c, rx, cx = eigensolver.get_eigenpair(0)
B = Function(space)
B.vector()[:] = rx
fileB = File("Eigenvektor.pvd")
fileB << B
