dolfin team mailing list archive

[Evan Lezar <mail@xxxxxxxxxxxxx>]

Pablo,

I have been using FEniCS for solving waveguide cutoff and dispersion
problems, and can send you a few code samples if you like. In my problems I
use homogeneous Dirichlet BCs, but this should not be too much of a problem
to extend to the non-homogeneous case.

As a start, I have attached a simple example for calculating the TE cutoff
modes of a hollow rectangular waveguide. To extend this to the
non-homogeneous case, you would have to use a non-zero expression when
instantiating the DirichletBC class.

Feel free to contact me if you have any questions.

Regards
--
Evan Lezar

Computational Electromagnetics Group
Department of Electrical and Electronic Engineering
University of Stellenbosch
Stellenbosch
South Africa

www.evanlezar.com

GoogleTalk: evanlezar
Skype: evanlezar


On Sun, May 8, 2011 at 12:11 AM, Pablo Suarez <
question156430@xxxxxxxxxxxxxxxxxxxxx> wrote:

> Question #156430 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/156430
>
> Description changed to:
> Hello,
>
>      I was wondering if anyone can answer this question.
> I need to solve for the follow PDE:
>
>        \nabla \times \nabla E - k^2 E = 0
>       (curl(curl(E)) - k^2 E = 0
>
> with BC   n \times E = F
> E is a vector field
> F is a vector field
> k^2 is a constant (k^2 is not an eigenvalue)
> this is in the interior of a 3D surface (presumably smooth)
> n is the outward normal  to the surface.
>
> This is the cavity problem for Maxwell Equation's  with non-homgenous
> Dirichlet BC.
> How can I incorporate this condition into my code?
> I have followed the example in undocumented demos curl curl. If anyone can
> help it will be greatly appreciated.
>
> Best wishes,
>
> Pablo
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>

__author__ = "Evan Lezar mail@xxxxxxxxxxxxx"
__date__ = "12 July 2010"
__copyright__ = "Copyright (C) 2010 Evan Lezar"

"""
A simple example for the use of pydolfin for the solution of the TE modes of a hollow rectangular waveguide.
"""
from time import time
import numpy as np
from dolfin import *

t0 = time()
# PEC walls for the application of dirichlet boundary conditions
class PECWalls(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary; 

# since the problem is constrained to the plane of the waveguide port, a transverse curl is required
def curl_t(w):
    return Dx(w[1], 0) - Dx(w[0], 1)

# Setup the mesh
width = 1.0
height = 0.5
mesh = Rectangle(0, 0, width, height, 8, 4, "right")


order = 2
# Use Nedelec basis functions of the first (curl conforming)
V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", order)
# Define the test and trial functions
v = TestFunction(V)
u = TrialFunction(V)

# Define the matrix elements
s = curl_t(v)*curl_t(u)*dx
t = inner(v, u)*dx

# Assemble the stiffness matrix (S) and mass matrix (T)
S = PETScMatrix()
T = PETScMatrix()
assemble(s, tensor=S)
assemble(t, tensor=T)

# setup the dirichlet boundary condition and apply this to the matrices
boundary_condition = DirichletBC(V, Expression(("0.0", "0.0"), element=V.ufl_element()), PECWalls())
b = PETScVector(S.size(0))
boundary_condition.apply(S, b)
boundary_condition.apply(T, b)


t_solver = time()
# Create the eigensolver and solve the eigensystem
esolver = SLEPcEigenSolver()
esolver.parameters["spectrum"] = "largest magnitude"
#esolver.parameters["spectral_transform"] = "shift-and-invert"
#esolver.parameters["spectral_shift"] = 10.0
esolver.parameters["solver"] = 'lapack'

# Process the computed eigenvalues and eigenvectors
num_modes_to_compute = 20 # The number of modes to be computed

# Solve the inverse problem: T*x = 1/k_o^2*S*x as spurious modes are shifted to infinity
esolver.solve(T, S)

t_solver = time() - t_solver



# store the results
mode_results = {'k_o_squared':[], 'coefficients':[]}


mode_count = 0
skip_count = 0;
for i in range(S.size(1)):
    if mode_count >= num_modes_to_compute:
        break
    (lr, lc) = esolver.get_eigenvalue(i)
    print lr + 1j*lc
#            check to see if it is a valid eigenpair
    if np.isfinite(lr) and np.isfinite(lc):
        if (lr == 1.0) and (lc == 0.0):
#                 Mode inserted by applying boudary conditions
            pass
        else:
#                    Is a valid mode
            (lr, lc, dofs_r_petsc, dofs_c_petsc) = esolver.get_eigenpair(i)
            
            k_o_squared = 1/(lr + 1j*lc)
            print k_o_squared
#            store the cutoff wavenumber squared as well as real and imaginary coefficients of the basis functions
            mode_results['k_o_squared'].append(np.complex(k_o_squared.real, k_o_squared.imag))
            mode_results['coefficients'].append((dofs_r_petsc, dofs_c_petsc))
            
            mode_count += 1
    else:
#                spurious mode
        pass

# theoretical cutoff squared values (only the first 4)
# k_o_squared/pi**2 = (m/a)**2 + (n/b)**2
# with a = 1.0m and b = 0.5m
# the first four modes are TE_10, TE_01, TE_20, and TE_11
theoretical_values = np.array([1,4,4,5])*np.pi**2
# Output the results
for current_mode in range(num_modes_to_compute):
    print "Square of the cutoff wavenumber for mode number %d: " % (current_mode+1),
    print mode_results['k_o_squared'][current_mode]
    if current_mode < 4:
        print "(theoretical value = %f)" % theoretical_values[current_mode]
        


t_total = time() - t0

print 'Total time: %f. Eigensolver time: %f' % (t_total, t_solver)

if False:
    # Plot a mode - use the real part of the computed coefficients ( the last 0 index )
    mode_to_plot = 3
    if mode_to_plot < num_modes_to_compute:
        # use only the real coefficients to instantiate the function
        mode_function = Function(V, mode_results['coefficients'][mode_to_plot][0])
        # this mode function can now be plotted
        plot(mode_function)
        interactive()