diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index f6975db2c4bd4f82494fd0c11199bb615fcc5f46..31b1d5061757709ab0257a93281bd987699b0796 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -359,7 +359,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
         with cl.CommandQueue(self.cl_context) as queue:
             return bind(self, sym.weights_and_area_elements(
                 self.ambient_dim,
-                where=sym.QBX_SOURCE_QUAD_STAGE2))(queue).with_queue(None)
+                dofdesc=sym.QBX_SOURCE_QUAD_STAGE2))(queue).with_queue(None)
 
     # }}}
 
diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py
index 83880aefe5e59bfe1279e3c7981c8881d13e143b..829e706c481227bbc49b3e524094555cb1f541b8 100644
--- a/pytential/qbx/refinement.py
+++ b/pytential/qbx/refinement.py
@@ -365,7 +365,7 @@ class RefinerWrangler(TreeWranglerBase):
         source_danger_zone_radii_by_panel = bind(lpot_source,
                 sym._source_danger_zone_radii(
                     lpot_source.ambient_dim,
-                    where=sym.GRANULARITY_ELEMENT))(self.queue)
+                    dofdesc=sym.GRANULARITY_ELEMENT))(self.queue)
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
         evt = knl(
@@ -607,7 +607,7 @@ def refine_for_global_qbx(lpot_source, wrangler,
                 from pytential import bind, sym
                 quad_resolution = bind(lpot_source, sym._quad_resolution(
                     lpot_source.ambient_dim,
-                    where=sym.GRANULARITY_ELEMENT))(wrangler.queue)
+                    dofdesc=sym.GRANULARITY_ELEMENT))(wrangler.queue)
 
                 violates_kernel_length_scale = \
                         wrangler.check_element_prop_threshold(
@@ -627,7 +627,7 @@ def refine_for_global_qbx(lpot_source, wrangler,
                 scaled_max_curv = bind(lpot_source,
                     sym.ElementwiseMax(
                         sym._scaled_max_curvature(lpot_source.ambient_dim),
-                        where=sym.GRANULARITY_ELEMENT))(wrangler.queue)
+                        dofdesc=sym.GRANULARITY_ELEMENT))(wrangler.queue)
 
                 violates_scaled_max_curv = \
                         wrangler.check_element_prop_threshold(
diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py
index f1f9fba13332df07925adb001c929e457a5da7cd..af8f80315cc2ac38ade83bd115de200628b89d28 100644
--- a/pytential/symbolic/execution.py
+++ b/pytential/symbolic/execution.py
@@ -144,11 +144,11 @@ class EvaluationMapper(EvaluationMapperBase):
 
             return result
 
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
         operand = self.rec(expr.operand)
         assert operand.shape == (discr.nnodes,)
 
-        where = sym.as_dofdesc(expr.where)
+        where = sym.as_dofdesc(expr.dofdesc)
         if where.granularity == sym.GRANULARITY_NODE:
             return _reduce(discr.nnodes,
                     node_knl(),
@@ -170,7 +170,7 @@ class EvaluationMapper(EvaluationMapperBase):
         return self._map_elementwise_reduction("max", expr)
 
     def map_ones(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
 
         result = (discr
                 .empty(queue=self.queue, dtype=discr.real_dtype)
@@ -180,12 +180,12 @@ class EvaluationMapper(EvaluationMapperBase):
         return result
 
     def map_node_coordinate_component(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
         return discr.nodes()[expr.ambient_axis] \
                 .with_queue(self.queue)
 
     def map_num_reference_derivative(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
 
         from pytools import flatten
         ref_axes = flatten([axis] * mult for axis, mult in expr.ref_axes)
@@ -195,7 +195,7 @@ class EvaluationMapper(EvaluationMapperBase):
                         .with_queue(self.queue)
 
     def map_q_weight(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
         return discr.quad_weights(self.queue) \
                 .with_queue(self.queue)
 
@@ -207,11 +207,11 @@ class EvaluationMapper(EvaluationMapperBase):
         except KeyError:
             bound_op = bind(
                     expr.expression,
-                    self.places[expr.where],
+                    self.places[expr.dofdesc],
                     self.bound_expr.iprec)
             bound_op_cache[expr] = bound_op
 
-        scipy_op = bound_op.scipy_op(expr.variable_name, expr.where,
+        scipy_op = bound_op.scipy_op(expr.variable_name, expr.dofdesc,
                 **dict((var_name, self.rec(var_expr))
                     for var_name, var_expr in six.iteritems(expr.extra_vars)))
 
diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py
index f7f004806890cf2c098e4fabe9fd706b1e92bf32..da2d449cf57d263e1fdf17e116fbc3ca51df5ff8 100644
--- a/pytential/symbolic/mappers.py
+++ b/pytential/symbolic/mappers.py
@@ -62,14 +62,14 @@ class IdentityMapper(IdentityMapperBase):
     map_node_max = map_node_sum
 
     def map_elementwise_sum(self, expr):
-        return type(expr)(self.rec(expr.operand), expr.where)
+        return type(expr)(self.rec(expr.operand), expr.dofdesc)
 
     map_elementwise_min = map_elementwise_sum
     map_elementwise_max = map_elementwise_sum
 
     def map_num_reference_derivative(self, expr):
         return type(expr)(expr.ref_axes, self.rec(expr.operand),
-                expr.where)
+                expr.dofdesc)
 
     # {{{ childless -- no need to rebuild
 
@@ -92,7 +92,7 @@ class IdentityMapper(IdentityMapperBase):
                 dict([
                     (name, self.rec(name_expr))
                     for name, name_expr in six.iteritems(expr.extra_vars)]),
-                expr.where)
+                expr.dofdesc)
 
     def map_int_g(self, expr):
         return expr.copy(
@@ -177,7 +177,7 @@ class EvaluationMapper(EvaluationMapperBase):
     def map_num_reference_derivative(self, expr):
         return componentwise(
                 lambda subexpr: type(expr)(
-                        expr.ref_axes, self.rec(subexpr), expr.where),
+                        expr.ref_axes, self.rec(subexpr), expr.dofdesc),
                 expr.operand)
 
     def map_int_g(self, expr):
@@ -224,7 +224,7 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
         return dofdesc
 
     def map_ones(self, expr):
-        return type(expr)(where=self._as_dofdesc(expr.where))
+        return type(expr)(dofdesc=self._as_dofdesc(expr.dofdesc))
 
     map_q_weight = map_ones
 
@@ -232,23 +232,23 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
         return type(expr)(
                 expr.ambient_axis,
                 expr.ref_axis,
-                self._as_dofdesc(expr.where))
+                self._as_dofdesc(expr.dofdesc))
 
     def map_node_coordinate_component(self, expr):
         return type(expr)(
                 expr.ambient_axis,
-                self._as_dofdesc(expr.where))
+                self._as_dofdesc(expr.dofdesc))
 
     def map_num_reference_derivative(self, expr):
         return type(expr)(
                 expr.ref_axes,
                 self.rec(expr.operand),
-                self._as_dofdesc(expr.where))
+                self._as_dofdesc(expr.dofdesc))
 
     def map_elementwise_sum(self, expr):
         return type(expr)(
                 self.rec(expr.operand),
-                self._as_dofdesc(expr.where))
+                self._as_dofdesc(expr.dofdesc))
 
     map_elementwise_min = map_elementwise_sum
     map_elementwise_max = map_elementwise_sum
@@ -272,8 +272,7 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
                     ))
 
     def map_inverse(self, expr):
-        dofdesc = prim.as_dofdesc(expr.where)
-
+        dofdesc = prim.as_dofdesc(expr.dofdesc)
         if dofdesc.geometry is None:
             dofdesc = dofdesc.copy(geometry=self.default_where)
 
@@ -585,7 +584,7 @@ def stringify_where(where):
 class StringifyMapper(BaseStringifyMapper):
 
     def map_ones(self, expr, enclosing_prec):
-        return "Ones[%s]" % stringify_where(expr.where)
+        return "Ones[%s]" % stringify_where(expr.dofdesc)
 
     def map_inverse(self, expr, enclosing_prec):
         return "Solve(%s = %s {%s})" % (
@@ -603,17 +602,17 @@ class StringifyMapper(BaseStringifyMapper):
 
     def map_elementwise_sum(self, expr, enclosing_prec):
         return "ElwiseSum[%s](%s)" % (
-                stringify_where(expr.where),
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_elementwise_min(self, expr, enclosing_prec):
         return "ElwiseMin[%s](%s)" % (
-                stringify_where(expr.where),
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_elementwise_max(self, expr, enclosing_prec):
         return "ElwiseMax[%s](%s)" % (
-                stringify_where(expr.where),
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_node_max(self, expr, enclosing_prec):
@@ -624,7 +623,7 @@ class StringifyMapper(BaseStringifyMapper):
 
     def map_node_coordinate_component(self, expr, enclosing_prec):
         return "x%d[%s]" % (expr.ambient_axis,
-                stringify_where(expr.where))
+                stringify_where(expr.dofdesc))
 
     def map_num_reference_derivative(self, expr, enclosing_prec):
         diff_op = " ".join(
@@ -635,7 +634,7 @@ class StringifyMapper(BaseStringifyMapper):
 
         result = "%s[%s] %s" % (
                 diff_op,
-                stringify_where(expr.where),
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_PRODUCT),
                 )
 
@@ -645,10 +644,10 @@ class StringifyMapper(BaseStringifyMapper):
             return result
 
     def map_parametrization_derivative(self, expr, enclosing_prec):
-        return "dx/dr[%s]" % (stringify_where(expr.where))
+        return "dx/dr[%s]" % (stringify_where(expr.dofdesc))
 
     def map_q_weight(self, expr, enclosing_prec):
-        return "w_quad[%s]" % stringify_where(expr.where)
+        return "w_quad[%s]" % stringify_where(expr.dofdesc)
 
     def _stringify_kernel_args(self, kernel_arguments):
         if not kernel_arguments:
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index a73e0c4d27f75c2bbf19bf2ae93439309f4cde7b..d163e39fcb98f0dbf7e857f6e466f4f90a5a95df 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -288,11 +288,11 @@ class MatrixBuilderBase(EvaluationMapperBase):
         op = sym.NumReferenceDerivative(
                 ref_axes=expr.ref_axes,
                 operand=sym.var("u"),
-                where=expr.where)
+                dofdesc=expr.dofdesc)
         return bind(self.places, op)(self.queue, u=rec_operand).get()
 
     def map_node_coordinate_component(self, expr):
-        op = sym.NodeCoordinateComponent(expr.ambient_axis, where=expr.where)
+        op = sym.NodeCoordinateComponent(expr.ambient_axis, dofdesc=expr.dofdesc)
         return bind(self.places, op)(self.queue).get()
 
     def map_call(self, expr):
diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py
index a8dc0c42390e7d56dc4220a7f13fdfe43e440264..d89d393c4948a2944d4efa4da2b89db65e55d4b1 100644
--- a/pytential/symbolic/pde/maxwell/__init__.py
+++ b/pytential/symbolic/pde/maxwell/__init__.py
@@ -70,7 +70,8 @@ def get_sym_maxwell_point_source(kernel, jxyz, k):
 
 # {{{ plane wave
 
-def get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, where=None):
+def get_sym_maxwell_plane_wave(amplitude_vec, v, omega,
+        epsilon=1, mu=1, dofdesc=None):
     r"""Return a symbolic expression that, when bound to a
     :class:`pytential.source.PointPotentialSource` will yield
     a field satisfying Maxwell's equations.
@@ -96,7 +97,7 @@ def get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, where=N
 
     # NOTE: for complex, need to ensure real(n).dot(imag(n)) = 0  (7.15)
 
-    x = sym.nodes(3, where).as_vector()
+    x = sym.nodes(3, dofdesc=dofdesc).as_vector()
 
     v_mag_squared = sym.cse(np.dot(v, v), "v_mag_squared")
     n = v/sym.sqrt(v_mag_squared)
diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py
index 840d804429344f16c62e8f388c654c65dffc6f5b..882b1b3c53cb53830bb94796b887a71494fdae32 100644
--- a/pytential/symbolic/pde/scalar.py
+++ b/pytential/symbolic/pde/scalar.py
@@ -48,15 +48,16 @@ class L2WeightedPDEOperator(object):
             warn("should use L2 weighting in %s" % type(self).__name__,
                     stacklevel=3)
 
-    def get_weight(self, where=None):
+    def get_weight(self, dofdesc=None):
         if self.use_l2_weighting:
-            return cse(area_element(self.kernel.dim, where=where)*QWeight(where))
+            return cse(area_element(self.kernel.dim, dofdesc=dofdesc)
+                    * QWeight(dofdesc=dofdesc))
         else:
             return 1
 
-    def get_sqrt_weight(self, where=None):
+    def get_sqrt_weight(self, dofdesc=None):
         if self.use_l2_weighting:
-            return sqrt_jac_q_weight(self.kernel.dim, where=where)
+            return sqrt_jac_q_weight(self.kernel.dim, dofdesc=dofdesc)
         else:
             return 1
 
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index aa730beb47cc6b320366ce087177f57e580ade4a..827be501b53fe853f2564f2d7cd081c177f1d936 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 import six
 from six.moves import intern
 from warnings import warn
+from functools import wraps
 
 import numpy as np
 from pymbolic.primitives import (  # noqa: F401,N813
@@ -41,9 +42,10 @@ from functools import partial
 
 
 __doc__ = """
-.. |where-blurb| replace:: A symbolic name for a geometric object (such
-    as a :class:`~meshmode.discretization.Discretization`) or a
-    :class:`DOFDescriptor`.
+.. |dofdesc-blurb| replace:: A
+    :class:`~pytential.symbolic.primitives.DOFDescriptor` or a symbolic
+    name for a geometric object (such as a
+    :class:`~meshmode.discretization.Discretization`).
 
 
 Object types
@@ -74,6 +76,23 @@ objects occur as part of a symbolic operator representation:
 those), which is not visible as an array axis in symbolic code. (They're
 visible only once evaluated.)
 
+DOF Description
+^^^^^^^^^^^^^^^
+
+.. autoclass:: DEFAULT_SOURCE
+.. autoclass:: DEFAULT_TARGET
+
+.. autoclass:: QBX_SOURCE_STAGE1
+.. autoclass:: QBX_SOURCE_STAGE2
+.. autoclass:: QBX_SOURCE_QUAD_STAGE2
+
+.. autoclass:: GRANULARITY_NODE
+.. autoclass:: GRANULARITY_CENTER
+.. autoclass:: GRANULARITY_ELEMENT
+
+.. autoclass:: DOFDescriptor
+.. autofunction:: as_dofdesc
+
 Placeholders
 ^^^^^^^^^^^^
 
@@ -201,6 +220,23 @@ Pretty-printing expressions
 """
 
 
+def _deprecate_kwargs(oldkey, newkey):
+    def super_wrapper(func):
+        @wraps(func)
+        def wrapper(*args, **kwargs):
+            if oldkey in kwargs:
+                warn("using `{}` keyword is deprecated. "
+                        "use `{}` instead".format(oldkey, newkey),
+                        DeprecationWarning, stacklevel=2)
+                kwargs[newkey] = kwargs[oldkey]
+                del kwargs[oldkey]
+
+            return func(*args, **kwargs)
+        return wrapper
+
+    return super_wrapper
+
+
 # {{{ dof descriptors
 
 class DEFAULT_SOURCE:  # noqa: N801
@@ -405,8 +441,9 @@ def make_sym_mv(name, num_components):
     return MultiVector(make_sym_vector(name, num_components))
 
 
-def make_sym_surface_mv(name, ambient_dim, dim, where=None):
-    par_grad = parametrization_derivative_matrix(ambient_dim, dim, where)
+@_deprecate_kwargs('where', 'dofdesc')
+def make_sym_surface_mv(name, ambient_dim, dim, dofdesc=None):
+    par_grad = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     return sum(
             var("%s%d" % (name, i))
@@ -471,17 +508,24 @@ class DiscretizationProperty(Expression):
     further arguments.
     """
 
-    init_arg_names = ("where",)
+    init_arg_names = ("dofdesc",)
 
-    def __init__(self, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, dofdesc=None):
         """
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
 
-        self.where = as_dofdesc(where)
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.where,)
+        return (self.dofdesc,)
 
 
 # {{{ discretization properties
@@ -494,29 +538,31 @@ class QWeight(DiscretizationProperty):
 
 class NodeCoordinateComponent(DiscretizationProperty):
 
-    init_arg_names = ("ambient_axis", "where")
+    init_arg_names = ("ambient_axis", "dofdesc")
 
-    def __init__(self, ambient_axis, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, ambient_axis, dofdesc=None):
         """
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
         self.ambient_axis = ambient_axis
-        DiscretizationProperty.__init__(self, where)
+        DiscretizationProperty.__init__(self, dofdesc)
 
     def __getinitargs__(self):
-        return (self.ambient_axis, self.where)
+        return (self.ambient_axis, self.dofdesc)
 
     mapper_method = intern("map_node_coordinate_component")
 
 
-def nodes(ambient_dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def nodes(ambient_dim, dofdesc=None):
     """Return a :class:`pymbolic.geometric_algebra.MultiVector` of node
     locations.
     """
 
     return MultiVector(
             make_obj_array([
-                NodeCoordinateComponent(i, where)
+                NodeCoordinateComponent(i, dofdesc)
                 for i in range(ambient_dim)]))
 
 
@@ -526,22 +572,24 @@ class NumReferenceDerivative(DiscretizationProperty):
     reference coordinates.
     """
 
-    init_arg_names = ("ref_axes", "operand", "where")
+    init_arg_names = ("ref_axes", "operand", "dofdesc")
 
-    def __new__(cls, ref_axes=None, operand=None, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __new__(cls, ref_axes=None, operand=None, dofdesc=None):
         # If the constructor is handed a multivector object, return an
         # object array of the operator applied to each of the
         # coefficients in the multivector.
 
         if isinstance(operand, np.ndarray):
             def make_op(operand_i):
-                return cls(ref_axes, operand_i, where=where)
+                return cls(ref_axes, operand_i, dofdesc=dofdesc)
 
             return componentwise(make_op, operand)
         else:
             return DiscretizationProperty.__new__(cls)
 
-    def __init__(self, ref_axes, operand, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, ref_axes, operand, dofdesc=None):
         """
         :arg ref_axes: a :class:`tuple` of tuples indicating indices of
             coordinate axes of the reference element to the number of derivatives
@@ -551,7 +599,7 @@ class NumReferenceDerivative(DiscretizationProperty):
 
             May also be a singile integer *i*, which is viewed as equivalent
             to ``((i, 1),)``.
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
 
         if isinstance(ref_axes, int):
@@ -569,15 +617,16 @@ class NumReferenceDerivative(DiscretizationProperty):
         self.ref_axes = ref_axes
 
         self.operand = operand
-        DiscretizationProperty.__init__(self, where)
+        DiscretizationProperty.__init__(self, dofdesc)
 
     def __getinitargs__(self):
-        return (self.ref_axes, self.operand, self.where)
+        return (self.ref_axes, self.operand, self.dofdesc)
 
     mapper_method = intern("map_num_reference_derivative")
 
 
-def reference_jacobian(func, output_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def reference_jacobian(func, output_dim, dim, dofdesc=None):
     """Return a :class:`np.array` representing the Jacobian of a vector function
     with respect to the reference coordinates.
     """
@@ -586,35 +635,38 @@ def reference_jacobian(func, output_dim, dim, where=None):
     for i in range(output_dim):
         func_component = func[i]
         for j in range(dim):
-            jac[i, j] = NumReferenceDerivative(j, func_component, where)
+            jac[i, j] = NumReferenceDerivative(j, func_component, dofdesc)
 
     return jac
 
 
-def parametrization_derivative_matrix(ambient_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def parametrization_derivative_matrix(ambient_dim, dim, dofdesc=None):
     """Return a :class:`np.array` representing the derivative of the
     reference-to-global parametrization.
     """
 
     return cse(
             reference_jacobian(
-                [NodeCoordinateComponent(i, where) for i in range(ambient_dim)],
-                ambient_dim, dim, where=where),
+                [NodeCoordinateComponent(i, dofdesc) for i in range(ambient_dim)],
+                ambient_dim, dim, dofdesc=dofdesc),
             "pd_matrix", cse_scope.DISCRETIZATION)
 
 
-def parametrization_derivative(ambient_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def parametrization_derivative(ambient_dim, dim, dofdesc=None):
     """Return a :class:`pymbolic.geometric_algebra.MultiVector` representing
     the derivative of the reference-to-global parametrization.
     """
 
-    par_grad = parametrization_derivative_matrix(ambient_dim, dim, where)
+    par_grad = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     from pytools import product
     return product(MultiVector(vec) for vec in par_grad.T)
 
 
-def pseudoscalar(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def pseudoscalar(ambient_dim, dim=None, dofdesc=None):
     """
     Same as the outer product of all parametrization derivative columns.
     """
@@ -622,34 +674,37 @@ def pseudoscalar(ambient_dim, dim=None, where=None):
         dim = ambient_dim - 1
 
     return cse(
-            parametrization_derivative(ambient_dim, dim, where)
+            parametrization_derivative(ambient_dim, dim, dofdesc)
             .project_max_grade(),
             "pseudoscalar", cse_scope.DISCRETIZATION)
 
 
-def area_element(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def area_element(ambient_dim, dim=None, dofdesc=None):
     return cse(
-            sqrt(pseudoscalar(ambient_dim, dim, where).norm_squared()),
+            sqrt(pseudoscalar(ambient_dim, dim, dofdesc).norm_squared()),
             "area_element", cse_scope.DISCRETIZATION)
 
 
-def sqrt_jac_q_weight(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def sqrt_jac_q_weight(ambient_dim, dim=None, dofdesc=None):
     return cse(
             sqrt(
-                area_element(ambient_dim, dim, where)
-                * QWeight(where)),
+                area_element(ambient_dim, dim, dofdesc)
+                * QWeight(dofdesc)),
             "sqrt_jac_q_weight", cse_scope.DISCRETIZATION)
 
 
-def normal(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def normal(ambient_dim, dim=None, dofdesc=None):
     """Exterior unit normals."""
 
     # Don't be tempted to add a sign here. As it is, it produces
     # exterior normals for positively oriented curves and surfaces.
 
     pder = (
-            pseudoscalar(ambient_dim, dim, where)
-            / area_element(ambient_dim, dim, where))
+            pseudoscalar(ambient_dim, dim, dofdesc)
+            / area_element(ambient_dim, dim, dofdesc))
     return cse(
             # Dorst Section 3.7.2
             pder << pder.I.inv(),
@@ -657,7 +712,8 @@ def normal(ambient_dim, dim=None, where=None):
             scope=cse_scope.DISCRETIZATION)
 
 
-def mean_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def mean_curvature(ambient_dim, dim=None, dofdesc=None):
     """(Numerical) mean curvature."""
 
     if dim is None:
@@ -665,16 +721,16 @@ def mean_curvature(ambient_dim, dim=None, where=None):
 
     if ambient_dim == 2 and dim == 1:
         # https://en.wikipedia.org/wiki/Curvature#Local_expressions
-        xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where)
+        xp, yp = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
         xpp, ypp = cse(
-                reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where),
+                reference_jacobian([xp[0], yp[0]], ambient_dim, dim, dofdesc),
                 "p2d_matrix", cse_scope.DISCRETIZATION)
 
         kappa = (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)
     elif ambient_dim == 3 and dim == 2:
         # https://en.wikipedia.org/wiki/Mean_curvature#Surfaces_in_3D_space
-        s_op = shape_operator(ambient_dim, dim=dim, where=where)
+        s_op = shape_operator(ambient_dim, dim=dim, dofdesc=dofdesc)
         kappa = -0.5 * sum(s_op[i, i] for i in range(s_op.shape[0]))
     else:
         raise NotImplementedError('not available in {}D for {}D surfaces'
@@ -683,21 +739,23 @@ def mean_curvature(ambient_dim, dim=None, where=None):
     return kappa
 
 
-def first_fundamental_form(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def first_fundamental_form(ambient_dim, dim=None, dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
     if ambient_dim != 3 and dim != 2:
         raise NotImplementedError("only available for surfaces in 3D")
 
-    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     return cse(
             np.dot(pd_mat.T, pd_mat),
             "fundform1")
 
 
-def second_fundamental_form(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def second_fundamental_form(ambient_dim, dim=None, dofdesc=None):
     """Compute the second fundamental form of a surface. This is in reference
     to the reference-to-global mapping in use for each element.
 
@@ -712,17 +770,17 @@ def second_fundamental_form(ambient_dim, dim=None, where=None):
     if ambient_dim != 3 and dim != 2:
         raise NotImplementedError("only available for surfaces in 3D")
 
-    r = nodes(ambient_dim, where=where).as_vector()
+    r = nodes(ambient_dim, dofdesc=dofdesc).as_vector()
 
     # https://en.wikipedia.org/w/index.php?title=Second_fundamental_form&oldid=821047433#Classical_notation
 
     from functools import partial
-    d = partial(NumReferenceDerivative, where=where)
+    d = partial(NumReferenceDerivative, dofdesc=dofdesc)
     ruu = d(((0, 2),), r)
     ruv = d(((0, 1), (1, 1)), r)
     rvv = d(((1, 2),), r)
 
-    nrml = normal(ambient_dim, dim, where).as_vector()
+    nrml = normal(ambient_dim, dim, dofdesc).as_vector()
 
     ff2_l = cse(np.dot(ruu, nrml), "fundform2_L")
     ff2_m = cse(np.dot(ruv, nrml), "fundform2_M")
@@ -736,7 +794,8 @@ def second_fundamental_form(ambient_dim, dim=None, where=None):
     return result
 
 
-def shape_operator(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def shape_operator(ambient_dim, dim=None, dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
@@ -744,8 +803,8 @@ def shape_operator(ambient_dim, dim=None, where=None):
         raise NotImplementedError("only available for surfaces in 3D")
 
     # https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_surfaces&oldid=833587563
-    (E, F), (F, G) = first_fundamental_form(ambient_dim, dim, where)
-    (e, f), (f, g) = second_fundamental_form(ambient_dim, dim, where)
+    (E, F), (F, G) = first_fundamental_form(ambient_dim, dim, dofdesc)
+    (e, f), (f, g) = second_fundamental_form(ambient_dim, dim, dofdesc)
 
     result = np.zeros((2, 2), dtype=object)
     result[0, 0] = e*G-f*F
@@ -758,7 +817,8 @@ def shape_operator(ambient_dim, dim=None, where=None):
             "shape_operator")
 
 
-def _panel_size(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _panel_size(ambient_dim, dim=None, dofdesc=None):
     # A broken quasi-1D approximation of 1D element size. Do not use.
 
     if dim is None:
@@ -805,12 +865,13 @@ def _small_mat_eigenvalues(mat):
                 "eigenvalue formula for %dx%d matrices" % (m, n))
 
 
+@_deprecate_kwargs('where', 'dofdesc')
 def _equilateral_parametrization_derivative_matrix(ambient_dim, dim=None,
-        where=None):
+        dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
-    pder_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pder_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     # The above procedure works well only when the 'reference' end of the
     # mapping is in equilateral coordinates.
@@ -823,7 +884,8 @@ def _equilateral_parametrization_derivative_matrix(ambient_dim, dim=None,
             "equilateral_pder_mat")
 
 
-def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
+@_deprecate_kwargs('where', 'dofdesc')
+def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, dofdesc=None,
         with_elementwise_max=True):
     """Return the largest factor by which the reference-to-global
     mapping stretches the bi-unit (i.e. :math:`[-1,1]`) reference
@@ -845,7 +907,7 @@ def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
     # reflects available quadrature resolution in that direction.
 
     equi_pder_mat = _equilateral_parametrization_derivative_matrix(
-            ambient_dim, dim, where)
+            ambient_dim, dim, dofdesc)
 
     # Compute eigenvalues of J^T to compute SVD.
     equi_pder_mat_jtj = cse(
@@ -865,21 +927,22 @@ def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
     result = Max(tuple(stretch_factors))
 
     if with_elementwise_max:
-        result = ElementwiseMax(result, where=where)
+        result = ElementwiseMax(result, dofdesc=dofdesc)
 
     return cse(result, "mapping_max_stretch", cse_scope.DISCRETIZATION)
 
 
-def _max_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _max_curvature(ambient_dim, dim=None, dofdesc=None):
     # An attempt at a 'max curvature' criterion.
 
     if dim is None:
         dim = ambient_dim - 1
 
     if ambient_dim == 2:
-        return abs(mean_curvature(ambient_dim, dim, where=where))
+        return abs(mean_curvature(ambient_dim, dim, dofdesc=dofdesc))
     elif ambient_dim == 3:
-        shape_op = shape_operator(ambient_dim, dim, where=where)
+        shape_op = shape_operator(ambient_dim, dim, dofdesc=dofdesc)
 
         abs_principal_curvatures = [
                 abs(x) for x in _small_mat_eigenvalues(shape_op)]
@@ -890,15 +953,17 @@ def _max_curvature(ambient_dim, dim=None, where=None):
                 "dimensions" % ambient_dim)
 
 
-def _scaled_max_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _scaled_max_curvature(ambient_dim, dim=None, dofdesc=None):
     """An attempt at a unit-less, scale-invariant quantity that characterizes
     'how much curviness there is on an element'. Values seem to hover around 1
     on typical meshes. Empirical evidence suggests that elements exceeding
     a threshold of about 0.8-1 will have high QBX truncation error.
     """
 
-    return _max_curvature(ambient_dim, dim, where=where) * \
-            _simplex_mapping_max_stretch_factor(ambient_dim, dim, where=where,
+    return _max_curvature(ambient_dim, dim, dofdesc=dofdesc) * \
+            _simplex_mapping_max_stretch_factor(ambient_dim, dim,
+                    dofdesc=dofdesc,
                     with_elementwise_max=False)
 
 # }}}
@@ -914,24 +979,27 @@ def _expansion_radii_factor(ambient_dim, dim):
     return 0.5 * dim_fudge_factor
 
 
-def _quad_resolution(ambient_dim, dim=None, granularity=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _quad_resolution(ambient_dim, dim=None, granularity=None, dofdesc=None):
     """This measures the quadrature resolution across the
     mesh. In a 1D uniform mesh of uniform 'parametrization speed', it
     should be the same as the panel length.
     """
 
-    source = as_dofdesc(where)
+    source = as_dofdesc(dofdesc)
     target = source.copy(granularity=granularity)
 
     stretch = _simplex_mapping_max_stretch_factor(ambient_dim,
-            dim=dim, where=source)
+            dim=dim, dofdesc=source)
     return Interpolation(source, target, stretch)
 
 
-def _source_danger_zone_radii(ambient_dim, dim=None, granularity=None, where=None):
-    where = as_dofdesc(where)
-    if where.discr_stage is None:
-        where = where.copy(discr_stage=QBX_SOURCE_STAGE2)
+@_deprecate_kwargs('where', 'dofdesc')
+def _source_danger_zone_radii(ambient_dim, dim=None,
+        granularity=None, dofdesc=None):
+    dofdesc = as_dofdesc(dofdesc)
+    if dofdesc.discr_stage is None:
+        dofdesc = dofdesc.copy(discr_stage=QBX_SOURCE_STAGE2)
 
     # This should be the expression of the expansion radii, but
     #
@@ -946,32 +1014,34 @@ def _source_danger_zone_radii(ambient_dim, dim=None, granularity=None, where=Non
 
     factor = 0.75 * _expansion_radii_factor(ambient_dim, dim)
     return factor * _quad_resolution(ambient_dim, dim=dim,
-            granularity=granularity, where=where)
+            granularity=granularity, dofdesc=dofdesc)
 
 
-def _close_target_tunnel_radii(ambient_dim, dim=None, granularity=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _close_target_tunnel_radii(ambient_dim, dim=None,
+        granularity=None, dofdesc=None):
     factor = 0.5 * _expansion_radii_factor(ambient_dim, dim)
 
     return factor * _quad_resolution(ambient_dim, dim=dim,
-            granularity=granularity, where=where)
+            granularity=granularity, dofdesc=dofdesc)
 
 
-def expansion_radii(ambient_dim, dim=None, granularity=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def expansion_radii(ambient_dim, dim=None, granularity=None, dofdesc=None):
     factor = _expansion_radii_factor(ambient_dim, dim)
 
     return cse(factor * _quad_resolution(ambient_dim, dim=dim,
-        granularity=granularity, where=where),
+        granularity=granularity, dofdesc=dofdesc),
         "expansion_radii",
         cse_scope.DISCRETIZATION)
 
 
-def expansion_centers(ambient_dim, side, dim=None, where=None):
-    where = as_dofdesc(where)
-
-    x = nodes(ambient_dim, where=where)
-    normals = normal(ambient_dim, dim=dim, where=where)
+@_deprecate_kwargs('where', 'dofdesc')
+def expansion_centers(ambient_dim, side, dim=None, dofdesc=None):
+    x = nodes(ambient_dim, dofdesc=dofdesc)
+    normals = normal(ambient_dim, dim=dim, dofdesc=dofdesc)
     radii = expansion_radii(ambient_dim, dim=dim,
-            granularity=GRANULARITY_NODE, where=where)
+            granularity=GRANULARITY_NODE, dofdesc=dofdesc)
 
     centers = x + side * radii * normals
     return cse(centers.as_vector(),
@@ -979,18 +1049,20 @@ def expansion_centers(ambient_dim, side, dim=None, where=None):
             cse_scope.DISCRETIZATION)
 
 
-def h_max(ambient_dim, dim=None, where=None):
-    where = as_dofdesc(where).copy(granularity=GRANULARITY_ELEMENT)
-    r = _quad_resolution(ambient_dim, dim=dim, where=where)
+@_deprecate_kwargs('where', 'dofdesc')
+def h_max(ambient_dim, dim=None, dofdesc=None):
+    dofdesc = as_dofdesc(dofdesc).copy(granularity=GRANULARITY_ELEMENT)
+    r = _quad_resolution(ambient_dim, dim=dim, dofdesc=dofdesc)
 
     return cse(NodeMax(r),
             "h_max",
             cse_scope.DISCRETIZATION)
 
 
-def weights_and_area_elements(ambient_dim, dim=None, where=None):
-    area = area_element(ambient_dim, dim=dim, where=where)
-    return cse(area * QWeight(where=where),
+@_deprecate_kwargs('where', 'dofdesc')
+def weights_and_area_elements(ambient_dim, dim=None, dofdesc=None):
+    area = area_element(ambient_dim, dim=dim, dofdesc=dofdesc)
+    return cse(area * QWeight(dofdesc=dofdesc),
             "weights_area_elements",
             cse_scope.DISCRETIZATION)
 
@@ -1035,38 +1107,47 @@ class NodeMax(SingleScalarOperandExpression):
     mapper_method = "map_node_max"
 
 
-def integral(ambient_dim, dim, operand, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def integral(ambient_dim, dim, operand, dofdesc=None):
     """A volume integral of *operand*."""
 
     return NodeSum(
-            area_element(ambient_dim, dim, where)
-            * QWeight(where)
+            area_element(ambient_dim, dim, dofdesc)
+            * QWeight(dofdesc)
             * operand)
 
 
 class SingleScalarOperandExpressionWithWhere(Expression):
 
-    init_arg_names = ("operand", "where")
+    init_arg_names = ("operand", "dofdesc")
 
-    def __new__(cls, operand=None, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __new__(cls, operand=None, dofdesc=None):
         # If the constructor is handed a multivector object, return an
         # object array of the operator applied to each of the
         # coefficients in the multivector.
 
         if isinstance(operand, (np.ndarray, MultiVector)):
             def make_op(operand_i):
-                return cls(operand_i, where)
+                return cls(operand_i, dofdesc)
 
             return componentwise(make_op, operand)
         else:
             return Expression.__new__(cls)
 
-    def __init__(self, operand, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, operand, dofdesc=None):
         self.operand = operand
-        self.where = where
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.operand, self.where)
+        return (self.operand, self.dofdesc)
 
 
 class ElementwiseSum(SingleScalarOperandExpressionWithWhere):
@@ -1098,54 +1179,71 @@ class Ones(Expression):
     discretization.
     """
 
-    init_arg_names = ("where",)
+    init_arg_names = ("dofdesc",)
 
-    def __init__(self, where=None):
-        self.where = where
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, dofdesc=None):
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.where,)
+        return (self.dofdesc,)
 
     mapper_method = intern("map_ones")
 
 
-def ones_vec(dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def ones_vec(dim, dofdesc=None):
     from pytools.obj_array import make_obj_array
     return MultiVector(
-                make_obj_array(dim*[Ones(where)]))
+                make_obj_array(dim*[Ones(dofdesc)]))
 
 
-def area(ambient_dim, dim, where=None):
-    return cse(integral(ambient_dim, dim, Ones(where), where), "area",
+@_deprecate_kwargs('where', 'dofdesc')
+def area(ambient_dim, dim, dofdesc=None):
+    return cse(integral(ambient_dim, dim, Ones(dofdesc), dofdesc), "area",
             cse_scope.DISCRETIZATION)
 
 
-def mean(ambient_dim, dim, operand, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def mean(ambient_dim, dim, operand, dofdesc=None):
     return (
-            integral(ambient_dim, dim, operand, where)
-            / area(ambient_dim, dim, where))
+            integral(ambient_dim, dim, operand, dofdesc)
+            / area(ambient_dim, dim, dofdesc))
 
 
 class IterativeInverse(Expression):
 
-    init_arg_names = ("expression", "rhs", "variable_name", "extra_vars", "where")
+    init_arg_names = ("expression", "rhs", "variable_name", "extra_vars", "dofdesc")
 
+    @_deprecate_kwargs('where', 'dofdesc')
     def __init__(self, expression, rhs, variable_name, extra_vars={},
-            where=None):
+            dofdesc=None):
         self.expression = expression
         self.rhs = rhs
         self.variable_name = variable_name
         self.extra_vars = extra_vars
-        self.where = where
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
         return (self.expression, self.rhs, self.variable_name,
-                self.extra_vars, self.where)
+                self.extra_vars, self.dofdesc)
 
     def get_hash(self):
         return hash((self.__class__,) + (self.expression,
             self.rhs, self.variable_name,
-            frozenset(six.iteritems(self.extra_vars)), self.where))
+            frozenset(six.iteritems(self.extra_vars)), self.dofdesc))
 
     mapper_method = intern("map_inverse")
 
@@ -1533,10 +1631,11 @@ def S(kernel, density,
             kernel_arguments, **kwargs)
 
 
-def tangential_derivative(ambient_dim, operand, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_derivative(ambient_dim, operand, dim=None, dofdesc=None):
     pder = (
-            pseudoscalar(ambient_dim, dim, where)
-            / area_element(ambient_dim, dim, where))
+            pseudoscalar(ambient_dim, dim, dofdesc)
+            / area_element(ambient_dim, dim, dofdesc))
 
     # FIXME: Should be formula (3.25) in Dorst et al.
     d = Derivative()
@@ -1544,15 +1643,16 @@ def tangential_derivative(ambient_dim, operand, dim=None, where=None):
             (d.dnabla(ambient_dim) * d(operand)) >> pder)
 
 
-def normal_derivative(ambient_dim, operand, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def normal_derivative(ambient_dim, operand, dim=None, dofdesc=None):
     d = Derivative()
     return d.resolve(
-            (normal(ambient_dim, dim, where).scalar_product(d.dnabla(ambient_dim)))
+            (normal(ambient_dim, dim, dofdesc).scalar_product(d.dnabla(ambient_dim)))
             * d(operand))
 
 
 def Sp(kernel, *args, **kwargs):
-    where = kwargs.get("target")
+    dofdesc = kwargs.get("target")
     if "qbx_forced_limit" not in kwargs:
         warn("not specifying qbx_forced_limit on call to 'Sp' is deprecated, "
                 "defaulting to 'avg'", DeprecationWarning, stacklevel=2)
@@ -1570,7 +1670,7 @@ def Sp(kernel, *args, **kwargs):
     return normal_derivative(
             ambient_dim,
             S(kernel, *args, **kwargs),
-            dim=dim, where=where)
+            dim=dim, dofdesc=dofdesc)
 
 
 def Spp(kernel, *args, **kwargs):
@@ -1583,11 +1683,11 @@ def Spp(kernel, *args, **kwargs):
                 "the kernel, or directly")
     dim = kwargs.pop("dim", None)
 
-    where = kwargs.get("target")
+    dofdesc = kwargs.get("target")
     return normal_derivative(
             ambient_dim,
             Sp(kernel, *args, **kwargs),
-            dim=dim, where=where)
+            dim=dim, dofdesc=dofdesc)
 
 
 def D(kernel, *args, **kwargs):
@@ -1600,7 +1700,7 @@ def D(kernel, *args, **kwargs):
                 "the kernel, or directly")
     dim = kwargs.pop("dim", None)
 
-    where = kwargs.get("source")
+    dofdesc = kwargs.get("source")
 
     if "qbx_forced_limit" not in kwargs:
         warn("not specifying qbx_forced_limit on call to 'D' is deprecated, "
@@ -1609,7 +1709,7 @@ def D(kernel, *args, **kwargs):
 
     return int_g_dsource(
             ambient_dim,
-            normal(ambient_dim, dim, where),
+            normal(ambient_dim, dim, dofdesc),
             kernel, *args, **kwargs).xproject(0)
 
 
@@ -1630,7 +1730,7 @@ def Dp(kernel, *args, **kwargs):
     return normal_derivative(
             ambient_dim,
             D(kernel, *args, **kwargs),
-            dim=dim, where=target)
+            dim=dim, dofdesc=target)
 
 # }}}
 
@@ -1641,15 +1741,16 @@ def Dp(kernel, *args, **kwargs):
 
 # {{{ conventional vector calculus
 
-def tangential_onb(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_onb(ambient_dim, dim=None, dofdesc=None):
     """Return a matrix of shape ``(ambient_dim, dim)`` with orthogonal columns
-    spanning the tangential space of the surface of *where*.
+    spanning the tangential space of the surface of *dofdesc*.
     """
 
     if dim is None:
         dim = ambient_dim - 1
 
-    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     # {{{ Gram-Schmidt
 
@@ -1668,30 +1769,34 @@ def tangential_onb(ambient_dim, dim=None, where=None):
     return orth_pd_mat
 
 
-def xyz_to_tangential(xyz_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def xyz_to_tangential(xyz_vec, dofdesc=None):
     ambient_dim = len(xyz_vec)
-    tonb = tangential_onb(ambient_dim, where=where)
+    tonb = tangential_onb(ambient_dim, dofdesc=dofdesc)
     return make_obj_array([
         np.dot(tonb[:, i], xyz_vec)
         for i in range(ambient_dim - 1)
         ])
 
 
-def tangential_to_xyz(tangential_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_to_xyz(tangential_vec, dofdesc=None):
     ambient_dim = len(tangential_vec) + 1
-    tonb = tangential_onb(ambient_dim, where=where)
+    tonb = tangential_onb(ambient_dim, dofdesc=dofdesc)
     return sum(
         tonb[:, i] * tangential_vec[i]
         for i in range(ambient_dim - 1))
 
 
-def project_to_tangential(xyz_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def project_to_tangential(xyz_vec, dofdesc=None):
     return tangential_to_xyz(
-            cse(xyz_to_tangential(xyz_vec, where), where))
+            cse(xyz_to_tangential(xyz_vec, dofdesc), dofdesc))
 
 
-def n_dot(vec, where=None):
-    nrm = normal(len(vec), where).as_vector()
+@_deprecate_kwargs('where', 'dofdesc')
+def n_dot(vec, dofdesc=None):
+    nrm = normal(len(vec), dofdesc).as_vector()
 
     return np.dot(nrm, vec)
 
@@ -1708,8 +1813,9 @@ def cross(vec_a, vec_b):
         for i in range(3)])
 
 
-def n_cross(vec, where=None):
-    return cross(normal(3, where).as_vector(), vec)
+@_deprecate_kwargs('where', 'dofdesc')
+def n_cross(vec, dofdesc=None):
+    return cross(normal(3, dofdesc).as_vector(), vec)
 
 
 def div(vec):
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index 2af1edcc811b6c5c165a2b66fe0223dc083b3ecb..ca276959ab3ab51b6bec494059a97b0952894283 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -130,11 +130,11 @@ def run_source_refinement_test(ctx_factory, mesh, order, helmholtz_k=None):
         sym.expansion_radii(lpot_source.ambient_dim))(queue).get()
     source_danger_zone_radii = bind(lpot_source, sym._source_danger_zone_radii(
         lpot_source.ambient_dim,
-        where=sym.GRANULARITY_ELEMENT))(queue).get()
+        dofdesc=sym.GRANULARITY_ELEMENT))(queue).get()
 
     quad_res = bind(lpot_source, sym._quad_resolution(
         lpot_source.ambient_dim,
-        where=sym.GRANULARITY_ELEMENT))(queue)
+        dofdesc=sym.GRANULARITY_ELEMENT))(queue)
 
     # {{{ check if satisfying criteria
 
diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py
index a884e7da7e8ef73cadb9201613138430d6912fc9..2b3dbfa3705bb05813613a19eac3d7d72b4a64af 100644
--- a/test/test_scalar_int_eq.py
+++ b/test/test_scalar_int_eq.py
@@ -622,7 +622,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize):
         bc = bind(
                 (point_source, density_discr),
                 sym.normal_derivative(
-                    qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET)
+                    qbx.ambient_dim, pot_src, dofdesc=sym.DEFAULT_TARGET)
                 )(queue, charges=source_charges_dev, **concrete_knl_kwargs)
 
     # }}}