diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0f0f630b023b09f998fb9fe2902939cb9a32ba9f..25e262d074551a616c5d65ad5874871afb936137 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -52,6 +52,21 @@ jobs:
                 build_py_project_in_conda_env
                 run_examples
 
+    docs:
+        name: Documentation
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v2
+        -   name: "Main Script"
+            run: |
+                sudo apt-get update
+                sudo apt-get install openmpi-bin libopenmpi-dev
+                CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+                curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/ci-support.sh
+                . ci-support.sh
+                build_py_project_in_conda_env
+                build_docs
+
 # vim: sw=4
 
 
diff --git a/.gitignore b/.gitignore
index 94648fab1309424c88e2836ab4034eab05049871..14be95d194bd112b5596328cf5bafcfc092e0cce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,3 +33,5 @@ run-debug-*
 
 .cache
 .pytest_cache
+
+*.json
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 66fc4e491c73a020bba5baee77f85048805696be..ebe4262476c11ae643380f4111a45f0d083f304c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -128,8 +128,6 @@ Documentation:
   - ". ./build-docs.sh"
   tags:
   - python3
-  only:
-  - master
 
 Flake8:
   script:
diff --git a/doc/conf.py b/doc/conf.py
index d75d2aa2ac850052f2c0ad196affaea5b0ec9e98..dc8f0e244b6a4212b817e4751f54e348bbe9661b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -170,7 +170,7 @@ html_sidebars = {
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+#html_static_path = ['_static']
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
@@ -327,4 +327,4 @@ intersphinx_mapping = {
     'https://documen.tician.de/meshmode/': None,
     'http://documen.tician.de/loopy/': None,
     }
-autoclass_content = "both"
+autoclass_content = "class"
diff --git a/doc/discretization.rst b/doc/discretization.rst
new file mode 100644
index 0000000000000000000000000000000000000000..414111348b5de7a8261b088bdf5b5a2c28680a24
--- /dev/null
+++ b/doc/discretization.rst
@@ -0,0 +1,12 @@
+Discretization
+==============
+
+.. module:: grudge
+
+.. automodule:: grudge.discretization
+
+Discretization with Eager Evaluation
+====================================
+
+.. automodule:: grudge.eager
+
diff --git a/doc/eager.rst b/doc/eager.rst
deleted file mode 100644
index 5d2bdd5b071a15623134a1f549589ce6d0178f13..0000000000000000000000000000000000000000
--- a/doc/eager.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Eagerly-evaluated Operation
-===========================
-
-.. automodule:: grudge.eager
-
diff --git a/doc/index.rst b/doc/index.rst
index 6d7782bc950b12431961f3bd0b8ee30deba4d050..7b6ce79d0a108cdcc8df3d05c00af34ee7213190 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -6,9 +6,9 @@ Contents:
 .. toctree::
     :maxdepth: 2
 
-    misc
+    discretization
     symbolic
-    eager
+    misc
 
 
 Indices and Tables
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 7c0573c2e5ec3515843b6f8c377ad10301256aaa..a5452a24b9e690ec14cd0b5bf455f711f3b7f2e7 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -30,12 +30,12 @@ import numpy as np  # noqa: F401
 from meshmode.array_context import ArrayContext
 
 
-# FIXME Naming not ideal
-class DiscretizationBase(object):
-    pass
+__doc__ = """
+.. autoclass:: DGDiscretizationWithBoundaries
+"""
 
 
-class DGDiscretizationWithBoundaries(DiscretizationBase):
+class DGDiscretizationWithBoundaries(object):
     """
     .. automethod :: discr_from_dd
     .. automethod :: connection_from_dds
@@ -53,10 +53,11 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         """
         :param quad_tag_to_group_factory: A mapping from quadrature tags (typically
             strings--but may be any hashable/comparable object) to a
-            :class:`meshmode.discretization.ElementGroupFactory` indicating with
-            which quadrature discretization the operations are to be carried out,
-            or *None* to indicate that operations with this quadrature tag should
-            be carried out with the standard volume discretization.
+            :class:`~meshmode.discretization.poly_element.ElementGroupFactory`
+            indicating with which quadrature discretization the operations are
+            to be carried out, or *None* to indicate that operations with this
+            quadrature tag should be carried out with the standard volume
+            discretization.
         """
 
         self._setup_actx = array_context
@@ -390,7 +391,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     def order(self):
         from warnings import warn
         warn("DGDiscretizationWithBoundaries.order is deprecated, "
-                "consider the orders of element groups instead. "
+                "consider using the orders of element groups instead. "
                 "'order' will go away in 2021.",
                 DeprecationWarning, stacklevel=2)
 
diff --git a/grudge/eager.py b/grudge/eager.py
index 87368161e5045b2a7c1b508b89cf60ebb336053a..4cff62c5efe9dda454cff20b06789f0db0412f44 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -45,15 +45,28 @@ __doc__ = """
 
 class EagerDGDiscretization(DGDiscretizationWithBoundaries):
     """
+    Inherits from :class:`~grudge.discretization.DGDiscretizationWithBoundaries`.
+
+    .. automethod:: __init__
     .. automethod:: project
     .. automethod:: nodes
+
     .. automethod:: grad
+    .. automethod:: d_dx
     .. automethod:: div
+
     .. automethod:: weak_grad
+    .. automethod:: weak_d_dx
     .. automethod:: weak_div
+
     .. automethod:: normal
     .. automethod:: inverse_mass
     .. automethod:: face_mass
+
+    .. automethod:: norm
+    .. automethod:: nodal_sum
+    .. automethod:: nodal_min
+    .. automethod:: nodal_max
     """
 
     def interp(self, src, tgt, vec):
@@ -64,6 +77,14 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         return self.project(src, tgt, vec)
 
     def project(self, src, tgt, vec):
+        """Project from one discretization to another, e.g. from the
+        volume to the boundary, or from the base to the an overintegrated
+        quadrature discretization.
+
+        :arg src: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one
+        :arg tgt: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one
+        :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+        """
         if (isinstance(vec, np.ndarray)
                 and vec.dtype.char == "O"
                 and not isinstance(vec, DOFArray)):
@@ -72,19 +93,80 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
 
         return self.connection_from_dds(src, tgt)(vec)
 
-    def nodes(self):
-        return self._volume_discr.nodes()
+    def nodes(self, dd=None):
+        r"""Return the nodes of a discretization.
+
+        :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+            Defaults to the base volume discretization.
+        :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+        """
+        if dd is None:
+            return self._volume_discr.nodes()
+        else:
+            return self.discr_from_dd(dd).nodes()
+
+    # {{{ derivatives
 
     @memoize_method
     def _bound_grad(self):
         return bind(self, sym.nabla(self.dim) * sym.Variable("u"), local_only=True)
 
     def grad(self, vec):
+        r"""Return the gradient of the volume function represented by *vec*.
+
+        :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+        :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+        """
         return self._bound_grad()(u=vec)
 
+    @memoize_method
+    def _bound_d_dx(self, xyz_axis):
+        return bind(self, sym.nabla(self.dim)[xyz_axis] * sym.Variable("u"),
+                local_only=True)
+
+    def d_dx(self, xyz_axis, vec):
+        r"""Return the derivative along axis *xyz_axis* of the volume function
+        represented by *vec*.
+
+        :arg xyz_axis: an integer indicating the axis along which the derivative
+            is taken
+        :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+        :returns: a :class:`~meshmode.dof_array.DOFArray`\ s
+        """
+        return self._bound_d_dx(xyz_axis)(u=vec)
+
+    def _div_helper(self, diff_func, vecs):
+        if not isinstance(vecs, np.ndarray):
+            raise TypeError("argument must be an object array")
+        assert vecs.dtype == np.object
+
+        if vecs.shape[-1] != self.ambient_dim:
+            raise ValueError("last dimension of *vecs* argument must match "
+                    "ambient dimension")
+
+        if len(vecs.shape) == 1:
+            return sum(diff_func(i, vec_i) for i, vec_i in enumerate(vecs))
+        else:
+            result = np.zeros(vecs.shape[:-1], dtype=object)
+            for idx in np.ndindex(vecs.shape[:-1]):
+                result[idx] = sum(
+                        diff_func(i, vec_i) for i, vec_i in enumerate(vecs[idx]))
+            return result
+
     def div(self, vecs):
-        return sum(
-                self.grad(vec_i)[i] for i, vec_i in enumerate(vecs))
+        r"""Return the divergence of the vector volume function
+        represented by *vecs*.
+
+        :arg vec: an object array of
+            a :class:`~meshmode.dof_array.DOFArray`\ s,
+            where the last axis of the array must have length
+            matching the volume dimension.
+        :returns: a :class:`~meshmode.dof_array.DOFArray`
+        """
+
+        return self._div_helper(
+                lambda i, subvec: self.d_dx(i, subvec),
+                vecs)
 
     @memoize_method
     def _bound_weak_grad(self, dd):
@@ -93,6 +175,16 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
                 local_only=True)
 
     def weak_grad(self, *args):
+        r"""Return the "weak gradient" of the volume function represented by
+        *vec*.
+
+        May be called with ``(vecs)`` or ``(dd, vecs)``.
+
+        :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+            Defaults to the base volume discretization if not provided.
+        :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+        :returns: an object array of :class:`~meshmode.dof_array.DOFArray`\ s
+        """
         if len(args) == 1:
             vec, = args
             dd = sym.DOFDesc("vol", sym.QTAG_NONE)
@@ -103,7 +195,48 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
 
         return self._bound_weak_grad(dd)(u=vec)
 
+    @memoize_method
+    def _bound_weak_d_dx(self, dd, xyz_axis):
+        return bind(self,
+                sym.stiffness_t(self.dim, dd_in=dd)[xyz_axis]
+                * sym.Variable("u", dd),
+                local_only=True)
+
+    def weak_d_dx(self, *args):
+        r"""Return the derivative along axis *xyz_axis* of the volume function
+        represented by *vec*.
+
+        May be called with ``(xyz_axis, vecs)`` or ``(dd, xyz_axis, vecs)``.
+
+        :arg xyz_axis: an integer indicating the axis along which the derivative
+            is taken
+        :arg vec: a :class:`~meshmode.dof_array.DOFArray`
+        :returns: a :class:`~meshmode.dof_array.DOFArray`\ s
+        """
+        if len(args) == 2:
+            xyz_axis, vec = args
+            dd = sym.DOFDesc("vol", sym.QTAG_NONE)
+        elif len(args) == 3:
+            dd, xyz_axis, vec = args
+        else:
+            raise TypeError("invalid number of arguments")
+
+        return self._bound_weak_d_dx(dd, xyz_axis)(u=vec)
+
     def weak_div(self, *args):
+        r"""Return the "weak divergence" of the vector volume function
+        represented by *vecs*.
+
+        May be called with ``(vecs)`` or ``(dd, vecs)``.
+
+        :arg dd: a :class:`~grudge.sym.DOFDesc`, or a value convertible to one.
+            Defaults to the base volume discretization if not provided.
+        :arg vec: a object array of
+            a :class:`~meshmode.dof_array.DOFArray`\ s,
+            where the last axis of the array must have length
+            matching the volume dimension.
+        :returns: a :class:`~meshmode.dof_array.DOFArray`
+        """
         if len(args) == 1:
             vecs, = args
             dd = sym.DOFDesc("vol", sym.QTAG_NONE)
@@ -112,8 +245,11 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         else:
             raise TypeError("invalid number of arguments")
 
-        return sum(
-                self.weak_grad(dd, vec_i)[i] for i, vec_i in enumerate(vecs))
+        return self._div_helper(
+                lambda i, subvec: self.weak_d_dx(dd, i, subvec),
+                vecs)
+
+    # }}}
 
     @memoize_method
     def normal(self, dd):
@@ -221,6 +357,8 @@ def interior_trace_pair(discrwb, vec):
     return TracePair("int_faces", i, e)
 
 
+# {{{ distributed-memory functionality
+
 class _RankBoundaryCommunication:
     base_tag = 1273
 
@@ -295,5 +433,7 @@ def cross_rank_trace_pairs(discrwb, vec, tag=None):
     else:
         return _cross_rank_trace_pairs_scalar_field(discrwb, vec, tag=tag)
 
+# }}}
+
 
 # vim: foldmethod=marker
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 4ecd920c8d82ee429855183b319fd191ed524a9f..df886d4571379570305e46087836d16c88ba28c3 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -94,12 +94,12 @@ class Operator(pymbolic.primitives.Expression):
     """
     .. attribute:: dd_in
 
-        an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
+        an instance of :class:`~grudge.sym.DOFDesc` describing the
         input discretization.
 
     .. attribute:: dd_out
 
-        an instance of :class:`~grudge.symbolic.primitives.DOFDesc` describing the
+        an instance of :class:`~grudge.sym.DOFDesc` describing the
         output discretization.
     """
 
@@ -477,7 +477,10 @@ class VandermondeOperator(ElementwiseLinearOperator):
 
 class MassOperatorBase(Operator):
     """
-    :attr:`dd_in` and :attr:`dd_out` may be surface or volume discretizations.
+    Inherits from :class:`Operator`.
+
+    :attr:`~Operator.dd_in` and :attr:`~Operator.dd_out` may be surface or volume
+    discretizations.
     """
 
     def __init__(self, dd_in=None, dd_out=None):
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 78261d658e2871379f042ee0d714e2d8cd9e3547..0a49034568f6ec45c25b45018e9451fb866fe361 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -55,6 +55,8 @@ class ExpressionBase(prim.Expression):
 
 __doc__ = """
 
+.. currentmodule:: grudge.sym
+
 DOF description
 ^^^^^^^^^^^^^^^