fenics team mailing list archive

[Kristian Ølgaard <k.b.oelgaard@xxxxxxxxx>]

On 31 August 2010 19:27, Anders Logg <logg@xxxxxxxxx> wrote:
> 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.

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

Kristian

> --
> Anders
>
>
> ---------- Forwarded message ----------
> From: noreply@xxxxxxxxxxxxx
> To: Anders Logg <logg@xxxxxxxxx>
> Date: Tue, 31 Aug 2010 17:07:24 -0000
> Subject: [Branch ~fenics-core/fenics-doc/main] Rev 138: Add python documentation for simple eigenvalue demo. Need to do stuff
> ------------------------------------------------------------
> 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
>  -------------
>
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to     : fenics@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~fenics
> More help   : https://help.launchpad.net/ListHelp
>
>