diff --git a/INSTALLATION.txt b/INSTALLATION.txt
index 6d24c0ff87..4f42ac5f5c 100644
--- a/INSTALLATION.txt
+++ b/INSTALLATION.txt
@@ -168,10 +168,8 @@ If necessary, you can download_ and install it for your platform
 .. note::
 
    :term:`FiPy` requires at least version 2.4.x of :term:`Python`. See
-   :ref:`RunningUnderPython3` for instructions on how to run
-   :term:`FiPy` with `Python 3.x`_.
+   the specialized instructions if you wish to :ref:`RunUnderPython3`.
 
-.. _Python 3.x:   http://www.nist.gov/cgi-bin/exit_nist.cgi?url=http://docs.python.org/py3k/
 .. _download: http://www.nist.gov/cgi-bin/exit_nist.cgi?url=http://www.python.org/download/
 
 :term:`Python` along with many of :term:`FiPy`'s required and optional
diff --git a/README.txt b/README.txt
index 2f5c6fdf26..ed1fa20656 100644
--- a/README.txt
+++ b/README.txt
@@ -54,25 +54,91 @@ necessary.
 
 The significant changes since version 2.1 are:
 
+- :ref:`CoupledEquations` are now supported.
+- A more robust mechanism for specifying :ref:`BoundaryConditions` is now 
+  used.
+- Most :class:`~fipy.meshes.mesh.Mesh`\es can be partitioned by 
+  :ref:`MeshingWithGmsh`.
+- :ref:`PYAMG` and :ref:`SCIPY` have been added to the :ref:`SOLVERS`.
+- FiPy is capable of :ref:`RunningUnderPython3`.
+- "getter" and "setter" methods have been pervasively changed to Python 
+  properties.
+- The test suite now runs much faster.
 - Tests can now be run on a full install using `fipy.test()`.
-- Grid classes now take an `Lx` argument. 
 - The functions of the :mod:`~fipy.tools.numerix` module are no longer 
   included in the :mod:`fipy` namespace. See :mod:`examples.updating.update2_0to3_0` 
   for details.
-- Support for Python 3. Please see :ref:`RunningUnderPython3` for details.
+- Equations containing a :class:`~fipy.terms.transientTerm.TransientTerm`,
+  must specify the timestep by passing a ``dt=`` argument when calling
+  :meth:`~fipy.terms.term.Term.solve` or :meth:`~fipy.terms.term.Term.sweep`.
 
 Tickets fixed in this release::
 
-    171	update the mayavi viewer to use  mayavi 2
-    286	'matplotlib: list index out of range' when no title given, but only sometimes
-    197	~binOp doesn't work on branches/version-2_0
-    194	`easy_install` instructions for MacOSX are broken
-    192	broken setuptools url with python 2.6
-    184	The FiPy webpage seems to be broken on Internet Explorer
-    168	Switch documentation to use `:math:` directive
-    198	FiPy2.0.2 LinearJORSolver.__init__  calls Solver rather than PysparseSolver
-    199	`gmshExport.exportAsMesh()` doesn't work
-    195	broken arithmetic face to cell distance calculations
+    45  Navier Stokes
+    85  CellVariable hasOld() should set self.old
+    101 Grids should take Lx, Ly, Lz arguments
+    145 tests should be run with fipy.tests()
+    177 remove ones and zeros from numerix.py
+    178 Default time steps should be infinite
+    291 term multiplication changes result
+    296 FAQ gives bad guidance for anisotropic diffusion
+    297 Use physical velocity in the manual/FAQ
+    298 mesh manipulation of periodic meshes leads to errors
+    299 Give helpfull error on - or / of meshes
+    301 wrong cell to cell normal in periodic meshes
+    302 gnuplot1d gives error on plot of facevariable
+    309 pypi is failing
+    312 Fresh FiPy gives ""ImportError: No viewers found"""
+    314 Absence of enthought.tvtk causes test failures
+    319 mesh in FiPy name space
+    324 --pysparse configuration should never attempt MPI imports
+    327 factoryMeshes.py not up to date with respect to keyword arguments
+    331 changed constraints don't propagate
+    332 anisotropic diffusion and constraints don't mix
+    333 `--Trilinos --no-pysparse` uses PySparse?!?
+    336 Profile and merge reconstrain branch
+    339 close out reconstrain branch
+    341 Fix fipy.terms._BinaryTerm test failure in parallel
+    343 diffusionTerm(var=var1).solver(var=var0) should fail sensibly
+    346 TeX is wrong in examples.phase.quaternary
+    348 Include Benny's improved interpolation patch
+    354 GmshExport is not tested and does not work
+    355 Introduce mesh.x as shorthand for mesh.cellCenters[0] etc
+    356 GmshImport should support all element types
+    357 GmshImport should read element colors
+    363 Reduce the run times for chemotaxis tests
+    366 tests take *too* long!!!
+    369 Make DiffusionTermNoCorrection the default
+    370 Epetra Norm2 failure in parallel
+    373 remove deprecated `steps=` from Solver
+    376 remove deprecated `diffusionTerm=` argument to ConvectionTerm
+    377 remove deprecated `NthOrderDiffusionTerm`
+    380 remove deprecated Variable.transpose()
+    381 remove deprecated viewers.make()
+    382 get running in Py3k
+    384 gmsh importer and gmsh tests don't clean up after themselves
+    385 `diffusionTerm._test()` requires PySparse
+    390 Improve test reporting to avoid inconsequential buildbot failures
+    391 efficiency_test chokes on liquidVapor2D.py
+    393 two `--scipy` failures
+    395 `--pysparse --inline` failures
+    417 Memory consumption growth with repeated meshing, especially with Gmsh
+    418 Viewers not working when plotting meshes with zero cells in parallel
+    419 examples/cahnHilliard/mesh2D.py broken with --trilinos
+    420 Epetra.PyComm() broken on Debian
+    421 cellVariable.min() broken in parallel
+    426 Add in parallel buildbot testing on more than 2 processors
+    427 Slow PyAMG solutions
+    434 Gmsh I/O
+    438 changes to gmshImport.py caused --inline problems
+    439 gmshImport tests fail on Windows due to shared file
+    441 Explicit convetion terms should fail when the equation has no TransientTerm (dt=None)
+    445 getFaceCenters() should return a FaceVariable
+    446 constraining values with ImplictSourceTerm not documented?
+    448 Gmsh2D does not respect background mesh
+    452 Gmsh background mesh doesn't work in parallel
+    453 faceValue as FaceCenters gives inline failures
+    454 Py3k and Windows test failures
 
 .. warning::
 
diff --git a/documentation/FAQ.txt b/documentation/FAQ.txt
index b1439ca739..7979681e62 100644
--- a/documentation/FAQ.txt
+++ b/documentation/FAQ.txt
@@ -97,7 +97,7 @@ How do I represent a `...` term that *doesn't* involve the dependent variable?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 It is important to realize that, even though an expression may
-superficially resemble one of those shown above, if the dependent variable
+superficially resemble one of those shown in :ref:`section:discretization`, if the dependent variable
 *for that PDE* does not appear in the appropriate place, then that
 term should be treated as a source.
 
@@ -193,6 +193,12 @@ For :term:`FiPy`'s purposes, however, this term represents the convection of
 >>> eq = TransientTerm() == (DiffusionTerm(coeff=D1) 
 ...                          + <Specific>ConvectionTerm(coeff=D2 * xi.faceGrad))
 
+.. note::
+
+   With the advent of :ref:`CoupledEquations` in FiPy 3.x, it is now
+   possible to represent both terms with
+   :class:`~fipy.terms.diffusionTerm.DiffusionTerm`.
+
 What if the coefficient of a term depends on the variable that I'm solving for?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
@@ -296,9 +302,11 @@ want, and then issue supplemental commands for the underlying plotting
 package. The better option is to make a "subclass" of the :term:`FiPy`
 :class:`Viewer <viewer.AbstractViewer>` that comes closest to producing the image you want. You can
 then override just the behavior you wan to change, while letting :term:`FiPy` do
-most of the heavy lifting. See :mod:`examples.phase.anisotropy` for an
-example of creating a custom :term:`Matplotlib` :class:`Viewer 
-<viewer.AbstractViewer>` class.
+most of the heavy lifting. See :mod:`examples.phase.anisotropy` and
+:mod:`examples.phase.polyxtal` for examples of creating a custom
+:term:`Matplotlib` :class:`Viewer <viewer.AbstractViewer>` class; see
+:mod:`examples.cahnHilliard.sphere` for an example of creating a custom
+:term:`Mayavi` :class:`Viewer <viewer.AbstractViewer>` class.
 
 .. _FAQ-IterationsTimestepsSweeps:
 
@@ -390,7 +398,7 @@ sweeps
   Sweeps are used to achieve better solutions in
   :mod:`examples.diffusion.mesh1D`, 
   :mod:`examples.phase.simple`, 
-  :mod:`examples.phase.binary`, and :mod:`examples.flow.stokesCavity`.
+  :mod:`examples.phase.binaryCoupled`, and :mod:`examples.flow.stokesCavity`.
 
 timesteps
   This outermost layer of repetition is of most practical interest to
@@ -436,7 +444,7 @@ timesteps
   chosen based on the expected interfacial velocity in
   :mod:`examples.phase.simple`. The timestep is gradually
   increased as the kinetics slow down in
-  :mod:`examples.cahnHilliard.mesh2D`.
+  :mod:`examples.cahnHilliard.mesh2DCoupled`.
 
 Finally, we can (and often do) combine all three layers of repetition:
 
@@ -496,111 +504,7 @@ How do I represent boundary conditions?
 
 .. currentmodule:: fipy.variables.cellVariable
 
-How do I represent a fixed value (Dirichlet) boundary condition?
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Use the :meth:`~CellVariable.constrain` method. For
-example, to fix `var` to have a value of `2` along the upper surface of a domain,
-use
-
->>> var.constrain(2., where=mesh.facesTop)
-
-.. note:: 
-
-   The old equivalent
-   :class:`~fipy.boundaryConditions.fixedValue.FixedValue` boundary
-   condition is now deprecated.
-
-How do I apply a (Neumann) fixed gradient boundary condition?
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Use the :attr:`~.CellVariable.faceGrad`.\ :meth:`~fipy.variables.variable.Variable.constrain`
-method. For example, to fix `var` to have a gradient of `(0,2)` along the upper
-surface of a 2D domain, use
-
->>> var.faceGrad.constrain(((0,),(2,)), where=mesh.facesTop)
-
-How do I apply a fixed flux boundary condition?
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Generally this can be implemented with a judicious use of
-:attr:`~.CellVariable.faceGrad`.\ :meth:`~fipy.variables.variable.Variable.constrain`.
-Failing that, an exterior flux term can be added to the equation. Firstly,
-set the terms' coefficients to be zero on the exterior faces,
-
->>> diffCoeff.constrain(0., mesh.exteriorFaces)
->>> convCoeff.constrain(0., mesh.exteriorFaces)
-
-then create an equation with an extra term to account for the exterior flux,
-
->>> eqn = (TransientTerm() + ConvectionTerm(convCoeff) 
-...        == DiffusionCoeff(diffCoeff)
-...        + (mesh.exteriorFaces * exteriorFlux).divergence)
-
-where `exteriorFlux` is a rank 1
-:class:`~fipy.variables.faceVariable.FaceVariable`.
-
-.. note:: 
-
-   The old equivalent :class:`~fipy.boundaryConditions.fixedFlux.FixedFlux`
-   boundary condition is now deprecated.
-
-How do I apply an outlet or inlet boundary condition?
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Convection terms default to a no flux boundary condition unless the
-exterior faces are associated with a constraint, in which case either
-an inlet or an outlet boundary condition is applied depending on the
-flow direction.
-
-How do I apply spatially varying boundary conditions?
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-The use of spatial varying boundary conditions is best demonstrated with an
-example. Given a 2D equation in the domain :math:`0 < x < 1` and :math:`0 < y < 1` with
-boundary conditions,
-
-.. math::
-
-  \phi = \left\{
-            \begin{aligned} 
-                xy &\quad \text{on $x>1/2$ and $y>1/2$} \\
-                \vec{n} \cdot \vec{F} = 0 &\quad \text{elsewhere}
-            \end{aligned}
-        \right.
-
-where :math:`\vec{F}` represents the flux. The boundary conditions in :term:`FiPy` can
-be written with the following code,
-
->>> x, y = mesh.faceCenters
->>> mask =  ((x < 0.5) | (y < 0.5))
->>> var.faceGrad.constrain(0, where=mesh.exteriorFaces & mask)
->>> var.constrain(x * y, where=mesh.exteriorFaces & ~mask)
-
-then
-
->>> eqn.solve(...)
-
-Further demonstrations of spatially varying boundary condition can be found
-in :mod:`examples.diffusion.mesh20x20`
-and :mod:`examples.diffusion.circle`
-
-.. %    http://thread.gmane.org/gmane.comp.python.fipy/726
-   %    http://thread.gmane.org/gmane.comp.python.fipy/846
-
-   %    \subsection{Fourth order boundary conditions}
-
-   %    http://thread.gmane.org/gmane.comp.python.fipy/923
-
-   %    \subsection{Periodic boundary conditions}
-
-   %    http://thread.gmane.org/gmane.comp.python.fipy/135
-
-   %    \subsection{Time dependent boundary conditions}
-
-   %    http://thread.gmane.org/gmane.comp.python.fipy/2
-
-   %    \subsection{Internal boundary conditions}
+See the :ref:`BoundaryConditions` section for more details.
 
 What does this error message mean?
 ----------------------------------
diff --git a/documentation/USAGE.txt b/documentation/USAGE.txt
index c1ef4c9c0c..c6942b68d1 100644
--- a/documentation/USAGE.txt
+++ b/documentation/USAGE.txt
@@ -223,7 +223,7 @@ Solving in Parallel
 
 :term:`FiPy` can use :term:`Trilinos` to solve equations in
 parallel. Most mesh classes in :mod:`fipy.meshes` can solve in
-paralled. This includes all "``...Grid...``" and "``...Gmsh...``"
+parallel. This includes all "``...Grid...``" and "``...Gmsh...``"
 class meshes. Currently, the only remaining serial-only meshes are
 :class:`~fipy.meshes.tri2D.Tri2D` and
 :class:`~fipy.meshes.skewedGrid2D.SkewedGrid2D`.
@@ -324,7 +324,7 @@ you need to do something with the entire solution, you can use
         Different solvers, different preconditioners, or a less restrictive
         tolerance may help.
 
-.. _RunningUnderPython3:
+.. _MeshingWithGmsh:
 
 -----------------
 Meshing with Gmsh
@@ -353,6 +353,8 @@ parallel including :class:`~fipy.meshes.gmshImport.Gmsh2D` and
     :term:`FiPy` solution accuracy can be compromised with highly
     non-orthogonal or non-conjunctional meshes.
 
+.. _CoupledEquations:
+
 ----------------------------
 Coupled and Vector Equations
 ----------------------------
@@ -360,13 +362,248 @@ Coupled and Vector Equations
 Equations can now be coupled together so that the contributions from
 all the equations appear in a single system matrix. This results in
 tighter coupling for equations with spatial and temporal derivatives
-in more than one variable. The use of coupled equation is described in
-:mod:`examples.diffusin.coupled`. Other examples that demonstrate the
-use of coupled equations are :mod:`examples.phase.binaryCoupled` and
+in more than one variable. In :term:`FiPy` equations are coupled
+together using the ``&`` operator::
+
+   >>> eqn0 = ...
+   >>> eqn1 = ...
+   >>> coupledEqn = eqn0 & eqn1
+
+The ``coupledEqn`` will use a combined system matrix that includes
+four quadrants for each of the different variable and equation
+combinations. In previous versions of :term:`FiPy` there has been no
+need to specify which variable a given term acts on when generating
+equations. The variable is simply specified when calling ``solve`` or
+``sweep`` and this functionality has been maintained in the case of
+single equations. However, for coupled equations the variable that a
+given term operates on now needs to be specified when the equation is
+generated. The syntax for generating coupled equations has the form::
+
+   >>> eqn0 = Term00(coeff=..., var=var0) + Term01(coeff..., var=var1) == source0
+   >>> eqn1 = Term10(coeff=..., var=var0) + Term11(coeff..., var=var1) == source1
+   >>> coupledEqn = eqn0 & eqn1
+
+and there is now no need to pass any variables when solving::
+
+   >>> coupledEqn.solve(dt=..., solver=...)
+
+In this case the matrix system will have the form
+
+.. math::
+
+   \left(
+   \begin{array}{c|c}
+   \text{\ttfamily Term00} & \text{\ttfamily Term01} \\ \hline
+   \text{\ttfamily Term10} & \text{\ttfamily Term11}
+   \end{array} \right)
+   \left(
+   \begin{array}{c}
+   \text{\ttfamily var0}  \\ \hline
+   \text{\ttfamily var1} 
+   \end{array} \right)
+   =	       
+   \left(
+   \begin{array}{c}
+   \text{\ttfamily source0}  \\ \hline
+   \text{\ttfamily source1} 
+   \end{array} \right)
+   
+:term:`FiPy` tries to make sensible decisions regarding each term's
+location in the matrix and the ordering of the variable column
+array. For example, if ``Term01`` is a transient term then ``Term01``
+would appear in the upper left diagonal and the ordering of the
+variable column array would be reversed.
+
+The use of coupled equation is described in detail in
+:mod:`examples.diffusion.coupled`. Other examples that demonstrate the
+use of coupled equations are :mod:`examples.phase.binaryCoupled`,
+:mod:`examples.phase.polyxtalCoupled` and
 :mod:`examples.cahnHilliard.mesh2DCoupled`. As well as coupling
 equations, true vector equations can now be written in :term:`FiPy`
 (see :mod:`examples.diffusion.coupled` for more details).
 
+.. _BoundaryConditions:
+
+-------------------
+Boundary Conditions
+-------------------
+
+.. currentmodule:: fipy.variables.cellVariable
+
+Applying fixed value (Dirichlet) boundary conditions
+====================================================
+
+To apply a fixed value boundary condition use the
+:meth:`~CellVariable.constrain` method. For example, to fix `var` to
+have a value of `2` along the upper surface of a domain, use
+
+>>> var.constrain(2., where=mesh.facesTop)
+
+.. note:: 
+
+   The old equivalent
+   :class:`~fipy.boundaryConditions.fixedValue.FixedValue` boundary
+   condition is now deprecated.
+
+Applying fixed gradient boundary conditions (Neumann)
+=====================================================
+
+To apply a fixed Gradient boundary condition use the
+:attr:`~.CellVariable.faceGrad`.\
+:meth:`~fipy.variables.variable.Variable.constrain` method. For
+example, to fix `var` to have a gradient of `(0,2)` along the upper
+surface of a 2D domain, use
+
+>>> var.faceGrad.constrain(((0,),(2,)), where=mesh.facesTop)
+
+Applying fixed flux boundary conditions
+=======================================
+
+Generally these can be implemented with a judicious use of
+:attr:`~.CellVariable.faceGrad`.\
+:meth:`~fipy.variables.variable.Variable.constrain`.  Failing that, an
+exterior flux term can be added to the equation. Firstly, set the
+terms' coefficients to be zero on the exterior faces,
+
+>>> diffCoeff.constrain(0., mesh.exteriorFaces)
+>>> convCoeff.constrain(0., mesh.exteriorFaces)
+
+then create an equation with an extra term to account for the exterior flux,
+
+>>> eqn = (TransientTerm() + ConvectionTerm(convCoeff) 
+...        == DiffusionCoeff(diffCoeff)
+...        + (mesh.exteriorFaces * exteriorFlux).divergence)
+
+where `exteriorFlux` is a rank 1
+:class:`~fipy.variables.faceVariable.FaceVariable`.
+
+.. note:: 
+
+   The old equivalent :class:`~fipy.boundaryConditions.fixedFlux.FixedFlux`
+   boundary condition is now deprecated.
+
+Applying outlet or inlet boundary conditions
+============================================
+
+Convection terms default to a no flux boundary condition unless the
+exterior faces are associated with a constraint, in which case either
+an inlet or an outlet boundary condition is applied depending on the
+flow direction.
+
+Applying spatially varying boundary conditions
+==============================================
+
+The use of spatial varying boundary conditions is best demonstrated with an
+example. Given a 2D equation in the domain :math:`0 < x < 1` and :math:`0 < y < 1` with
+boundary conditions,
+
+.. math::
+
+  \phi = \left\{
+            \begin{aligned} 
+                xy &\quad \text{on $x>1/2$ and $y>1/2$} \\
+                \vec{n} \cdot \vec{F} = 0 &\quad \text{elsewhere}
+            \end{aligned}
+        \right.
+
+where :math:`\vec{F}` represents the flux. The boundary conditions in :term:`FiPy` can
+be written with the following code,
+
+>>> X, Y = mesh.faceCenters
+>>> mask =  ((X < 0.5) | (Y < 0.5))
+>>> var.faceGrad.constrain(0, where=mesh.exteriorFaces & mask)
+>>> var.constrain(X * Y, where=mesh.exteriorFaces & ~mask)
+
+then
+
+>>> eqn.solve(...)
+
+Further demonstrations of spatially varying boundary condition can be found
+in :mod:`examples.diffusion.mesh20x20`
+and :mod:`examples.diffusion.circle`
+
+Applying internal boundary conditions
+=====================================
+
+Applying internal boundary conditions can be achieved through the use
+of implicit and explicit sources. An equation of the form
+
+>>> eqn = TransientTerm() == DiffusionTerm()
+
+can be constrained to have a fixed internal ``value`` at a position
+given by ``mask`` with the following alterations
+
+>>> eqn = TransientTerm() == DiffusionTerm() - ImplicitSourceTerm(mask * largeValue) + mask * largeValue * value 
+
+The parameter ``largeValue`` must be chosen to be large enough to
+completely dominate the matrix diagonal and the RHS vector in cells
+that are masked. The ``mask`` variable would typically be a
+``CellVariable`` boolean constructed using the cell center values.
+
+One must be careful to distinguish between constraining internal cell
+values during the solve step and simply applying arbitrary constraints
+to a ``CellVariable``. Applying a constraint,
+
+>>> var.constrain(value, where=mask)
+
+simply fixes the returned value of ``var`` at ``mask`` to be
+``value``. It does not have any effect on the implicit value of ``var`` at the
+``mask`` location during the linear solve so it is not a substitute
+for the source term machinations described above. Future releases of
+:term:`FiPy` may implicitly deal with this discrepancy, but the current
+release does not. A simple example can be used to demonstrate this::
+
+>>> m = Grid1D(nx=2, dx=1.)
+>>> var = CellVariable(mesh=m)
+
+Apply a constraint to the faces for a right side boundary condition
+(which works).
+
+>>> var.constrain(1., where=m.facesRight)
+
+Create the equation with the source term constraint described above
+
+>>> mask = m.x < 1.
+>>> largeValue = 1e+10
+>>> value = 0.25
+>>> eqn = DiffusionTerm() - ImplicitSourceTerm(largeValue * mask) + largeValue * mask * value
+
+and the expected value is obtained.
+
+>>> eqn.solve(var)
+>>> print var
+[ 0.25  0.75]
+
+However, if a constraint is used without the source term constraint an
+unexpected value is obtained
+
+>>> var.constrain(0.25, where=mask)
+>>> eqn = DiffusionTerm()
+>>> eqn.solve(var)
+>>> print var
+[ 0.25  1.  ]
+
+although the left cell has the expected value as it is constrained.
+
+.. %    http://thread.gmane.org/gmane.comp.python.fipy/726
+   %    http://thread.gmane.org/gmane.comp.python.fipy/846
+
+   %    \subsection{Fourth order boundary conditions}
+
+   %    http://thread.gmane.org/gmane.comp.python.fipy/923
+
+   %    \subsection{Periodic boundary conditions}
+
+   %    http://thread.gmane.org/gmane.comp.python.fipy/135
+
+   %    \subsection{Time dependent boundary conditions}
+
+   %    http://thread.gmane.org/gmane.comp.python.fipy/2
+
+   %    \subsection{Internal boundary conditions}
+
+.. _RunningUnderPython3:
+
 ----------------------
 Running under Python 3     
 ----------------------
diff --git a/documentation/numerical/discret.txt b/documentation/numerical/discret.txt
index 41a41fa9f5..34dc5545bb 100644
--- a/documentation/numerical/discret.txt
+++ b/documentation/numerical/discret.txt
@@ -182,14 +182,10 @@ This term is represented in :term:`FiPy` as
 
 >>> DiffusionTerm(coeff=Gamma1)
 
-which is synonymous with
-
->>> ImplicitDiffusionTerm(coeff=Gamma1)
-    
 or 
 
 >>> ExplicitDiffusionTerm(coeff=Gamma1)
-    
+
 :class:`~explicitDiffusionTerm.ExplicitDiffusionTerm` is provided primarily for
 illustrative purposes, although :mod:`examples.diffusion.mesh1D` 
 demonstrates its use in Crank-Nicolson time stepping.
@@ -272,10 +268,9 @@ equation to a set of discrete linear equations that can then be solved
 to obtain the value of the dependent variable at each CV center. This
 results in a sparse linear system that requires an efficient iterative
 scheme to solve. The iterative schemes available to :term:`FiPy` are
-currently encapsulated in the :term:`PySparse` suite of solvers
-and include most common solvers such as the conjugate gradient method
-and LU decomposition. There are plans to include other solver suites
-that are compatible with :term:`Python`.
+currently encapsulated in the :term:`PySparse` and :term:`PyTrilinos` 
+suites of solvers and include most common solvers such as the conjugate
+gradient method and LU decomposition.
 
 Combining Equations :eq:`num:tra`, :eq:`num:con`,
 :eq:`num:dif` and :eq:`num:sou`, the complete
diff --git a/documentation/numerical/equation.txt b/documentation/numerical/equation.txt
index da213e558d..e3356b48dc 100644
--- a/documentation/numerical/equation.txt
+++ b/documentation/numerical/equation.txt
@@ -21,7 +21,7 @@ Examples of such systems are wide ranging, but include problems that
 exhibit a combination of diffusing and reacting species, as well as
 such diverse problems as determination of the electric potential in
 heart tissue, of fluid flow, stress evolution, and even the
-Schr\"odinger equation.
+Schroedinger equation.
 
 A general conservation equation, solved using :term:`FiPy`, can include any
 combination of the following terms,
diff --git a/documentation/numerical/index.txt b/documentation/numerical/index.txt
index 7253e0a763..71193e9cd5 100644
--- a/documentation/numerical/index.txt
+++ b/documentation/numerical/index.txt
@@ -24,7 +24,7 @@ The FVM can be thought of as a subset of the Finite Element Method
 system of equations fully equivalent to the FVM can be obtained with
 the FEM using as weighting functions the characteristic functions of
 FV cells, i.e., functions equal to unity [Mattiussi:1997]_. Analogously,
-the the discretization of equations with the FVM reduces to the FDM on
+the discretization of equations with the FVM reduces to the FDM on
 Cartesian grids.
 
 .. toctree::
diff --git a/documentation/tutorial/index.txt b/documentation/tutorial/index.txt
index 07e3c4814f..f22a483a16 100644
--- a/documentation/tutorial/index.txt
+++ b/documentation/tutorial/index.txt
@@ -14,5 +14,5 @@ for their containing package.
 .. toctree::
    :maxdepth: 4
 
-   package/generated/subpackage
+   package/generated/package.subpackage
 
diff --git a/examples/cahnHilliard/mesh2DCoupled.py b/examples/cahnHilliard/mesh2DCoupled.py
index 8767ccd13a..f917d42759 100755
--- a/examples/cahnHilliard/mesh2DCoupled.py
+++ b/examples/cahnHilliard/mesh2DCoupled.py
@@ -31,7 +31,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the Cahn-Hilliard problem in two dimensions.
+
 The spinodal decomposition phenomenon is a spontaneous separation of
 an initially homogenous mixture into two distinct regions of different
 properties (spin-up/spin-down, component A/component B). It is a
@@ -171,12 +172,16 @@
 has a shape `(2,)`
 
 >>> source = (- d2fdphi2 * v0 + dfdphi) * (0, 1)
->>> impCoeff = -d2fdphi2 * ((0, 0), (1., 0)) + ((0, 0), (0, -1.))
+>>> impCoeff = -d2fdphi2 * ((0, 0), 
+...                         (1., 0)) + ((0, 0), 
+...                                     (0, -1.))
 
 This is the same equation as the previous definition of `eq`, but now in
 a vector format.
 
->>> eq = TransientTerm(((1., 0.), (0., 0.))) == DiffusionTerm([((0., D), (-epsilon**2, 0.))]) + ImplicitSourceTerm(impCoeff) + source
+>>> eq = TransientTerm(((1., 0.), 
+...                     (0., 0.))) == DiffusionTerm([((0.,          D), 
+...                                                   (-epsilon**2, 0.))]) + ImplicitSourceTerm(impCoeff) + source
 
 >>> dexp = -5
 >>> elapsed = 0.
diff --git a/examples/cahnHilliard/sphere.py b/examples/cahnHilliard/sphere.py
index f4e7fb19e7..116388f9aa 100755
--- a/examples/cahnHilliard/sphere.py
+++ b/examples/cahnHilliard/sphere.py
@@ -31,9 +31,9 @@
  # ###################################################################
  ##
 
-r""" 
-Solves the Cahn-Hilliard problem on the surface of a sphere, such as
-may occur on vesicles (http://www.youtube.com/watch?v=kDsFP67_ZSE).
+r"""Solves the Cahn-Hilliard problem on the surface of a sphere.
+
+This phenomenon canoccur on vesicles (http://www.youtube.com/watch?v=kDsFP67_ZSE).
 
 >>> from fipy import *
 
diff --git a/examples/convection/exponential1D/mesh1D.py b/examples/convection/exponential1D/mesh1D.py
index fd3f8d36ae..dc1a566860 100755
--- a/examples/convection/exponential1D/mesh1D.py
+++ b/examples/convection/exponential1D/mesh1D.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the steady-state convection-diffusion equation in one dimension.
 
 This example solves the steady-state convection-diffusion equation
 given by
@@ -41,7 +41,7 @@
 
    \nabla \cdot \left(D \nabla \phi + \vec{u} \phi \right) = 0
 
-with coefficients :math:`D = 1` and :math:`\vec{u} = (10,)`, or
+with coefficients :math:`D = 1` and :math:`\vec{u} = 10\hat{\i}`, or
 
 >>> diffCoeff = 1.
 >>> convCoeff = (10.,)
@@ -61,7 +61,7 @@
 
 The solution variable is initialized to ``valueLeft``:
 
->>> var = CellVariable(mesh=mesh, name = "variable")
+>>> var = CellVariable(mesh=mesh, name="variable")
 
 and impose the boundary conditions
 
diff --git a/examples/convection/exponential1DSource/mesh1D.py b/examples/convection/exponential1DSource/mesh1D.py
index 94dbf34907..dc0e0e382f 100755
--- a/examples/convection/exponential1DSource/mesh1D.py
+++ b/examples/convection/exponential1DSource/mesh1D.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the steady-state convection-diffusion equation with a constant source.
 
 Like :mod:`examples.convection.exponential1D.mesh1D`
 this example solves a steady-state convection-diffusion equation, but adds a constant source, 
diff --git a/examples/convection/robin.py b/examples/convection/robin.py
index dcb9d16591..dc1ac6803e 100644
--- a/examples/convection/robin.py
+++ b/examples/convection/robin.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve an advection-diffusion equation with a Robin boundary condition.
 
 This example demonstrates how to apply a Robin boundary condition to
 an advection-diffusion equation. The equation we wish to solve is
diff --git a/examples/convection/source.py b/examples/convection/source.py
index 705e6db73d..9b5ab4f7c4 100644
--- a/examples/convection/source.py
+++ b/examples/convection/source.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a convection problem with a source.
+
 This example solves the equation
 
 .. math::
@@ -44,7 +45,7 @@
 implementation of an outflow boundary condition, which is not
 currently implemented in FiPy. An :class:`~fipy.terms.implicitSourceTerm.ImplicitSourceTerm` object
 will be used to represent this term. The derivative of :math:`\phi` can be
-represented by a :class:`~fipy.terms.abstractConvectionTerm._AbstractConvectionTerm` with a constant unitary velocity
+represented by a :class:`~fipy.terms.ConvectionTerm` with a constant unitary velocity
 field from left to right. The following is an example code that includes
 a test against the analytical result.
 
diff --git a/examples/diffusion/anisotropy.py b/examples/diffusion/anisotropy.py
index afbf042b57..df00d72820 100644
--- a/examples/diffusion/anisotropy.py
+++ b/examples/diffusion/anisotropy.py
@@ -32,10 +32,9 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the diffusion equation with an anisotropic diffusion coefficient.
 
-This example demonstrates how to solve diffusion with an anisotropic
-coefficient.  We wish to solve the problem
+We wish to solve the problem
 
 .. math::
 
@@ -86,7 +85,7 @@
 >>> import os
 >>> mesh = Gmsh2D(os.path.splitext(__file__)[0] + '.msh', communicator=serial) # doctest: +GMSH
 
-Set the center most cell to have a value.
+Set the centermost cell to have a value.
 
 >>> var = CellVariable(mesh=mesh, hasOld=1) # doctest: +GMSH
 >>> x, y = mesh.cellCenters # doctest: +GMSH
diff --git a/examples/diffusion/circle.py b/examples/diffusion/circle.py
index db73a9888d..2497285a71 100755
--- a/examples/diffusion/circle.py
+++ b/examples/diffusion/circle.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the diffusion equation in a circular domain meshed with triangles.
+
 This example demonstrates how to solve a simple diffusion problem on a
 non-standard mesh with varying boundary conditions. The :term:`Gmsh` package
 is used to create the mesh. Firstly, define some parameters for the
@@ -54,20 +55,20 @@
 
 >>> from fipy import *
 >>> mesh = Gmsh2D('''
-...                       cellSize = %(cellSize)g;
-...                       radius = %(radius)g;
-...                       Point(1) = {0, 0, 0, cellSize};
-...                       Point(2) = {-radius, 0, 0, cellSize};
-...                       Point(3) = {0, radius, 0, cellSize};
-...                       Point(4) = {radius, 0, 0, cellSize};
-...                       Point(5) = {0, -radius, 0, cellSize};
-...                       Circle(6) = {2, 1, 3};
-...                       Circle(7) = {3, 1, 4};
-...                       Circle(8) = {4, 1, 5};
-...                       Circle(9) = {5, 1, 2};
-...                       Line Loop(10) = {6, 7, 8, 9};
-...                       Plane Surface(11) = {10};
-...                       ''' % locals()) # doctest: +GMSH
+...               cellSize = %(cellSize)g;
+...               radius = %(radius)g;
+...               Point(1) = {0, 0, 0, cellSize};
+...               Point(2) = {-radius, 0, 0, cellSize};
+...               Point(3) = {0, radius, 0, cellSize};
+...               Point(4) = {radius, 0, 0, cellSize};
+...               Point(5) = {0, -radius, 0, cellSize};
+...               Circle(6) = {2, 1, 3};
+...               Circle(7) = {3, 1, 4};
+...               Circle(8) = {4, 1, 5};
+...               Circle(9) = {5, 1, 2};
+...               Line Loop(10) = {6, 7, 8, 9};
+...               Plane Surface(11) = {10};
+...               ''' % locals()) # doctest: +GMSH
 
 Using this mesh, we can construct a solution variable
 
diff --git a/examples/diffusion/circleQuad.py b/examples/diffusion/circleQuad.py
index dd221a88bd..c0f80e3685 100755
--- a/examples/diffusion/circleQuad.py
+++ b/examples/diffusion/circleQuad.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the diffusion equation in a circular domain meshed with quadrangles.
+
 This example demonstrates how to solve a simple diffusion problem on a
 non-standard mesh with varying boundary conditions. The :term:`Gmsh` package
 is used to create the mesh. Firstly, define some parameters for the
diff --git a/examples/diffusion/coupled.py b/examples/diffusion/coupled.py
index 8ea65778c5..6791caab70 100644
--- a/examples/diffusion/coupled.py
+++ b/examples/diffusion/coupled.py
@@ -32,78 +32,125 @@
  # ###################################################################
  ##
 
-r"""
-How to use coupled equations.
+r"""Solve the biharmonic equation as a coupled pair of diffusion equations.
 
+:term:`FiPy` has only first order time derivatives so equations such
+as the biharmonic wave equation written as
 
-## Coupled
-
-
-from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
-
-m = Grid1D(nx=100, Lx=1.)
-
-v0 = CellVariable(mesh=m, hasOld=True, value=0.5)
-v1 = CellVariable(mesh=m, hasOld=True, value=0.5)
-
-v0.constrain(0, m.facesLeft)
-v0.constrain(1, m.facesRight)
-
-v1.constrain(1, m.facesLeft)
-v1.constrain(0, m.facesRight)
-
-eqn0 = TransientTerm(var=v0) == DiffusionTerm(-1, var=v1) + DiffusionTerm(0.01, var=v0)
-eqn1 = TransientTerm(var=v1) == DiffusionTerm(1, var=v0) + DiffusionTerm(0.01, var=v1)
-
-eqn = eqn0 & eqn1
-
-vi = Viewer((v0, v1))
-
-for t in range(1): 
-    v0.updateOld()
-    v1.updateOld()
-    eqn.solve(dt=1.)
-    vi.plot()
-
-## Vector
-
-v = CellVariable(mesh=m, hasOld=True, value=0.5, elementshape=(2,))
-
-v.constrain(0, [m.facesLeft, m.facesRight])
-v.constrain(1, [m.facesRight, m.facesLeft])
-
-eqn = TransientTerm(((1, 0), (0, 1))) == DiffusionTerm([((0.01, -1), (1, 0.01))])
-
-vi = Viewer((v[0], v[1]))
-
-for t in range(1): 
-    v.updateOld()
-    eqn.solve(var=v, dt=1.)
-    vi.plot()
-
-## Uncoupled
-
-m = Grid1D(nx=100, Lx=1.)
-
-v0 = CellVariable(mesh=m, hasOld=True, value=0.5)
-v1 = CellVariable(mesh=m, hasOld=True, value=0.5)
-
-v0.constrain(0, m.facesLeft)
-v0.constrain(1, m.facesRight)
-
-v1.constrain(1, m.facesLeft)
-v1.constrain(0, m.facesRight)
-
-eq0 = TransientTerm() == -v1.faceGrad.divergence + DiffusionTerm(0.01)
-eq1 = TransientTerm() == v0.faceGrad.divergence + DiffusionTerm(0.01)
+.. math::
+    
+   \frac{\partial^4 v}{\partial x^4} + \frac{\partial^2 v}{\partial t^2} &= 0
+   
+cannot be represented as a single equation. We need to decompose the
+biharmonic equation into two equations that are first order in time in
+the following way,
+
+.. math::
+    
+   \frac{\partial^2 v_0}{\partial x^2} + \frac{\partial v_1}{\partial t} &= 0 \\
+   \frac{\partial^2 v_1}{\partial x^2} - \frac{\partial v_0}{\partial t} &= 0 
+
+Historically, :term:`FiPy` required systems of coupled equations to be
+solved successively, "sweeping" the equations to convergence. As a
+practical example, we use the following system
+
+.. math::
+    
+   \frac{\partial v_0}{\partial t} &= 0.01 \nabla^2 v_0 - \nabla^2 v_1 \\
+   \frac{\partial v_1}{\partial t} &= \nabla^2 v_0 + 0.01 \nabla^2 v_1
+   
+subject to the boundary conditions 
+
+.. math::
+   :nowrap:
+    
+   \begin{align*}
+   v_0|_{x=0} &= 0 & v_0|_{x=1} &= 1 \\
+   v_1|_{x=0} &= 1 & v_1|_{x=1} &= 0
+   \end{align*}
+
+This system closely resembles the pure biharmonic equation, but has an
+additional diffusion contribution to improve numerical stability.  The
+example system is solved with the following block of code using
+explicit coupling for the cross-coupled terms.
+
+>>> from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
+
+>>> m = Grid1D(nx=100, Lx=1.)
+
+>>> v0 = CellVariable(mesh=m, hasOld=True, value=0.5)
+>>> v1 = CellVariable(mesh=m, hasOld=True, value=0.5)
+
+>>> v0.constrain(0, m.facesLeft)
+>>> v0.constrain(1, m.facesRight)
+
+>>> v1.constrain(1, m.facesLeft)
+>>> v1.constrain(0, m.facesRight)
+
+>>> eq0 = TransientTerm() == DiffusionTerm(coeff=0.01) - v1.faceGrad.divergence
+>>> eq1 = TransientTerm() == v0.faceGrad.divergence + DiffusionTerm(coeff=0.01)
+
+>>> vi = Viewer((v0, v1))
+
+>>> for t in range(100): 
+...     v0.updateOld()
+...     v1.updateOld()
+...     res0 = res1 = 1e100
+...     while max(res0, res1) > 0.1:
+...         res0 = eq0.sweep(var=v0, dt=1e-5)
+...         res1 = eq1.sweep(var=v1, dt=1e-5)
+...     if t % 10 == 0:
+...         vi.plot()
+
+The uncoupled method still works, but it can be advantageous to solve
+the two equations simultaneously. In this case, by coupling the
+equations, we can eliminate the explicit sources and dramatically
+increase the time steps:
+
+>>> v0.value = 0.5
+>>> v1.value = 0.5
+
+>>> eqn0 = TransientTerm(var=v0) == DiffusionTerm(0.01, var=v0) - DiffusionTerm(1, var=v1)
+>>> eqn1 = TransientTerm(var=v1) == DiffusionTerm(1, var=v0) + DiffusionTerm(0.01, var=v1)
+
+>>> eqn = eqn0 & eqn1
+
+>>> for t in range(1): 
+...     v0.updateOld()
+...     v1.updateOld()
+...     eqn.solve(dt=1.e-3)
+...     vi.plot()
+
+It is also possible to pose the same equations in vector form:
+
+>>> v = CellVariable(mesh=m, hasOld=True, value=[[0.5], [0.5]], elementshape=(2,))
+
+>>> v.constrain([[0], [1]], m.facesLeft)
+>>> v.constrain([[1], [0]], m.facesRight)
+
+>>> eqn = TransientTerm([[1, 0], 
+...                      [0, 1]]) == DiffusionTerm([[[0.01, -1], 
+...                                                  [1, 0.01]]])
+
+>>> vi = Viewer((v[0], v[1]))
+
+>>> for t in range(1): 
+...     v.updateOld()
+...     eqn.solve(var=v, dt=1.e-3)
+...     vi.plot()
+
+Whether you pose your problem in coupled or vector form should be dictated by
+the underlying physics. If :math:`v_0` and :math:`v_1` represent the
+concentrations of two conserved species, then it is natural to write two
+seperate governing equations and to couple them. If they represent two
+components of a vector field, then the vector formulation is obviously more
+natural. FiPy will solve the same matrix system either way.
+"""
+__docformat__ = 'restructuredtext'
 
-vi = Viewer((v0, v1))
+if __name__ == '__main__':
+    import fipy.tests.doctestPlus
+    exec(fipy.tests.doctestPlus._getScript())
 
-for t in range(10000): 
-    v0.updateOld()
-    v1.updateOld()
-    eq0.solve(var=v0, dt=1e-5)
-    eq1.solve(var=v1, dt=1e-5)
-    vi.plot()
+    raw_input('finished')
 
-"""
diff --git a/examples/diffusion/electrostatics.py b/examples/diffusion/electrostatics.py
index 8b9321a26b..a4fac62047 100755
--- a/examples/diffusion/electrostatics.py
+++ b/examples/diffusion/electrostatics.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the Poisson equation in one dimension.
+
 The Poisson equation is a particular example of the steady-state diffusion
 equation. We examine a few cases in one dimension.
 
diff --git a/examples/diffusion/mesh1D.py b/examples/diffusion/mesh1D.py
index adc5ad698d..497593d777 100755
--- a/examples/diffusion/mesh1D.py
+++ b/examples/diffusion/mesh1D.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a one-dimensional diffusion equation under different conditions.
+
 To run this example from the base :term:`FiPy` directory, type::
     
     $ python examples/diffusion/mesh1D.py
@@ -57,7 +58,7 @@
 
 >>> nx = 50
 >>> dx = 1.
->>> mesh = Grid1D(nx = nx, dx = dx)
+>>> mesh = Grid1D(nx=nx, dx=dx)
 
 :term:`FiPy` solves all equations at the centers of the cells of the mesh. We
 thus need a :class:`~fipy.variables.cellVariable.CellVariable` object to hold the values of the
@@ -389,8 +390,8 @@
 .. index:: FaceVariable
 
 >>> D = FaceVariable(mesh=mesh, value=1.0)
->>> x = mesh.faceCenters[0]
->>> D.setValue(0.1, where=(L / 4. <= x) & (x < 3. * L / 4.))
+>>> X = mesh.faceCenters[0]
+>>> D.setValue(0.1, where=(L / 4. <= X) & (X < 3. * L / 4.))
 
 The boundary conditions are a fixed value of 
 
diff --git a/examples/diffusion/mesh20x20.py b/examples/diffusion/mesh20x20.py
index 6d79f8a125..d81acbdb36 100755
--- a/examples/diffusion/mesh20x20.py
+++ b/examples/diffusion/mesh20x20.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a two-dimensional diffusion problem in a square domain.
 
 This example solves a diffusion problem and demonstrates the use of
 applying boundary condition patches.
@@ -68,11 +68,11 @@
 to the top-left and bottom-right corners.  Neumann boundary conditions
 are automatically applied to the top-right and bottom-left corners.
 
->>> x, y = mesh.faceCenters
->>> facesTopLeft = ((mesh.facesLeft & (y > L / 2))
-...                 | (mesh.facesTop & (x < L / 2)))
->>> facesBottomRight = ((mesh.facesRight & (y < L / 2))
-...                     | (mesh.facesBottom & (x > L / 2)))
+>>> X, Y = mesh.faceCenters
+>>> facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
+...                 | (mesh.facesTop & (X < L / 2)))
+>>> facesBottomRight = ((mesh.facesRight & (Y < L / 2))
+...                     | (mesh.facesBottom & (X > L / 2)))
 
 >>> phi.constrain(valueTopLeft, facesTopLeft)
 >>> phi.constrain(valueBottomRight, facesBottomRight)
diff --git a/examples/diffusion/mesh20x20Coupled.py b/examples/diffusion/mesh20x20Coupled.py
index 49ee67f2bc..1eb2ff18a1 100755
--- a/examples/diffusion/mesh20x20Coupled.py
+++ b/examples/diffusion/mesh20x20Coupled.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a coupled set of diffusion equations in two dimensions.
 
 This example solves a diffusion problem and demonstrates the use of
 applying boundary condition patches.
diff --git a/examples/diffusion/nthOrder/input4thOrder1D.py b/examples/diffusion/nthOrder/input4thOrder1D.py
index fd4aa8a80b..ccac613af0 100755
--- a/examples/diffusion/nthOrder/input4thOrder1D.py
+++ b/examples/diffusion/nthOrder/input4thOrder1D.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a fourth-order diffusion problem.
 
 This example uses the :class:`~fipy.terms.diffusionTerm.DiffusionTerm` class to solve the equation
 
diff --git a/examples/flow/stokesCavity.py b/examples/flow/stokesCavity.py
index a5f2ada8c1..23bad3d406 100755
--- a/examples/flow/stokesCavity.py
+++ b/examples/flow/stokesCavity.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the Navier-Stokes equation in the viscous limit.
 
 Many thanks to Benny Malengier <bm@cage.ugent.be> for reworking this example and
 actually making it work correctly...see changeset:3799
diff --git a/examples/levelSet/advection/circle.py b/examples/levelSet/advection/circle.py
index 7593cbe45c..da3a32de4b 100755
--- a/examples/levelSet/advection/circle.py
+++ b/examples/levelSet/advection/circle.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a circular distance function equation and then advect it.
+
 This example first imposes a circular distance function:
 
 .. math::
diff --git a/examples/levelSet/advection/mesh1D.py b/examples/levelSet/advection/mesh1D.py
index 3781fb9976..18a9e3d6d2 100755
--- a/examples/levelSet/advection/mesh1D.py
+++ b/examples/levelSet/advection/mesh1D.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the distance function equation in one dimension and then advect it.
+
 This example first solves the distance function equation in one dimension:
 
 .. math::
diff --git a/examples/levelSet/distanceFunction/circle.py b/examples/levelSet/distanceFunction/circle.py
index ec0fbcc5f6..03f81db867 100755
--- a/examples/levelSet/distanceFunction/circle.py
+++ b/examples/levelSet/distanceFunction/circle.py
@@ -32,10 +32,9 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve the level set equation in two dimensions for a circle. 
 
-Here we solve the level set equation in two dimensions for a circle. The 2D
-level set equation can be written,
+The 2D level set equation can be written,
 
 .. math::
 
diff --git a/examples/levelSet/distanceFunction/mesh1D.py b/examples/levelSet/distanceFunction/mesh1D.py
index 2926429264..ebd7bd4f57 100755
--- a/examples/levelSet/distanceFunction/mesh1D.py
+++ b/examples/levelSet/distanceFunction/mesh1D.py
@@ -32,9 +32,9 @@
  # ###################################################################
  ##
 
-r"""
+r"""Create a level set variable in one dimension. 
 
-Here we create a level set variable in one dimension. The level set
+The level set
 variable calculates its value over the domain to be the distance from
 the zero level set. This can be represented succinctly in the
 following equation with a boundary condition at the zero level set
diff --git a/examples/levelSet/electroChem/gold.py b/examples/levelSet/electroChem/gold.py
index 67f78943ee..b371e47339 100644
--- a/examples/levelSet/electroChem/gold.py
+++ b/examples/levelSet/electroChem/gold.py
@@ -31,7 +31,8 @@
  # ###################################################################
  ##
  
-r"""
+r"""Model electrochemical superfill of gold using the CEAC mechanism.
+
 This input file
 is a demonstration of the use of :term:`FiPy` for
 modeling gold superfill. The material properties and experimental
diff --git a/examples/levelSet/electroChem/howToWriteAScript.py b/examples/levelSet/electroChem/howToWriteAScript.py
index c6968f820e..e18d8d4810 100644
--- a/examples/levelSet/electroChem/howToWriteAScript.py
+++ b/examples/levelSet/electroChem/howToWriteAScript.py
@@ -31,7 +31,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Tutorial for writing an electrochemical superfill script.
 
 This input file demonstrates how to
 create a new superfill script if the existing suite of scripts do
diff --git a/examples/levelSet/electroChem/leveler.py b/examples/levelSet/electroChem/leveler.py
index 3c98c420e0..57263fc115 100755
--- a/examples/levelSet/electroChem/leveler.py
+++ b/examples/levelSet/electroChem/leveler.py
@@ -31,7 +31,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Model electrochemical superfill of copper with leveler and accelerator additives.
+
 This input file is a demonstration of the use of :term:`FiPy` for
 modeling copper superfill with leveler and accelerator
 additives. The material properties and experimental parameters
diff --git a/examples/levelSet/electroChem/simpleTrenchSystem.py b/examples/levelSet/electroChem/simpleTrenchSystem.py
index 0d78699d73..587b3a6f8c 100644
--- a/examples/levelSet/electroChem/simpleTrenchSystem.py
+++ b/examples/levelSet/electroChem/simpleTrenchSystem.py
@@ -31,7 +31,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Model electrochemical superfill using the CEAC mechanism.
+
 This input file
 is a demonstration of the use of :term:`FiPy`
 for modeling electrodeposition using the CEAC mechanism. The
diff --git a/examples/phase/anisotropy.py b/examples/phase/anisotropy.py
index 4b3f48d446..e457aed975 100755
--- a/examples/phase/anisotropy.py
+++ b/examples/phase/anisotropy.py
@@ -30,7 +30,8 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a dendritic solidification problem.
+
 To convert a liquid material to a solid,  it must be cooled to a 
 temperature below its melting point (known as "undercooling" or "supercooling"). The rate of 
 solidification is often assumed (and experimentally found) to be proportional to the
diff --git a/examples/phase/binaryCoupled.py b/examples/phase/binaryCoupled.py
index f4627b5dd9..8589e051f4 100755
--- a/examples/phase/binaryCoupled.py
+++ b/examples/phase/binaryCoupled.py
@@ -32,7 +32,8 @@
  # ########################################################################
  ##
 
-r"""
+r"""Simultaneously solve a phase-field evolution and solute diffusion problem in one-dimension.
+
 It is straightforward to extend a phase field model to include binary alloys.
 As in :mod:`examples.phase.simple`, we will examine a 1D problem
 
diff --git a/examples/phase/generated/polyxtal.png b/examples/phase/generated/polyxtal.png
new file mode 100644
index 0000000000..8b35a0f541
Binary files /dev/null and b/examples/phase/generated/polyxtal.png differ
diff --git a/examples/phase/impingement/mesh20x20.py b/examples/phase/impingement/mesh20x20.py
index b5302c88f1..a33162fa90 100755
--- a/examples/phase/impingement/mesh20x20.py
+++ b/examples/phase/impingement/mesh20x20.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve for the impingement of four grains in two dimensions.
 
 In the following examples, we solve the same set of equations as in
 :mod:`examples.phase.impingement.mesh40x1`
diff --git a/examples/phase/impingement/mesh40x1.py b/examples/phase/impingement/mesh40x1.py
index d9df63519a..69215a3a65 100755
--- a/examples/phase/impingement/mesh40x1.py
+++ b/examples/phase/impingement/mesh40x1.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve for the impingement of two grains in one dimension.
 
 In this example we solve a coupled phase and orientation equation on a one
 dimensional grid. This is another aspect of the model of Warren, Kobayashi,
diff --git a/examples/phase/index.txt b/examples/phase/index.txt
index 9b7bb3e453..329844a033 100644
--- a/examples/phase/index.txt
+++ b/examples/phase/index.txt
@@ -12,4 +12,6 @@ Phase Field Examples
    examples.phase.anisotropy
    examples.phase.impingement.mesh40x1
    examples.phase.impingement.mesh20x20
+   examples.phase.polyxtal
+   examples.phase.polyxtalCoupled
 
diff --git a/examples/phase/polyxtal.py b/examples/phase/polyxtal.py
new file mode 100755
index 0000000000..f20a6b353c
--- /dev/null
+++ b/examples/phase/polyxtal.py
@@ -0,0 +1,376 @@
+#!/usr/bin/env python
+
+## 
+ # ###################################################################
+ #  FiPy - Python-based finite volume PDE solver
+ # 
+ #  Author: Jonathan Guyer <guyer@nist.gov>
+ #  Author: Daniel Wheeler <daniel.wheeler@nist.gov>
+ #  Author: James Warren   <jwarren@nist.gov>
+ #    mail: NIST
+ #     www: http://www.ctcms.nist.gov/fipy/
+ #  
+ # ========================================================================
+ # This software was developed at the National Institute of Standards
+ # and Technology by employees of the Federal Government in the course
+ # of their official duties.  Pursuant to title 17 Section 105 of the
+ # United States Code this software is not subject to copyright
+ # protection and is in the public domain.  FiPy is an experimental
+ # system.  NIST assumes no responsibility whatsoever for its use by
+ # other parties, and makes no guarantees, expressed or implied, about
+ # its quality, reliability, or any other characteristic.  We would
+ # appreciate acknowledgement if the software is used.
+ # 
+ # This software can be redistributed and/or modified freely
+ # provided that any derivative works bear some notice that they are
+ # derived from it, and any modified versions bear some notice that
+ # they have been modified.
+ # ========================================================================
+ #  
+ # ###################################################################
+ ##
+
+r"""Solve the dendritic growth of nuclei and subsequent grain impingement.
+
+To convert a liquid material to a solid,  it must be cooled to a 
+temperature below its melting point (known as "undercooling" or "supercooling"). The rate of 
+solidification is often assumed (and experimentally found) to be proportional to the
+undercooling. Under the right circumstances, the
+solidification front can become unstable, leading to dendritic
+patterns. 
+Warren, Kobayashi, Lobkovsky and Carter [WarrenPolycrystal]_
+have described a phase field model ("Allen-Cahn", "non-conserved
+Ginsberg-Landau", or "model A" of Hohenberg & Halperin) of such a system,
+including the effects of discrete crystalline orientations (anisotropy).
+
+We start with a regular 2D Cartesian mesh
+
+>>> from fipy import *
+>>> dx = dy = 0.025
+>>> if __name__ == "__main__":
+...     nx = ny = 200
+... else:
+...     nx = ny = 200
+>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
+
+and we'll take fixed timesteps
+
+>>> dt = 5e-4
+
+We consider the simultaneous evolution of a "phase field" variable
+:math:`\phi` (taken to be 0 in the liquid phase and 1 in the solid)
+
+>>> phase = CellVariable(name=r'$\phi$', mesh=mesh, hasOld=True)
+
+a dimensionless undercooling 
+:math:`\Delta T` (:math:`\Delta T = 0` at the melting point)
+
+>>> dT = CellVariable(name=r'$\Delta T$', mesh=mesh, hasOld=True)
+
+and an orientation :math:`-\pi < \theta \le \pi`
+
+>>> theta = ModularVariable(name=r'$\theta$', mesh=mesh, hasOld=True)
+>>> theta.value = -numerix.pi + 0.0001
+
+The ``hasOld`` flag causes the storage of the value of variable from the
+previous timestep. This is necessary for solving equations with
+non-linear coefficients or for coupling between PDEs.
+
+The governing equation for the temperature field is the heat flux
+equation, with a source due to the latent heat of solidification
+
+.. math::
+   
+   \frac{\partial \Delta T}{\partial t} 
+   = D_T \nabla^2 \Delta T 
+   + \frac{\partial \phi}{\partial t}
+   + c\left(T_0 - T\right)
+
+>>> DT = 2.25
+>>> q = Variable(0.)
+>>> T_0 = -0.1
+>>> heatEq = (TransientTerm()
+...           == DiffusionTerm(DT)
+...           + (phase - phase.old) / dt
+...           + q * T_0 - ImplicitSourceTerm(q))
+
+The governing equation for the phase field is
+
+.. math::
+    
+   \tau_{\phi} \frac{\partial \phi}{\partial t} 
+   = \nabla \cdot \mathsf{D} \nabla \phi
+   +   \phi ( 1 - \phi ) m ( \phi , \Delta T)
+   - 2 s \phi | \nabla \theta | - \epsilon^2 \phi | \nabla \theta |^2
+   
+where
+
+.. math::
+
+   m(\phi, \Delta T) 
+   = \phi - \frac{1}{2} 
+   - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 \Delta T \right)
+   
+represents a source of anisotropy. The coefficient 
+:math:`\mathsf{D}` 
+is an anisotropic diffusion tensor in two dimensions
+
+.. math::
+
+   \mathsf{D} = \alpha^2 \left( 1 + c \beta \right)
+   \left[
+   \begin{matrix}
+       1 + c \beta & -c \frac{\partial \beta}{\partial \psi} \\
+       c \frac{\partial \beta}{\partial \psi} & 1 + c \beta
+   \end{matrix}
+   \right]
+   
+where :math:`\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}`,
+:math:`\Phi = \tan \left( \frac{ N } { 2 } \psi \right)`, 
+:math:`\psi = \theta 
++ \arctan \frac{\partial \phi / \partial y}{\partial \phi / \partial x}`, 
+:math:`\theta` is the orientation, and :math:`N` is the symmetry.
+
+>>> alpha = 0.015
+>>> c = 0.02
+>>> N = 4.
+
+>>> psi = theta.arithmeticFaceValue + numerix.arctan2(phase.faceGrad[1], 
+...                                                   phase.faceGrad[0])
+>>> Phi = numerix.tan(N * psi / 2)
+>>> PhiSq = Phi**2
+>>> beta = (1. - PhiSq) / (1. + PhiSq)
+>>> DbetaDpsi = -N * 2 * Phi / (1 + PhiSq)
+>>> Ddia = (1.+ c * beta)
+
+>>> Doff = c * DbetaDpsi
+>>> I0 = Variable(value=((1,0), (0,1)))
+>>> I1 = Variable(value=((0,-1), (1,0)))
+>>> D = alpha**2 * Ddia * (Ddia * I0 + Doff * I1)
+
+With these expressions defined, we can construct the phase field equation
+as
+
+>>> tau_phase = 3e-4
+>>> kappa1 = 0.9
+>>> kappa2 = 20.
+>>> epsilon = 0.008
+>>> s = 0.01
+>>> thetaMag = theta.grad.mag
+>>> phaseEq = (TransientTerm(tau_phase)
+...            == DiffusionTerm(D)
+...            + ImplicitSourceTerm((phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * dT))
+...                                 * (1 - phase)
+...                                 - (2 * s + epsilon**2 * thetaMag) * thetaMag))
+
+The governing equation for orientation is given by
+
+.. math::
+
+   P(\epsilon | \nabla \theta |) \tau_{\theta} \phi^2 
+   \frac{\partial \theta}{\partial t} 
+   = \nabla \cdot \left[ \phi^2 \left( \frac{s}{| \nabla \theta |} 
+   + \epsilon^2 \right) \nabla \theta \right]
+
+where
+
+.. math::
+
+   P(w) = 1 - \exp{(-\beta w)} + \frac{\mu}{\epsilon} \exp{(-\beta w)}
+   
+The ``theta`` equation is built in the following way. The details for
+this equation are fairly involved, see J.A. Warren *et al.*. The main
+detail is that a source must be added to correct for the
+discretization of ``theta`` on the circle.
+
+>>> tau_theta = 3e-3
+>>> mu = 1e3
+>>> gamma = 1e3
+>>> thetaSmallValue = 1e-6
+>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
+>>> beta_theta = 1e5
+>>> expo = epsilon * beta_theta * theta.grad.mag
+>>> expo = (expo < 100.) * (expo - 100.) + 100.
+>>> Pfunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)
+
+>>> gradMagTheta = theta.faceGrad.mag
+>>> eps = 1. / gamma / 10.
+>>> gradMagTheta += (gradMagTheta < eps) * eps
+>>> IGamma = (gradMagTheta > 1. / gamma) * (1 / gradMagTheta - gamma) + gamma
+>>> v_theta = phase.arithmeticFaceValue * (s * IGamma + epsilon**2)
+>>> D_theta = phase.arithmeticFaceValue**2 * (s * IGamma + epsilon**2)
+
+The source term requires the evaluation of the face gradient without
+the modular operator. ``theta``:meth:`~fipy.variables.modularVariable.ModularVariable.getFaceGradNoMod`
+evaluates the gradient without modular arithmetic.
+
+>>> thetaEq = (TransientTerm(tau_theta * phaseMod**2 * Pfunc) 
+...            == DiffusionTerm(D_theta)
+...            + (D_theta * (theta.faceGrad - theta.faceGradNoMod)).divergence)
+
+We seed a circular solidified region in the center
+
+>>> x, y = mesh.cellCenters
+>>> numSeeds = 10
+>>> numerix.random.seed(12345)
+>>> for Cx, Cy, orientation in numerix.random.random([numSeeds, 3]):
+...     radius = dx * 5.
+...     seed = ((x - Cx * nx * dx)**2 + (y - Cy * ny * dy)**2) < radius**2
+...     phase[seed] = 1.
+...     theta[seed] = numerix.pi * (2 * orientation - 1)
+
+and quench the entire simulation domain below the melting point
+
+>>> dT.setValue(-0.5)
+
+In a real solidification process, dendritic branching is induced by small thermal
+fluctuations along an otherwise smooth surface, but the granularity of the
+:class:`~fipy.meshes.mesh.Mesh` is enough "noise" in this case, so we don't need to explicitly
+introduce randomness, the way we did in the Cahn-Hilliard problem.  
+
+FiPy's viewers are utilitarian, striving to let the user see *something*,
+regardless of their operating system or installed packages, so you the default
+color scheme of grain orientation won't be very informative "out of the box".
+Because all of Python is accessible and FiPy is object oriented, it is not hard
+to adapt one of the existing viewers to create a specialized display:
+
+>>> if __name__ == "__main__":
+...     try:
+...         class OrientationViewer(Matplotlib2DGridViewer):
+...             def __init__(self, phase, orientation, title=None, limits={}, **kwlimits):
+...                 self.phase = phase
+...                 Matplotlib2DGridViewer.__init__(self, vars=(orientation,), title=title,
+...                                                 limits=limits, colorbar=None, **kwlimits)
+...                                         
+...                 # make room for non-existent colorbar
+...                 # stolen from matplotlib.colorbar.make_axes
+...                 # https://github.com/matplotlib/matplotlib/blob
+...                 #   /ec1cd2567521c105a451ce15e06de10715f8b54d/lib
+...                 #   /matplotlib/colorbar.py#L838
+...                 fraction = 0.15
+...                 pb = self.axes.get_position(original=True).frozen()
+...                 pad = 0.05
+...                 x1 = 1.0-fraction
+...                 pb1, pbx, pbcb = pb.splitx(x1-pad, x1)
+...                 panchor = (1.0, 0.5)
+...                 self.axes.set_position(pb1)
+...                 self.axes.set_anchor(panchor)
+...         
+...                 # make the gnomon
+...                 fig = self.axes.get_figure()
+...                 self.gnomon = fig.add_axes([0.85, 0.425, 0.15, 0.15], polar=True) 
+...                 self.gnomon.set_thetagrids([180, 270, 0, 90], 
+...                                            [r"$\pm\pi$", r"$-\frac{\pi}{2}$", "$0$", r"$+\frac{\pi}{2}$"],
+...                                            frac=1.3)
+...                 self.gnomon.set_theta_zero_location("N")
+...                 self.gnomon.set_theta_direction(-1)
+...                 self.gnomon.set_rgrids([1.], [""])
+...                 N = 100
+...                 theta = numerix.arange(-numerix.pi, numerix.pi, 2 * numerix.pi / N)
+...                 radii = numerix.ones((N,))
+...                 bars = self.gnomon.bar(theta, radii, width=2 * numerix.pi / N, bottom=0.0)
+...                 colors = self._orientation_and_phase_to_rgb(orientation=numerix.array([theta]), phase=1.)
+...                 for c, t, bar in zip(colors[0], theta, bars):
+...                     bar.set_facecolor(c)
+...                     bar.set_edgecolor(c)
+...              
+...             def _reshape(self, var):
+...                 '''return values of var in an 2D array'''
+...                 return numerix.reshape(numerix.array(var), 
+...                                        var.mesh.shape[::-1])[::-1]
+...     
+...             @staticmethod
+...             def _orientation_and_phase_to_rgb(orientation, phase):
+...                 from matplotlib import colors
+...         
+...                 hsv = numerix.empty(orientation.shape + (3,))
+...                 hsv[..., 0] = (orientation / numerix.pi + 1) / 2.
+...                 hsv[..., 1] = 1.
+...                 hsv[..., 2] = phase
+...                 
+...                 return colors.hsv_to_rgb(hsv)
+...                 
+...             @property
+...             def _data(self):
+...                 '''convert phase and orientation to rgb image array
+...                 
+...                 orientation (-pi, pi) -> hue (0, 1)
+...                 phase (0, 1) -> value (0, 1)
+...                 '''
+...                 orientation = self._reshape(self.vars[0])
+...                 phase = self._reshape(self.phase)
+...                 
+...                 return self._orientation_and_phase_to_rgb(orientation, phase)
+...                                         
+...             def _plot(self):
+...                 self.image.set_data(self._data)
+...                 
+...         from matplotlib import pyplot
+...         pyplot.ion()
+...         w, h = pyplot.figaspect(1.)
+...         fig = pyplot.figure(figsize=(2*w, h))
+...         timer = fig.text(0.1, 0.9, "t = %.3f" % 0, fontsize=18)
+...         
+...         viewer = MultiViewer(viewers=(MatplotlibViewer(vars=dT,
+...                                                        cmap=pyplot.cm.hot,
+...                                                        datamin=-0.5, 
+...                                                        datamax=0.5,
+...                                                        axes=fig.add_subplot(121)),
+...                                       OrientationViewer(phase=phase, 
+...                                                         orientation=theta, 
+...                                                         title=theta.name, 
+...                                                         axes=fig.add_subplot(122))))
+...     except ImportError:
+...         viewer = MultiViewer(viewers=(Viewer(vars=dT,
+...                                              datamin=-0.5, 
+...                                              datamax=0.5),
+...                                       Viewer(vars=phase,
+...                                              datamin=0.,
+...                                              datamax=1.),
+...                                       Viewer(vars=theta,
+...                                              datamin=-numerix.pi,
+...                                              datamax=numerix.pi)))
+>>> viewer.plot()
+
+and iterate the solution in time, plotting as we go,
+
+>>> if __name__ == "__main__":
+...     total_time = 2.
+... else:
+...     total_time = dt * 10
+>>> elapsed = 0.
+>>> save_interval = 0.002
+>>> save_at = save_interval
+
+>>> while elapsed < total_time:
+...     if elapsed > 0.3:
+...         q.value = 100
+...     phase.updateOld()
+...     dT.updateOld()
+...     theta.updateOld()
+...     thetaEq.solve(theta, dt=dt)
+...     phaseEq.solve(phase, dt=dt)
+...     heatEq.solve(dT, dt=dt)
+...     elapsed += dt
+...     if __name__ == "__main__" and elapsed >= save_at:
+...         timer.set_text("t = %.3f" % elapsed)
+...         viewer.plot()
+...         save_at += save_interval
+
+.. image:: polyxtal.*
+   :width: 90%
+   :align: center
+   :alt: undercooling and grain orientation during solidification of a collection of anisotropic seeds
+
+The non-uniform temperature results from the release of latent 
+heat at the solidifying interface. The dendrite arms grow fastest 
+where the temperature gradient is steepest.
+"""
+__docformat__ = 'restructuredtext'
+
+if __name__ == '__main__':
+    import fipy.tests.doctestPlus
+    exec(fipy.tests.doctestPlus._getScript())
+
+    raw_input('finished')
+
diff --git a/examples/phase/polyxtalCoupled.py b/examples/phase/polyxtalCoupled.py
new file mode 100755
index 0000000000..ec86cd34c0
--- /dev/null
+++ b/examples/phase/polyxtalCoupled.py
@@ -0,0 +1,377 @@
+#!/usr/bin/env python
+
+## 
+ # ###################################################################
+ #  FiPy - Python-based finite volume PDE solver
+ # 
+ #  Author: Jonathan Guyer <guyer@nist.gov>
+ #  Author: Daniel Wheeler <daniel.wheeler@nist.gov>
+ #  Author: James Warren   <jwarren@nist.gov>
+ #    mail: NIST
+ #     www: http://www.ctcms.nist.gov/fipy/
+ #  
+ # ========================================================================
+ # This software was developed at the National Institute of Standards
+ # and Technology by employees of the Federal Government in the course
+ # of their official duties.  Pursuant to title 17 Section 105 of the
+ # United States Code this software is not subject to copyright
+ # protection and is in the public domain.  FiPy is an experimental
+ # system.  NIST assumes no responsibility whatsoever for its use by
+ # other parties, and makes no guarantees, expressed or implied, about
+ # its quality, reliability, or any other characteristic.  We would
+ # appreciate acknowledgement if the software is used.
+ # 
+ # This software can be redistributed and/or modified freely
+ # provided that any derivative works bear some notice that they are
+ # derived from it, and any modified versions bear some notice that
+ # they have been modified.
+ # ========================================================================
+ #  
+ # ###################################################################
+ ##
+
+r"""Simultaneously solve the dendritic growth of nuclei and subsequent grain impingement.
+
+To convert a liquid material to a solid,  it must be cooled to a 
+temperature below its melting point (known as "undercooling" or "supercooling"). The rate of 
+solidification is often assumed (and experimentally found) to be proportional to the
+undercooling. Under the right circumstances, the
+solidification front can become unstable, leading to dendritic
+patterns. 
+Warren, Kobayashi, Lobkovsky and Carter [WarrenPolycrystal]_
+have described a phase field model ("Allen-Cahn", "non-conserved
+Ginsberg-Landau", or "model A" of Hohenberg & Halperin) of such a system,
+including the effects of discrete crystalline orientations (anisotropy).
+
+We start with a regular 2D Cartesian mesh
+
+>>> from fipy import *
+>>> dx = dy = 0.025
+>>> if __name__ == "__main__":
+...     nx = ny = 200
+... else:
+...     nx = ny = 200
+>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
+
+and we'll take fixed timesteps
+
+>>> dt = 5e-4
+
+We consider the simultaneous evolution of a "phase field" variable
+:math:`\phi` (taken to be 0 in the liquid phase and 1 in the solid)
+
+>>> phase = CellVariable(name=r'$\phi$', mesh=mesh, hasOld=True)
+
+a dimensionless undercooling 
+:math:`\Delta T` (:math:`\Delta T = 0` at the melting point)
+
+>>> dT = CellVariable(name=r'$\Delta T$', mesh=mesh, hasOld=True)
+
+and an orientation :math:`-\pi < \theta \le \pi`
+
+>>> theta = ModularVariable(name=r'$\theta$', mesh=mesh, hasOld=True)
+>>> theta.value = -numerix.pi + 0.0001
+
+The ``hasOld`` flag causes the storage of the value of variable from the
+previous timestep. This is necessary for solving equations with
+non-linear coefficients or for coupling between PDEs.
+
+The governing equation for the temperature field is the heat flux
+equation, with a source due to the latent heat of solidification
+
+.. math::
+   
+   \frac{\partial \Delta T}{\partial t} 
+   = D_T \nabla^2 \Delta T 
+   + \frac{\partial \phi}{\partial t}
+   + c\left(T_0 - T\right)
+
+>>> DT = 2.25
+>>> q = Variable(0.)
+>>> T_0 = -0.1
+>>> heatEq = (TransientTerm(var=dT)
+...           == DiffusionTerm(coeff=DT, var=dT)
+...           + TransientTerm(var=phase)
+...           + q * T_0 - ImplicitSourceTerm(coeff=q, var=dT))
+
+The governing equation for the phase field is
+
+.. math::
+    
+   \tau_{\phi} \frac{\partial \phi}{\partial t} 
+   = \nabla \cdot \mathsf{D} \nabla \phi
+   +   \phi ( 1 - \phi ) m ( \phi , \Delta T)
+   - 2 s \phi | \nabla \theta | - \epsilon^2 \phi | \nabla \theta |^2
+   
+where
+
+.. math::
+
+   m(\phi, \Delta T) 
+   = \phi - \frac{1}{2} 
+   - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 \Delta T \right)
+   
+represents a source of anisotropy. The coefficient 
+:math:`\mathsf{D}` 
+is an anisotropic diffusion tensor in two dimensions
+
+.. math::
+
+   \mathsf{D} = \alpha^2 \left( 1 + c \beta \right)
+   \left[
+   \begin{matrix}
+       1 + c \beta & -c \frac{\partial \beta}{\partial \psi} \\
+       c \frac{\partial \beta}{\partial \psi} & 1 + c \beta
+   \end{matrix}
+   \right]
+   
+where :math:`\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}`,
+:math:`\Phi = \tan \left( \frac{ N } { 2 } \psi \right)`, 
+:math:`\psi = \theta 
++ \arctan \frac{\partial \phi / \partial y}{\partial \phi / \partial x}`, 
+:math:`\theta` is the orientation, and :math:`N` is the symmetry.
+
+>>> alpha = 0.015
+>>> c = 0.02
+>>> N = 4.
+
+>>> psi = theta.arithmeticFaceValue + numerix.arctan2(phase.faceGrad[1], 
+...                                                   phase.faceGrad[0])
+>>> Phi = numerix.tan(N * psi / 2)
+>>> PhiSq = Phi**2
+>>> beta = (1. - PhiSq) / (1. + PhiSq)
+>>> DbetaDpsi = -N * 2 * Phi / (1 + PhiSq)
+>>> Ddia = (1.+ c * beta)
+
+>>> Doff = c * DbetaDpsi
+>>> I0 = Variable(value=((1,0), (0,1)))
+>>> I1 = Variable(value=((0,-1), (1,0)))
+>>> D = alpha**2 * Ddia * (Ddia * I0 + Doff * I1)
+
+With these expressions defined, we can construct the phase field equation
+as
+
+>>> tau_phase = 3e-4
+>>> kappa1 = 0.9
+>>> kappa2 = 20.
+>>> epsilon = 0.008
+>>> s = 0.01
+>>> thetaMag = theta.grad.mag
+>>> phaseEq = (TransientTerm(coeff=tau_phase, var=phase)
+...            == DiffusionTerm(coeff=D, var=phase)
+...            + ImplicitSourceTerm(coeff=((phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * dT))
+...                                        * (1 - phase)
+...                                        - (2 * s + epsilon**2 * thetaMag) * thetaMag),
+...                                 var=phase))
+
+The governing equation for orientation is given by
+
+.. math::
+
+   P(\epsilon | \nabla \theta |) \tau_{\theta} \phi^2 
+   \frac{\partial \theta}{\partial t} 
+   = \nabla \cdot \left[ \phi^2 \left( \frac{s}{| \nabla \theta |} 
+   + \epsilon^2 \right) \nabla \theta \right]
+
+where
+
+.. math::
+
+   P(w) = 1 - \exp{(-\beta w)} + \frac{\mu}{\epsilon} \exp{(-\beta w)}
+   
+The ``theta`` equation is built in the following way. The details for
+this equation are fairly involved, see J.A. Warren *et al.*. The main
+detail is that a source must be added to correct for the
+discretization of ``theta`` on the circle.
+
+>>> tau_theta = 3e-3
+>>> mu = 1e3
+>>> gamma = 1e3
+>>> thetaSmallValue = 1e-6
+>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
+>>> beta_theta = 1e5
+>>> expo = epsilon * beta_theta * theta.grad.mag
+>>> expo = (expo < 100.) * (expo - 100.) + 100.
+>>> Pfunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)
+
+>>> gradMagTheta = theta.faceGrad.mag
+>>> eps = 1. / gamma / 10.
+>>> gradMagTheta += (gradMagTheta < eps) * eps
+>>> IGamma = (gradMagTheta > 1. / gamma) * (1 / gradMagTheta - gamma) + gamma
+>>> v_theta = phase.arithmeticFaceValue * (s * IGamma + epsilon**2)
+>>> D_theta = phase.arithmeticFaceValue**2 * (s * IGamma + epsilon**2)
+
+The source term requires the evaluation of the face gradient without
+the modular operator. ``theta``:meth:`~fipy.variables.modularVariable.ModularVariable.getFaceGradNoMod`
+evaluates the gradient without modular arithmetic.
+
+>>> thetaEq = (TransientTerm(coeff=tau_theta * phaseMod**2 * Pfunc, var=theta) 
+...            == DiffusionTerm(coeff=D_theta, var=theta)
+...            + PowerLawConvectionTerm(coeff=v_theta * (theta.faceGrad - theta.faceGradNoMod), var=phase))
+
+We seed a circular solidified region in the center
+
+>>> x, y = mesh.cellCenters
+>>> numSeeds = 10
+>>> numerix.random.seed(12345)
+>>> for Cx, Cy, orientation in numerix.random.random([numSeeds, 3]):
+...     radius = dx * 5.
+...     seed = ((x - Cx * nx * dx)**2 + (y - Cy * ny * dy)**2) < radius**2
+...     phase[seed] = 1.
+...     theta[seed] = numerix.pi * (2 * orientation - 1)
+
+and quench the entire simulation domain below the melting point
+
+>>> dT.setValue(-0.5)
+
+In a real solidification process, dendritic branching is induced by small thermal
+fluctuations along an otherwise smooth surface, but the granularity of the
+:class:`~fipy.meshes.mesh.Mesh` is enough "noise" in this case, so we don't need to explicitly
+introduce randomness, the way we did in the Cahn-Hilliard problem.  
+
+FiPy's viewers are utilitarian, striving to let the user see *something*,
+regardless of their operating system or installed packages, so you the default
+color scheme of grain orientation won't be very informative "out of the box".
+Because all of Python is accessible and FiPy is object oriented, it is not hard
+to adapt one of the existing viewers to create a specialized display:
+
+>>> if __name__ == "__main__":
+...     try:
+...         class OrientationViewer(Matplotlib2DGridViewer):
+...             def __init__(self, phase, orientation, title=None, limits={}, **kwlimits):
+...                 self.phase = phase
+...                 Matplotlib2DGridViewer.__init__(self, vars=(orientation,), title=title,
+...                                                 limits=limits, colorbar=None, **kwlimits)
+...                                         
+...                 # make room for non-existent colorbar
+...                 # stolen from matplotlib.colorbar.make_axes
+...                 # https://github.com/matplotlib/matplotlib/blob
+...                 #   /ec1cd2567521c105a451ce15e06de10715f8b54d/lib
+...                 #   /matplotlib/colorbar.py#L838
+...                 fraction = 0.15
+...                 pb = self.axes.get_position(original=True).frozen()
+...                 pad = 0.05
+...                 x1 = 1.0-fraction
+...                 pb1, pbx, pbcb = pb.splitx(x1-pad, x1)
+...                 panchor = (1.0, 0.5)
+...                 self.axes.set_position(pb1)
+...                 self.axes.set_anchor(panchor)
+...         
+...                 # make the gnomon
+...                 fig = self.axes.get_figure()
+...                 self.gnomon = fig.add_axes([0.85, 0.425, 0.15, 0.15], polar=True) 
+...                 self.gnomon.set_thetagrids([180, 270, 0, 90], 
+...                                            [r"$\pm\pi$", r"$-\frac{\pi}{2}$", "$0$", r"$+\frac{\pi}{2}$"],
+...                                            frac=1.3)
+...                 self.gnomon.set_theta_zero_location("N")
+...                 self.gnomon.set_theta_direction(-1)
+...                 self.gnomon.set_rgrids([1.], [""])
+...                 N = 100
+...                 theta = numerix.arange(-numerix.pi, numerix.pi, 2 * numerix.pi / N)
+...                 radii = numerix.ones((N,))
+...                 bars = self.gnomon.bar(theta, radii, width=2 * numerix.pi / N, bottom=0.0)
+...                 colors = self._orientation_and_phase_to_rgb(orientation=numerix.array([theta]), phase=1.)
+...                 for c, t, bar in zip(colors[0], theta, bars):
+...                     bar.set_facecolor(c)
+...                     bar.set_edgecolor(c)
+...              
+...             def _reshape(self, var):
+...                 '''return values of var in an 2D array'''
+...                 return numerix.reshape(numerix.array(var), 
+...                                        var.mesh.shape[::-1])[::-1]
+...     
+...             @staticmethod
+...             def _orientation_and_phase_to_rgb(orientation, phase):
+...                 from matplotlib import colors
+...         
+...                 hsv = numerix.empty(orientation.shape + (3,))
+...                 hsv[..., 0] = (orientation / numerix.pi + 1) / 2.
+...                 hsv[..., 1] = 1.
+...                 hsv[..., 2] = phase
+...                 
+...                 return colors.hsv_to_rgb(hsv)
+...                 
+...             @property
+...             def _data(self):
+...                 '''convert phase and orientation to rgb image array
+...                 
+...                 orientation (-pi, pi) -> hue (0, 1)
+...                 phase (0, 1) -> value (0, 1)
+...                 '''
+...                 orientation = self._reshape(self.vars[0])
+...                 phase = self._reshape(self.phase)
+...                 
+...                 return self._orientation_and_phase_to_rgb(orientation, phase)
+...                                         
+...             def _plot(self):
+...                 self.image.set_data(self._data)
+...                 
+...         from matplotlib import pyplot
+...         pyplot.ion()
+...         w, h = pyplot.figaspect(1.)
+...         fig = pyplot.figure(figsize=(2*w, h))
+...         timer = fig.text(0.1, 0.9, "t = %.3f" % 0, fontsize=18)
+...         
+...         viewer = MultiViewer(viewers=(MatplotlibViewer(vars=dT,
+...                                                        cmap=pyplot.cm.hot,
+...                                                        datamin=-0.5, 
+...                                                        datamax=0.5,
+...                                                        axes=fig.add_subplot(121)),
+...                                       OrientationViewer(phase=phase, 
+...                                                         orientation=theta, 
+...                                                         title=theta.name, 
+...                                                         axes=fig.add_subplot(122))))
+...     except ImportError:
+...         viewer = MultiViewer(viewers=(Viewer(vars=dT,
+...                                              datamin=-0.5, 
+...                                              datamax=0.5),
+...                                       Viewer(vars=phase,
+...                                              datamin=0.,
+...                                              datamax=1.),
+...                                       Viewer(vars=theta,
+...                                              datamin=-numerix.pi,
+...                                              datamax=numerix.pi)))
+>>> viewer.plot()
+
+and iterate the solution in time, plotting as we go,
+
+>>> eq = thetaEq & phaseEq & heatEq
+
+>>> if __name__ == "__main__":
+...     total_time = 2.
+... else:
+...     total_time = dt * 10
+>>> elapsed = 0.
+>>> save_interval = 0.002
+>>> save_at = save_interval
+
+>>> while elapsed < total_time:
+...     if elapsed > 0.3:
+...         q.value = 100
+...     phase.updateOld()
+...     dT.updateOld()
+...     theta.updateOld()
+...     eq.solve(dt=dt)
+...     elapsed += dt
+...     if __name__ == "__main__" and elapsed >= save_at:
+...         timer.set_text("t = %.3f" % elapsed)
+...         viewer.plot()
+...         save_at += save_interval
+
+.. image:: polyxtal.*
+   :width: 90%
+   :align: center
+   :alt: undercooling and grain orientation during solidification of a collection of anisotropic seeds
+
+The non-uniform temperature results from the release of latent 
+heat at the solidifying interface. The dendrite arms grow fastest 
+where the temperature gradient is steepest.
+"""
+__docformat__ = 'restructuredtext'
+
+if __name__ == '__main__':
+    import fipy.tests.doctestPlus
+    exec(fipy.tests.doctestPlus._getScript())
+
+    raw_input('finished')
+
diff --git a/examples/phase/quaternary.py b/examples/phase/quaternary.py
index 3384c7d3bd..2a9142bce5 100755
--- a/examples/phase/quaternary.py
+++ b/examples/phase/quaternary.py
@@ -32,7 +32,8 @@
  # ###################################################################
  ##
 
-r""" 
+r""" Solve a phase-field evolution and diffusion of four species in one-dimension.
+
 The same procedure used to construct the two-component phase field
 diffusion problem in :mod:`examples.phase.binary` can be used to build up a
 system of multiple components. Once again, we'll focus on 1D.
diff --git a/examples/phase/simple.py b/examples/phase/simple.py
index f87aa8f2d8..4e9c37faf1 100755
--- a/examples/phase/simple.py
+++ b/examples/phase/simple.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a phase-field (Allen-Cahn) problem in one-dimension.
 
 To run this example from the base FiPy directory, type
 ``python examples/phase/simple/input.py`` at the command line.  A viewer
diff --git a/examples/phase/test.py b/examples/phase/test.py
index 1fe7678a6e..5c47ba48fe 100755
--- a/examples/phase/test.py
+++ b/examples/phase/test.py
@@ -45,7 +45,9 @@ def _suite():
                                        'quaternary',
                                        'simple',
                                        'symmetry',
-                                       'binaryCoupled'
+                                       'binaryCoupled',
+                                       'polyxtal',
+                                       'polyxtalCoupled'
                                    ), 
                                    base = __name__)
     
diff --git a/examples/reactiveWetting/liquidVapor1D.py b/examples/reactiveWetting/liquidVapor1D.py
index 2dd453c298..f72b0323e1 100755
--- a/examples/reactiveWetting/liquidVapor1D.py
+++ b/examples/reactiveWetting/liquidVapor1D.py
@@ -32,7 +32,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""Solve a single-component, liquid-vapor, van der Waals system.
 
 This example solves a single-component, liquid-vapor, van der Waals system as
 described by Wheeler et *al.* [PhysRevE.82.051601]_. The free energy for this
diff --git a/examples/updating/update0_1to1_0.py b/examples/updating/update0_1to1_0.py
index 87f1d7ad42..34210dd7d9 100755
--- a/examples/updating/update0_1to1_0.py
+++ b/examples/updating/update0_1to1_0.py
@@ -30,7 +30,7 @@
  # ###################################################################
  ##
 
-r"""
+r"""How to update scripts from version 0.1 to 1.0.
 
 It seems unlikely that many users are still running :term:`FiPy` 0.1, but for those
 that are, the syntax of :term:`FiPy` scripts changed considerably between version 0.1
@@ -90,10 +90,10 @@
 >>> from fipy.boundaryConditions.fixedValue import FixedValue
 >>> from fipy.boundaryConditions.fixedFlux import FixedFlux
 >>> boundaryConditions = (
-...     FixedValue(mesh.facesLeft, valueLeft),
-...     FixedValue(mesh.facesRight, valueRight),
-...     FixedFlux(mesh.facesTop, 0.),
-...     FixedFlux(mesh.facesBottom, 0.)
+...     FixedValue(mesh.getFacesLeft(), valueLeft),
+...     FixedValue(mesh.getFacesRight(), valueRight),
+...     FixedFlux(mesh.getFacesTop(), 0.),
+...     FixedFlux(mesh.getFacesBottom(), 0.)
 ...     )
 
 The solution variable is initialized to `valueLeft`:
@@ -161,7 +161,7 @@
 .. index:: numerix
 
 >>> axis = 0
->>> x = mesh.cellCenters[:,axis]
+>>> x = mesh.getCellCenters()[:,axis]
 >>> from fipy.tools import numerix
 >>> CC = 1. - numerix.exp(-convCoeff[axis] * x / diffCoeff)
 >>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
@@ -227,8 +227,8 @@
 >>> valueRight = 1.
 >>> from fipy.boundaryConditions.fixedValue import FixedValue
 >>> boundaryConditions = (
-...     FixedValue(mesh.facesLeft, valueLeft),
-...     FixedValue(mesh.facesRight, valueRight))
+...     FixedValue(mesh.getFacesLeft(), valueLeft),
+...     FixedValue(mesh.getFacesRight(), valueRight))
 
 The creation of the solution variable is unchanged:
     
diff --git a/examples/updating/update1_0to2_0.py b/examples/updating/update1_0to2_0.py
index 66868c3b76..a73a650a72 100644
--- a/examples/updating/update1_0to2_0.py
+++ b/examples/updating/update1_0to2_0.py
@@ -55,21 +55,21 @@
 
  * The dimension axis of a :class:`~fipy.variables.variable.Variable` is now first, not last
    
-   >>> x = mesh.cellCenters[0]
+   >>> x = mesh.getCellCenters()[0]
 
    instead of
    
-   >>> x = mesh.cellCenters[...,0]
+   >>> x = mesh.getCellCenters()[...,0]
        
    This seemingly arbitrary change simplifies a great many things in :term:`FiPy`, but
    the one most noticeable to the user is that you can now write
    
-   >>> x, y = mesh.cellCenters
+   >>> x, y = mesh.getCellCenters()
   
    instead of
    
-   >>> x = mesh.cellCenters[...,0]
-   >>> y = mesh.cellCenters[...,1]
+   >>> x = mesh.getCellCenters()[...,0]
+   >>> y = mesh.getCellCenters()[...,1]
 
    Unfortunately, we cannot reliably automate this conversion, but we find that
    searching for "``...,``" and "``:,``" finds almost everything. Please don't
@@ -110,12 +110,12 @@
       Because vector fields are properly supported, use vector operations to
       manipulate them, such as
       
-      >>> phase.faceGrad.dot((( 0, 1),
+      >>> phase.getFaceGrad().dot((( 0, 1),
       ...                          (-1, 0)))
 
       instead of the hackish
       
-      >>> phase.faceGrad._take((1, 0), axis=1) * (-1, 1)
+      >>> phase.getFaceGrad()._take((1, 0), axis=1) * (-1, 1)
 
  * For internal reasons, :term:`FiPy` now supports 
    :class:`~fipy.variables.cellVariable.CellVariable` and 
@@ -137,13 +137,13 @@
    :class:`~fipy.boundaryConditions.boundaryCondition.BoundaryCondition` now takes a mask, 
    instead of a list of :class:`~fipy.meshes.face.Face` IDs. Now you write
    
-   >>> X, Y = mesh.faceCenters
-   >>> FixedValue(faces=mesh.exteriorFaces & (X**2 < 1e-6), value=...)
+   >>> X, Y = mesh.getFaaceCenters()
+   >>> FixedValue(faces=mesh.getExteriorFaces() & (X**2 < 1e-6), value=...)
        
    instead of
    
-   >>> exteriorFaces = mesh.exteriorFaces
-   >>> X = exteriorFaces.centers[...,0]
+   >>> exteriorFaces = mesh.getExteriorFaces()
+   >>> X = exteriorFaces.getCenters()[...,0]
    >>> FixedValue(faces=exteriorFaces.where(X**2 < 1e-6), value=...)
        
    With the old syntax, a different call to
@@ -152,11 +152,11 @@
    difficult to specify boundary conditions that depended both on position in
    space and on the current values of any other :class:`~fipy.variables.variable.Variable`.
    
-   >>> FixedValue(faces=(mesh.exteriorFaces 
+   >>> FixedValue(faces=(mesh.getExteriorFaces() 
    ...                   & (((X**2 < 1e-6) 
    ...                       & (Y > 3.)) 
-   ...                      | (phi.arithmeticFaceValue 
-   ...                         < sin(gamma.arithmeticFaceValue)))), value=...)
+   ...                      | (phi.getArithmeticFaceValue()
+   ...                         < sin(gamma.getArithmeticFaceValue())))), value=...)
 
    although it probably could have been done with a rather convoluted (and
    slow!) ``filter`` function passed to ``where``. There no longer are any ``filter``
@@ -177,7 +177,7 @@
    >>> for cell in positiveCells:
    ...     initialArray[cell.ID] = 1.
 
-   Although they still exist, we find very lille cause to ever call
+   Although they still exist, we find very little cause to ever call
    :meth:`~fipy.meshes.mesh.Mesh.getCells` 
    or :meth:`fipy.meshes.mesh.Mesh.getFaces`.
    
diff --git a/examples/updating/update2_0to3_0.py b/examples/updating/update2_0to3_0.py
index 5728fcb463..06d75eff66 100644
--- a/examples/updating/update2_0to3_0.py
+++ b/examples/updating/update2_0to3_0.py
@@ -69,7 +69,40 @@
    
    .. note:: the old behavior can be obtained, at least for now, by setting 
              the :envvar:`FIPY_INCLUDE_NUMERIX_ALL` environment variable.
+             
+ * If your equation contains a :class:`~fipy.terms.transientTerm.TransientTerm`,
+   then you must specify the timestep by passing a ``dt=`` argument when calling
+   :meth:`~fipy.terms.term.Term.solve` or :meth:`~fipy.terms.term.Term.sweep`.
    
+The remaining changes are not *required*, but they make scripts easier to read
+and we recommend them. :term:`FiPy` may issue a :exc:`DeprecationWarning` for some cases,
+to indicate that we may not maintain the old syntax indefinitely.
+
+ * "getter" and "setter" methods have been replaced with properties, e.g., use
+ 
+   >>> x, y = mesh.cellCenters
+   
+   instead of
+   
+   >>> x, y = mesh.getCellCenters()
+   
+ * Boundary conditions are better applied with the
+   :meth:`~fipy.variables.cellVariable.CellVariable.constrain` method than with
+   the old :class:`~fipy.boundaryConditions.fixedValue.FixedValue` and
+   :class:`~fipy.boundaryConditions.fixedFlux.FixedFlux` classes. See
+   :ref:`BoundaryConditions`.
+  
+ * Individual :class:`~fipy.meshes.mesh.Mesh` classes should be imported
+   directly from :mod:`fipy.meshes` and not :mod:`fipy.meshes.numMesh`.
+   
+ * The :term:`Gmsh` meshes now have simplified names:
+   :class:`~fipy.meshes.gmshMesh.Gmsh2D` instead of
+   :class:`~fipy.meshes.gmshMesh.GmshImporter2D`,
+   :class:`~fipy.meshes.gmshMesh.Gmsh3D` instead of
+   :class:`~fipy.meshes.gmshMesh.GmshImporter3D`, and
+   :class:`~fipy.meshes.gmshMesh.Gmsh2DIn3DSpace` instead of
+   :class:`~fipy.meshes.gmshMesh.GmshImporter2DIn3DSpace`.
+
 .. _mailing list:         http://www.ctcms.nist.gov/fipy/mail.html
 """
 __docformat__ = 'restructuredtext'