1DMarkingProblem.py

from dolfin import *
import numpy as np
import math
from petsc4py import PETSc


#BCSType = 'Dirichlet'
#BCSType = 'Neumann'
BCSType = 'Mixed'
print BCSType
nx = 100
Z = Constant(1.)
dV = 0.0
mesh = IntervalMesh(nx, dV, 1.)

V = FunctionSpace(mesh, "CG", 2)


class boundaryLeft(SubDomain):
	def inside (self, x, on_boundary):
		return on_boundary and near(x[0],dV)
class boundaryRight(SubDomain):
	def inside (self, x, on_boundary):
		return on_boundary and near(x[0],1.)

boundaryLeft = boundaryLeft();
boundaryRight = boundaryRight(); 

def MarkBoundaries(mesh):
	sub_domains = MeshFunction("int", mesh, mesh.topology().dim()- 1)
	boundaries = FacetFunction("size_t", mesh)
	boundaryLeft.mark(sub_domains,1)
	boundaryRight.mark(sub_domains,2)
	
	ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)

def eigenvalues(V, bcs):

	p = TrialFunction(V)
	q = TestFunction(V)
	a = -inner(grad(p), grad(q))*dx
	if (BCSType == 'Mixed'):
		a += (-Z*q*p)*ds
	b = inner(q, p)*dx

	A = PETScMatrix()
	assemble(a, tensor = A)
	B = PETScMatrix()
	assemble(b, tensor = B)

	[bc.apply(A) for bc in bcs]
	[bc.zero(B) for bc in bcs]

	solver = SLEPcEigenSolver(A, B)
	solver.parameters["solver"] = "arpack"
	eps	 = solver.eps(); st	 = eps.getST();
	st.setType('sinvert'); ksp = st.getKSP(); ksp.setType('preonly');
	pc  = ksp.getPC(); pc.setType('lu'); pc.setFactorSolverPackage('mumps');
	
	solver.parameters["spectrum"] = "smallest magnitude"
	solver.parameters["spectral_transform"] = "shift-and-invert"
	solver.parameters["spectral_shift"] = 0.0

	neigs = 1
	solver.solve(neigs)

	computed_eigenvalues = []
	for i in range(min(neigs, solver.get_number_converged())):
		(r,im,rx,ix) = solver.get_eigenpair(i)
		#computed_eigenvalues.append(math.sqrt(abs(r))/math.pi)
		computed_eigenvalues.append([r,im])

	realParts = Function(V, rx)
	imagParts = Function(V, ix)

	return np.sort(np.array(computed_eigenvalues))

def print_eigenvalues(mesh):
	MarkBoundaries(mesh)

	if (BCSType == 'Dirichlet'):
		bcs = [DirichletBC(V, Constant(0.), boundaryLeft),
			   DirichletBC(V, Constant(0.), boundaryRight)]
	elif (BCSType == 'Neumann'):
		bcs = []
	elif (BCSType == 'Mixed'):
		bcs = [DirichletBC(V, Constant(0.), boundaryRight)]
	eig = eigenvalues(V, bcs); print eig

print_eigenvalues(mesh)