Question #141904 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/141904
Status: Answered => Open
Raphael Kruse is still having a problem:
Hi,
I wrote the following code to solve the problem. I would be very
grateful if someone with experience could go over it. Very likely it is
possible to make the code more efficient.
"""
This scripts builds a FEM-approximation of a Carleman operator Q, which is
given by an integral kernel q(x,y). Then it solves the corresponding
eigenproblem of finding all eigenpairs (lambda,u) such that Qu = lambda*u.
Q is given by [Qu](y) = \int q(x,y) u(x) dx
Author: Raphael Kruse (Univerity Bielefeld), 26 Jan 2011
"""
from dolfin import *
import numpy as np
from numpy import linalg as LA
# parameters for the kernel q
sigma = 1 # standard deviation
gamma = 0.2 # correlation length
dof = 200 # degrees of freedom
tol = 1E-14 # tolerance for eigenvalues
# Defining the integral kernel
q_1d = Expression('sigma*sigma*exp(-pow(x[0]-y,2)/(gamma*gamma))')
q_1d.sigma = sigma
q_1d.gamma = gamma
# Test for PETSc and SLEPc
if not has_la_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 = UnitInterval(dof-1)
coor = mesh.coordinates()
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)
# building the mass matrix M as PETScMatrix
print "Building the mass matrix M...",
M = PETScMatrix()
a = u*v*dx
assemble(a, tensor=M)
print "Done"
# Building the Qu matrix for u running through the testfunctions
print "Building the Q matrix...",
Qu = np.zeros( (dof,dof) ) # holds all values of convolutions of q with
# testfunctions
for i in range(dof):
q_1d.y = coor[i][0] # y runs through the coordinates of the mesh
L = q_1d*v*dx
b = assemble(L, mesh=mesh)
Qu[i,:] = b.array()
Qu = Qu.transpose()
# Computing the values of \int Qu(y) v(y) dy for each u,v
Q = np.zeros( (dof, dof) ) # the matrix Q
Fct_Qu = Function(V)
for k in range(dof):
quk = np.array(Qu[k][:]) # runs through the values of u
Fct_Qu.vector()[:] = quk
L = Fct_Qu*v*dx
b = assemble(L, mesh=mesh)
Q[:,k] = b.array()
print "Done"
# transforming Q into PETScMatrix
print "Transforming Q into PETScMatrix format...",
Qh = PETScMatrix(M)
for i in range(dof):
for k in range(dof):
Qh[i,k] = Q[i,k]
print "Done"
#print Qh.array()
# Create eigensolver
eigensolver = SLEPcEigenSolver()
eigensolver.parameters["solver"] = "lapack"
# Compute all eigenvalues of Q x = \lambda M x
--
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@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
