--- Begin Message ---
------------------------------------------------------------
revno: 138
committer: Marie E. Rognes <meg@xxxxxxxxx>
branch nick: fenics-doc
timestamp: Tue 2010-08-31 19:03:34 +0200
message:
Add python documentation for simple eigenvalue demo. Need to do stuff
with mesh. Where should we put meshes that are not built-in?
modified:
source/demos/la/eigenvalue/simple/common.txt
source/demos/la/eigenvalue/simple/python/demo.py
source/demos/la/eigenvalue/simple/python/documentation.rst
--
lp:fenics-doc
https://code.launchpad.net/~fenics-core/fenics-doc/main
You are subscribed to branch lp:fenics-doc.
To unsubscribe from this branch go to https://code.launchpad.net/~fenics-core/fenics-doc/main/+edit-subscription
=== modified file 'source/demos/la/eigenvalue/simple/common.txt'
--- source/demos/la/eigenvalue/simple/common.txt 2010-08-31 15:26:56 +0000
+++ source/demos/la/eigenvalue/simple/common.txt 2010-08-31 17:03:34 +0000
@@ -1,7 +1,33 @@
This demo illustrates how to:
-*
-
-Equation and problem definition
--------------------------------
+* Load a mesh from a file
+* Solve an eigenvalue problem
+* Use a specific linear algebra backend (PETSc)
+* Initialize a finite element function with a coefficient vector
+
+Problem definition
+------------------
+
+Sometimes one wants to solve an eigenvalue problem such as this one:
+find the eigenvalues :math:`\lambda \in \mathbb{R}` and the
+corresponding eigenvectors :math:`x \in \mathbb{R^n}` such that
+
+.. math::
+
+ A x = \lambda x
+
+In the finite element world, the matrix :math:`A` often originates
+from some partial differential operator. For instance, :math:`A` can
+be the stiffness matrix corresponding to this bilinear form:
+
+.. math::
+
+ a(u, v) = \int_{\Omega} \nabla u \cdot \nabla v \ {\rm d} x.
+
+Here, we will let the space :math:`V` (of dimension :math:`n`) consist
+of continuous piecewise linear polynomials defined relative to some
+mesh (Lagrange finite elements). For this example, we will consider a
+3-D mesh of tetrahedra generated elsewhere.
+
+In the following, we show how this eigenvalue problem can be solved.
=== modified file 'source/demos/la/eigenvalue/simple/python/demo.py'
--- source/demos/la/eigenvalue/simple/python/demo.py 2010-08-31 15:26:56 +0000
+++ source/demos/la/eigenvalue/simple/python/demo.py 2010-08-31 17:03:34 +0000
@@ -23,32 +23,33 @@
print "DOLFIN has not been configured with SLEPc. Exiting."
exit()
-# Make sure we use the PETSc backend
-parameters["linear_algebra_backend"] = "PETSc"
-
-# Define mesh, function space and basis functions
-mesh = UnitSquare(4, 4)
+# Define mesh, function space
+mesh = Mesh("mesh.xml.gz")
V = FunctionSpace(mesh, "CG", 1)
+
+# Define basis and then form the stiffness matrix
+u = TrialFunction(V)
v = TestFunction(V)
-u = TrialFunction(V)
-
-# Define form for stiffness matrix
a = dot(grad(v), grad(u))*dx
# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)
+# Create eigensolver
+eigensolver = SLEPcEigenSolver()
+
# Compute all eigenvalues of A x = \lambda x
-esolver = SLEPcEigenSolver()
-esolver.solve(A)
-
-# Get largest (first) eigenpair
-a, b, x, y = esolver.get_eigenpair(0)
-print "Largest eigenvalue: ", a
+print "Computing eigenvalues. This can take a minute."
+eigensolver.solve(A)
+
+# Extract largest (first) eigenpair
+r, c, rx, cx = eigensolver.get_eigenpair(0)
+
+print "Largest eigenvalue: ", r
# Initialize function with eigenvector
-u = Function(V, x)
+u = Function(V, rx)
# Plot eigenfunction
plot(u)
=== modified file 'source/demos/la/eigenvalue/simple/python/documentation.rst'
--- source/demos/la/eigenvalue/simple/python/documentation.rst 2010-08-31 15:26:56 +0000
+++ source/demos/la/eigenvalue/simple/python/documentation.rst 2010-08-31 17:03:34 +0000
@@ -3,14 +3,115 @@
A simple eigenvalue solver
==========================
+We recommend that you are familiar with demo for the Poisson equation
+before looking at this demo.
+
+
.. include:: ../common.txt
+
+If you want a more complex problem, we suggest that you look at the
+other eigenvalue demo.
+
+
+
Implementation
--------------
This demo is implemented in a single Python file, :download:`demo.py`,
which contains both the variational forms and the solver.
+The eigensolver functionality in DOLFIN relies on the library SLEPc
+which in turn relies on the linear algebra library PETSc. Therefore,
+both PETSc and SLEPc are required for this demo. We can test whether
+PETSc and SLEPc are available, and exit if not, as follows:
+
+.. code-block:: python
+
+ from dolfin import *
+
+ # 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()
+
+First, we need to construct the matrix :math:`A`. This will be done in
+three steps: define the mesh and the function space associated with
+it; constructing the variational form defining the matrix; and then
+assembling this form. The code is shown below
+
+.. code-block:: python
+
+ # Define mesh, function space
+ mesh = Mesh("mesh.xml.gz")
+ V = FunctionSpace(mesh, "CG", 1)
+
+ # Define basis and then form the stiffness matrix
+ u = TrialFunction(V)
+ v = TestFunction(V)
+ a = dot(grad(v), grad(u))*dx
+
+ # Assemble stiffness form
+ A = PETScMatrix()
+ assemble(a, tensor=A)
+
+Note that we (in this example) first define the matrix ``A`` as a
+``PETScMatrix`` and then assemble the form into it. This is an easy
+way to ensure that the matrix has the right type.
+
+In order to solve the eigenproblem, we need an eigensolver:
+
+.. code-block:: python
+
+ # Create eigensolver
+ eigensolver = SLEPcEigenSolver()
+
+Now, we ready solve the eigenproblem by calling the ``solve`` method
+of the eigensolver, giving the matrix ``A`` as a single argument. Note
+that eigenvalue problems tend to be computationally intensive and may
+hence take a while.
+
+.. code-block:: python
+
+ # Compute all eigenvalues of A x = \lambda x
+ print "Computing eigenvalues. This can take a minute."
+ eigensolver.solve(A)
+
+The result is kept by the eigensolver, but can fortunately be
+extracted. Here, we have computed all eigenvalues, and they will be
+sorted by largest magnitude. We can extract the real and complex part
+(``r`` and ``c``) of the largest eigenvalue and the real and complex
+part of the corresponding eigenvector (``ru`` and ``cu``) by asking
+for the first eigenpair as follows:
+
+.. code-block:: python
+
+ # Extract largest (first) eigenpair
+ r, c, rx, cx = eigensolver.get_eigenpair(0)
+
+Finally, we want to examine the results. The eigenvalue can easily be
+printed. But, the real part of eigenvector is probably most easily
+visualized by constructing the corresponding eigenfunction. This can
+be done by creating a ``Function`` in the function space ``V`` with
+the eigenvector ``rx`` as an additional argument. Then the
+eigenfunction can be manipulated as any other ``Function``, and in
+particular plotted:
+
+.. code-block:: python
+
+ print "Largest eigenvalue: ", r
+
+ # Initialize function with eigenvector
+ u = Function(V, rx)
+
+ # Plot eigenfunction
+ plot(u)
+ interactive()
+
Complete code
-------------
--- End Message ---
