Subject: [Fenics] [[email-masked]: [Branch ~fenics-core/fenics-doc/main] Rev 138: Add python documentation for simple eigenvalue demo. Need to do stuff]

From: Anders Logg <[email-masked]>
Date: Tue, 31 Aug 2010 19:27:53 +0200

I just put the meshes in the same directory, one copy in each of the
C++ and Python directories. I don't think we need to worry about the
size of the meshes.

--
Anders

Subject: [Branch ~fenics-core/fenics-doc/main] Rev 138: Add python documentation for simple eigenvalue demo. Need to do stuff
From: [email-masked]
To: Anders Logg <[email-masked]>
Date: Tue, 31 Aug 2010 17:07:24 -0000

------------------------------------------------------------
revno: 138
committer: Marie E. Rognes <[email-masked]>
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
 -------------

---

From: Marie Rognes <[email-masked]>
Date: Tue, 31 Aug 2010 19:29:06 +0200

> C++ and Python directories. I don't think we need to worry about the

Ok!

--
Marie

---

From: Kristian Ølgaard <[email-masked]>
Date: Tue, 31 Aug 2010 21:01:23 +0200

> C++ and Python directories. I don't think we need to worry about the

A bit late, but I also vote for adding meshes locally, to make the
demos self contained.

Kristian

> --
> Anders
>
>
> ---------- Forwarded message ----------
> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~fenics
> More help   : https://help.launchpad.net/ListHelp
>
>