Contour plot of SLEPcEigenSolver eigenvector solution

Question

In order to do a mode analysis on step-index fibers and later PCFs I wanted to check rectangular waveguides to see that everything is working. For a rectangular waveguide I found several code snippets, but I cannot get the final contour plotting of the fundamental mode to work. I assign the eigenvector to my combined function space but I am wondering how I can access the vector components maybe even outside of dolfin with some standard python routines.

A running example is:

from dolfin import *

parameters["linear_algebra_backend"] = "PETSc"

mesh used for the rectangular hollow guides

a = 1.0

b = 0.5

create a rectangular mesh with origin (0,0) extending to (a,b) with 8 edges along the long side and 4 elements along the short side

mesh = RectangleMesh(0, 0, a, b, 8, 4)

e_r = 1.0

u_r = 1.0

vector_order = 2

nodal_order = 2

define the functions spaces

V = FunctionSpace (mesh, "Nedelec 1st kind H(curl)", vector_order)

Q = FunctionSpace(mesh, "Lagrange", nodal_order)

W = V * Q 

define the test and trial functions from the combined space here N_v and N_u# are Nedelec basis functions and L_v and L_u are Lagrange basis functions

(N_i, L_i) = TestFunctions(W)

(N_j, L_j) = TrialFunctions(W)

define the forms (matrix elements) for cutoff analysis for the basis functio# s

def curl_t(w):

    return Dx(w[1],0) - Dx(w[0],1)

s_tt_ij = 1.0/u_r * dot (curl_t(N_i), curl_t(N_j))

t_tt_ij = e_r * dot (N_i, N_j)

s_zz_ij = 1.0/u_r* dot (grad(L_i), grad(L_j))

t_zz_ij = e_r * L_i * L_j

post-multiplication by dx will result in integration over the domain of the mesh at assembly time

s_ij = (s_tt_ij + s_zz_ij) * dx

t_ij = (t_tt_ij + t_zz_ij) * dx

assemble the system matrices. DOLFIN automatically evaluates each of the forms for all the relevant test and trial function combinations ie. all possible v

alues of i and j

S = PETScMatrix()

T = PETScMatrix()

assemble (s_ij, tensor=S)

assemble (t_ij, tensor=T)

define the subdomains with BC

class PECWall(SubDomain):

        def inside(self, x, on_boundary):

            return on_boundary;

electric_wall = DirichletBC(W, Expression (("0.0", "0.0", "0.0") ) , PECWall())

apply the BC to assembled matrices for the cutoff problem:

electric_wall.apply(S)

electric_wall.apply(T)

esolver = SLEPcEigenSolver(S,T)

esolver.parameters["spectrum"] = "smallest real"

esolver.parameters["solver"] = "lapack"

esolver.solve()

cutoff = None

print S.size(1)

for i in range(S.size(1)):

        (lr,lc,vr,vc) = esolver.get_eigenpair(i)

        if lr > 1 and lc == 0:

            cutoff = sqrt(lr)

            break

if cutoff is None:

        print "Unable to find dominant mode"

else:

        print "Cutoff frequency: ", cutoff

f = Function(W, vr)

(transverse,axial) = f.split()

print axial

plot(transverse)

plot(axial)

Answer 1

Take a look at the demo in demo/documented/eigenvalue. Write the eigenvector 'function' to a PVD or XDMF file and use ParaView to open it.

Comments

I wrote the eigenvector 'function' u (as in the demo_eigenvalue.py) to a pvd file, but that file an the associated vtu are empty.

Below

u = Function(V)
u.vector()[:] = rx

I added

file = File('eigenvalue.pvd')
file << u

Any idea what I am doing wrong here?

---

Post the shortest possible (but complete) code that leads to an 'empty' pvd file.

---

based on the demo_eigenvalue.py

from dolfin import *
# Test for PETSc and SLEPc
if not has_linear_algebra_backend("PETSc"):
    print "DOLFIN has not been configured with PETSc. Exiting."
    exit()

if not has_slepc():
    print "DOLFIN has not been configured with SLEPc. Exiting."
    exit()

# Define mesh, function space
mesh = Mesh("box_with_dent.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 1)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A)

# Compute all eigenvalues of A x = \lambda x
print "Computing eigenvalues. This can take a minute."
eigensolver.solve()

# Extract largest (first) eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(0)

print "Largest eigenvalue: ", r

# Initialize function and assign eigenvector
u = Function(V)
u.vector()[:] = rx

file = File('eigenvalue.pvd')
file << u