This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Setting Boundary Condition on Eigenvalue Problem.

0 votes

Hi,

Below is code to plot the first three eigenvalues in an infinite potential well
(adapted from http://en.wikiversity.org/wiki/User:Tclamb/FEniCS#Quantum_Harmonic_Oscillator )

While the boundary condition (solution = 0 at the ends) is defined, it doesn't seem to be implemented. Indeed, when the eigenfunctions are plotted by the code, not all of them obey the boundary condition.

from dolfin import *

#define mesh and function space
mesh = IntervalMesh(100, -10, 10)
V = FunctionSpace(mesh, 'Lagrange', 3)

#define boundary
def boundary(x, on_boundary):
    return on_boundary

#apply essential boundary conditions
bc = DirichletBC(V, 0, boundary)

#define functions
u = TrialFunction(V)
v = TestFunction(V)

Pot = Expression('0.0')

#define problem
a = (inner(grad(u), grad(v)) \
     + Pot*u*v)*dx
m = u*v*dx

#assemble stiffness matrix
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)

#create eigensolver
eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters['solver']   = 'lapack'
eigensolver.parameters['tolerance'] = 1.e-15

#solve for eigenvalues
eigensolver.solve()

u = Function(V)
for i in range(0,3):
    #extract next eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print 'eigenvalue:', r

    #assign eigenvector to function
    u.vector()[:] = rx

    #plot eigenfunction
    plot(u)
    interactive()
asked Apr 11, 2014 by sixtysymbols FEniCS User (2,280 points)

1 Answer

+3 votes
 
Best answer

Hi, you have created object for Dirichlet bcs but haven't really used it.

#assemble stiffness matrix
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)
bc.apply(A)          # apply the boundary conditions
bc.apply(M)
answered Apr 11, 2014 by MiroK FEniCS Expert (80,920 points)
selected Apr 11, 2014 by sixtysymbols

Thanks for this.

P.S. It looks like the solver would fail for higher eigenvalues.

If I only apply the boundary conditions to A, however, it works fine.

...