← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to demo/pde/waveguide/python/demo.py

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/waveguide/python/demo.py
 
1
""" This demo demonstrates the calculation and visualization of a TM
 
2
(Transverse Magnetic) cutoff mode of a rectangular waveguide.
 
3
 
 
4
For more information regarding waveguides see 
 
5
 
 
6
http://www.ee.bilkent.edu.tr/~microwave/programs/magnetic/rect/info.htm
 
7
 
 
8
See the pdf in the parent folder and the following reference
 
9
 
 
10
The Finite Element in Electromagnetics (2nd Ed)
 
11
Jianming Jin [7.2.1 - 7.2.2]
 
12
 
 
13
"""
 
14
__author__ = "Evan Lezar evanlezar@gmail.com"
 
15
__date__ = "2008-08-22 -- 2008-12-17"
 
16
__copyright__ = "Copyright (C) 2008 Evan Lezar"
 
17
__license__  = "GNU LGPL Version 2.1"
 
18
 
 
19
# Modified by Anders Logg, 2008.
 
20
 
 
21
from dolfin import *
 
22
 
 
23
# Test for PETSc and SLEPc
 
24
try:
 
25
    PETScMatrix()
 
26
except:
 
27
    print "PyDOLFIN has not been configured with PETSc. Exiting."
 
28
    exit()
 
29
try:
 
30
    SLEPcEigenSolver()
 
31
except:
 
32
    print "PyDOLFIN has not been configured with SLEPc. Exiting."
 
33
    exit()
 
34
 
 
35
# Make sure we use the PETSc backend
 
36
dolfin_set("linear algebra backend", "PETSc")
 
37
 
 
38
# Create mesh
 
39
width = 1.0
 
40
height = 0.5
 
41
mesh = Rectangle(0, 0, width, height, 4, 2)
 
42
 
 
43
# Define the function space
 
44
V = FunctionSpace(mesh, "Nedelec", 2)
 
45
 
 
46
# Define the test and trial functions
 
47
v = TestFunction(V)
 
48
u = TrialFunction(V)
 
49
 
 
50
# Define the forms - generates an generalized eigenproblem of the form 
 
51
# [S]{h} = k_o^2[T]{h}
 
52
# with the eigenvalues k_o^2 representing the square of the cutoff wavenumber 
 
53
# and the corresponding right-eigenvector giving the coefficients of the 
 
54
# discrete system used to obtain the approximate field anywhere in the domain   
 
55
s = dot(curl_t(v), curl_t(u))*dx
 
56
t = dot(v, u)*dx
 
57
 
 
58
# Assemble the stiffness matrix (S) and mass matrix (T)
 
59
S = PETScMatrix()
 
60
T = PETScMatrix()
 
61
assemble(s, tensor=S)
 
62
assemble(t, tensor=T)
 
63
 
 
64
# Solve the eigensystem
 
65
esolver = SLEPcEigenSolver()
 
66
esolver.set("eigenvalue spectrum", "smallest real")
 
67
esolver.set("eigenvalue solver", "lapack")
 
68
esolver.solve(S, T)
 
69
 
 
70
# The result should have real eigenvalues but due to rounding errors, some of 
 
71
# the resultant eigenvalues may be small complex values. 
 
72
# only consider the real part
 
73
 
 
74
# Now, the system contains a number of zero eigenvalues (near zero due to 
 
75
# rounding) which are eigenvalues corresponding to the null-space of the curl 
 
76
# operator and are a mathematical construct and do not represent physically 
 
77
# realizable modes.  These are called spurious modes.  
 
78
# So, we need to identify the smallest, non-zero eigenvalue of the system - 
 
79
# which corresponds with cutoff wavenumber of the the dominant cutoff mode.
 
80
cutoff = None
 
81
for i in range(S.size(1)):
 
82
    (lr, lc) = esolver.getEigenvalue(i)
 
83
    if lr > 1 and lc == 0:
 
84
        cutoff = sqrt(lr)
 
85
        break
 
86
 
 
87
if cutoff is None:
 
88
    print "Unable to find dominant mode"
 
89
else:
 
90
    print "Cutoff frequency:", cutoff