On 31 August 2010 19:27, Anders Logg <[email protected]> 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: [email protected] > To: Anders Logg <[email protected]> > 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 <[email protected]> > 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 : [email protected] > Unsubscribe : https://launchpad.net/~fenics > More help : https://help.launchpad.net/ListHelp > > _______________________________________________ Mailing list: https://launchpad.net/~fenics Post to : [email protected] Unsubscribe : https://launchpad.net/~fenics More help : https://help.launchpad.net/ListHelp
