diff --git a/.editorconfig b/.editorconfig
new file mode 100644
index 0000000000000000000000000000000000000000..dcbc21d86f9e4b17ea7e8803d538c4c0f0b6276a
--- /dev/null
+++ b/.editorconfig
@@ -0,0 +1,32 @@
+# https://editorconfig.org/
+# https://github.com/editorconfig/editorconfig-vim 
+# https://github.com/editorconfig/editorconfig-emacs 
+
+root = true
+
+[*]
+indent_style = space
+end_of_line = lf
+charset = utf-8
+trim_trailing_whitespace = true
+insert_final_newline = true
+
+[*.py]
+indent_size = 4
+
+[*.rst]
+indent_size = 4
+
+[*.cpp]
+indent_size = 2
+
+[*.hpp]
+indent_size = 2
+
+# There may be one in doc/
+[Makefile]
+indent_style = tab
+
+# https://github.com/microsoft/vscode/issues/1679
+[*.md]
+trim_trailing_whitespace = false
diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index a96923dbb792ed9874595a3d84177a4f4d429635..df8ce3cefb07663abbb75529c7106c9e1873b35a 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -460,7 +460,7 @@ def dg_flux(c, tpair):
             np.dot(v.int, normal),
             u.int * normal) - flux_weak
 
-    return sym.interp(tpair.dd, "all_faces")(c*flux_strong)
+    return sym.project(tpair.dd, "all_faces")(c*flux_strong)
 
 
 def get_strong_wave_op_with_discr_direct(cl_ctx, dims=2, order=4):
@@ -499,8 +499,8 @@ def get_strong_wave_op_with_discr_direct(cl_ctx, dims=2, order=4):
 
     rad_normal = sym.normal(BTAG_ALL, dims)
 
-    rad_u = sym.cse(sym.interp("vol", BTAG_ALL)(u))
-    rad_v = sym.cse(sym.interp("vol", BTAG_ALL)(v))
+    rad_u = sym.cse(sym.project("vol", BTAG_ALL)(u))
+    rad_v = sym.cse(sym.project("vol", BTAG_ALL)(v))
 
     rad_bc = sym.cse(join_fields(
             0.5*(rad_u - sign*np.dot(rad_normal, rad_v)),
@@ -627,7 +627,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper):
             profile_data["bytes_written"] = (
                     profile_data.get("bytes_written", 0) + result.nbytes)
 
-            if op.mapper_method == "map_interpolation":
+            if op.mapper_method == "map_projection":
                 profile_data["interp_bytes_read"] = (
                         profile_data.get("interp_bytes_read", 0) + field.nbytes)
                 profile_data["interp_bytes_written"] = (
@@ -649,7 +649,7 @@ class ExecutionMapperWithMemOpCounting(ExecutionMapperWrapper):
                     expr.op,
                     (
                         # TODO: Not comprehensive.
-                        sym_op.InterpolationOperator,
+                        sym_op.ProjectionOperator,
                         sym_op.RefFaceMassOperator,
                         sym_op.RefMassOperator,
                         sym_op.RefInverseMassOperator,
diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py
index aa9e87ffa0a7885adf9cc8d67d82d1be15a44e93..f060c364b558273d7ab97d678f021ac6e30eabd6 100644
--- a/examples/wave/wave-eager.py
+++ b/examples/wave/wave-eager.py
@@ -38,7 +38,7 @@ from grudge.shortcuts import make_visualizer
 
 
 def interior_trace_pair(discr, vec):
-    i = discr.interp("vol", "int_faces", vec)
+    i = discr.project("vol", "int_faces", vec)
     e = with_object_array_or_scalar(
             lambda el: discr.opposite_face_connection()(el.queue, el),
             i)
@@ -69,15 +69,15 @@ def wave_flux(discr, c, w_tpair):
             0.5*normal_times(v_jump),
             )
 
-    return discr.interp(w_tpair.dd, "all_faces", c*flux_weak)
+    return discr.project(w_tpair.dd, "all_faces", c*flux_weak)
 
 
 def wave_operator(discr, c, w):
     u = w[0]
     v = w[1:]
 
-    dir_u = discr.interp("vol", BTAG_ALL, u)
-    dir_v = discr.interp("vol", BTAG_ALL, v)
+    dir_u = discr.project("vol", BTAG_ALL, u)
+    dir_v = discr.project("vol", BTAG_ALL, v)
     dir_bval = join_fields(dir_u, dir_v)
     dir_bc = join_fields(-dir_u, dir_v)
 
diff --git a/grudge/eager.py b/grudge/eager.py
index e64e2fa8344a848f4d0a218cfd0286ec9f896e2a..2000ef0cd0de44f6d7d74b3061179d0a9016dc93 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -46,9 +46,16 @@ def without_queue(ary):
 
 class EagerDGDiscretization(DGDiscretizationWithBoundaries):
     def interp(self, src, tgt, vec):
+        from warnings import warn
+        warn("using 'interp' is deprecated, use 'project' instead.",
+                DeprecationWarning, stacklevel=2)
+
+        return self.project(src, tgt, vec)
+
+    def project(self, src, tgt, vec):
         if is_obj_array(vec):
             return with_object_array_or_scalar(
-                    lambda el: self.interp(src, tgt, el), vec)
+                    lambda el: self.project(src, tgt, el), vec)
 
         return self.connection_from_dds(src, tgt)(vec.queue, vec)
 
diff --git a/grudge/execution.py b/grudge/execution.py
index b3faae20d4cc6862f81f3d2c2071ad1e53fa8a8b..8a4f777a4b0bf2abba8b27ca1e7cbc061ad9ab67 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -250,7 +250,7 @@ class ExecutionMapper(mappers.Evaluator,
 
         return result
 
-    def map_interpolation(self, op, field_expr):
+    def map_projection(self, op, field_expr):
         conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out)
         return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
 
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 88c09668bf6d61fb26ca4349f075843d3ba76495..659c10f8519f1397aa8609cb94401e1bbc257d7c 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -94,7 +94,7 @@ class StrongAdvectionOperator(AdvectionOperatorBase):
         u = sym.var("u")
 
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(
+            return sym.project(pair.dd, "all_faces")(
                     self.flux(pair))
 
         return (
@@ -118,11 +118,11 @@ class WeakAdvectionOperator(AdvectionOperatorBase):
         u = sym.var("u")
 
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(
+            return sym.project(pair.dd, "all_faces")(
                     self.flux(pair))
 
         bc_in = self.inflow_u
-        # bc_out = sym.interp("vol", self.outflow_tag)(u)
+        # bc_out = sym.project(sym.DD_VOLUME, self.outflow_tag)(u)
 
         return sym.InverseMassOperator()(
                 np.dot(
@@ -149,20 +149,20 @@ class VariableCoefficientAdvectionOperator(AdvectionOperatorBase):
         self.quad_tag = quad_tag
 
     def flux(self, u):
-        surf_v = sym.interp(sym.DD_VOLUME, u.dd)(self.v)
+        surf_v = sym.project(sym.DD_VOLUME, u.dd)(self.v)
         return advection_weak_flux(self.flux_type, u, surf_v)
 
     def sym_operator(self):
         u = sym.var("u")
 
         def flux(pair):
-            return sym.interp(pair.dd, face_dd)(self.flux(pair))
+            return sym.project(pair.dd, face_dd)(self.flux(pair))
 
         face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
         boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
         quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag)
 
-        to_quad = sym.interp(sym.DD_VOLUME, quad_dd)
+        to_quad = sym.project(sym.DD_VOLUME, quad_dd)
         stiff_t_op = sym.stiffness_t(self.ambient_dim,
                 dd_in=quad_dd, dd_out=sym.DD_VOLUME)
 
@@ -227,26 +227,26 @@ class SurfaceAdvectionOperator(AdvectionOperatorBase):
         self.quad_tag = quad_tag
 
     def flux(self, u):
-        surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
+        surf_v = sym.project(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
         return surface_advection_weak_flux(self.flux_type, u, surf_v)
 
     def penalty(self, u):
-        surf_v = sym.interp(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
+        surf_v = sym.project(sym.DD_VOLUME, u.dd.with_qtag(None))(self.v)
         return surface_penalty_flux(u, surf_v, tau=0.0)
 
     def sym_operator(self):
         u = sym.var("u")
 
         def flux(pair):
-            return sym.interp(pair.dd, face_dd)(self.flux(pair))
+            return sym.project(pair.dd, face_dd)(self.flux(pair))
 
         def penalty(pair):
-            return sym.interp(pair.dd, face_dd)(self.penalty(pair))
+            return sym.project(pair.dd, face_dd)(self.penalty(pair))
 
         face_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
         quad_dd = sym.DOFDesc(sym.DTAG_VOLUME_ALL, self.quad_tag)
 
-        to_quad = sym.interp(sym.DD_VOLUME, quad_dd)
+        to_quad = sym.project(sym.DD_VOLUME, quad_dd)
         stiff_t_op = sym.stiffness_t(self.ambient_dim,
                 dd_in=quad_dd, dd_out=sym.DD_VOLUME)
 
diff --git a/grudge/models/em.py b/grudge/models/em.py
index 225b3b4a411beae96359df0548ff5366f00ea675..f564196711a1b2749a1987d93cb6f2eeaf85eb74 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -191,8 +191,8 @@ class MaxwellOperator(HyperbolicOperator):
         "Construct part of the flux operator template for PEC boundary conditions"
         e, h = self.split_eh(w)
 
-        pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e))
-        pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h))
+        pec_e = sym.cse(sym.project("vol", self.pec_tag)(e))
+        pec_h = sym.cse(sym.project("vol", self.pec_tag)(h))
 
         return join_fields(-pec_e, pec_h)
 
@@ -200,8 +200,8 @@ class MaxwellOperator(HyperbolicOperator):
         "Construct part of the flux operator template for PMC boundary conditions"
         e, h = self.split_eh(w)
 
-        pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e))
-        pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h))
+        pmc_e = sym.cse(sym.project("vol", self.pmc_tag)(e))
+        pmc_h = sym.cse(sym.project("vol", self.pmc_tag)(h))
 
         return join_fields(pmc_e, -pmc_h)
 
@@ -221,8 +221,8 @@ class MaxwellOperator(HyperbolicOperator):
         absorb_Z = (mu/epsilon)**0.5  # noqa: N806
         absorb_Y = 1/absorb_Z  # noqa: N806
 
-        absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e))
-        absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h))
+        absorb_e = sym.cse(sym.project("vol", self.absorb_tag)(e))
+        absorb_h = sym.cse(sym.project("vol", self.absorb_tag)(h))
 
         bc = join_fields(
                 absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
@@ -278,7 +278,7 @@ class MaxwellOperator(HyperbolicOperator):
                 ]
 
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(self.flux(pair))
+            return sym.project(pair.dd, "all_faces")(self.flux(pair))
 
         return (
                 - self.local_derivatives(w)
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 6defafc878f856fe83449b12dc055cfdf0b1f432..116cd329b0d648722b2537bb99545d783301b7b0 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -113,8 +113,8 @@ class StrongWaveOperator(HyperbolicOperator):
         # boundary conditions -------------------------------------------------
 
         # dirichlet BCs -------------------------------------------------------
-        dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
-        dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
+        dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
+        dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
         if self.dirichlet_bc_f:
             # FIXME
             from warnings import warn
@@ -129,15 +129,15 @@ class StrongWaveOperator(HyperbolicOperator):
         dir_bc = sym.cse(dir_bc, "dir_bc")
 
         # neumann BCs ---------------------------------------------------------
-        neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
-        neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+        neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
+        neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
         neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
 
         # radiation BCs -------------------------------------------------------
         rad_normal = sym.normal(self.radiation_tag, d)
 
-        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
-        rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
+        rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
+        rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
 
         rad_bc = sym.cse(join_fields(
                 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
@@ -148,7 +148,7 @@ class StrongWaveOperator(HyperbolicOperator):
         nabla = sym.nabla(d)
 
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(
+            return sym.project(pair.dd, "all_faces")(
                     self.flux(pair))
 
         result = (
@@ -253,8 +253,8 @@ class WeakWaveOperator(HyperbolicOperator):
         # boundary conditions -------------------------------------------------
 
         # dirichlet BCs -------------------------------------------------------
-        dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
-        dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
+        dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
+        dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
         if self.dirichlet_bc_f:
             # FIXME
             from warnings import warn
@@ -269,15 +269,15 @@ class WeakWaveOperator(HyperbolicOperator):
         dir_bc = sym.cse(dir_bc, "dir_bc")
 
         # neumann BCs ---------------------------------------------------------
-        neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
-        neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+        neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
+        neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
         neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
 
         # radiation BCs -------------------------------------------------------
         rad_normal = sym.normal(self.radiation_tag, d)
 
-        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
-        rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
+        rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
+        rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
 
         rad_bc = sym.cse(join_fields(
                 0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
@@ -286,7 +286,7 @@ class WeakWaveOperator(HyperbolicOperator):
 
         # entire operator -----------------------------------------------------
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(self.flux(pair))
+            return sym.project(pair.dd, "all_faces")(self.flux(pair))
 
         result = sym.InverseMassOperator()(
                 join_fields(
@@ -395,9 +395,9 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         # boundary conditions -------------------------------------------------
 
         # dirichlet BCs -------------------------------------------------------
-        dir_c = sym.cse(sym.interp("vol", self.dirichlet_tag)(self.c))
-        dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
-        dir_v = sym.cse(sym.interp("vol", self.dirichlet_tag)(v))
+        dir_c = sym.cse(sym.project("vol", self.dirichlet_tag)(self.c))
+        dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
+        dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
         if self.dirichlet_bc_f:
             # FIXME
             from warnings import warn
@@ -412,28 +412,28 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         dir_bc = sym.cse(dir_bc, "dir_bc")
 
         # neumann BCs ---------------------------------------------------------
-        neu_c = sym.cse(sym.interp("vol", self.neumann_tag)(self.c))
-        neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
-        neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+        neu_c = sym.cse(sym.project("vol", self.neumann_tag)(self.c))
+        neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
+        neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
         neu_bc = sym.cse(join_fields(neu_c, neu_u, -neu_v), "neu_bc")
 
         # radiation BCs -------------------------------------------------------
         rad_normal = sym.normal(self.radiation_tag, d)
 
-        rad_c = sym.cse(sym.interp("vol", self.radiation_tag)(self.c))
-        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
-        rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
+        rad_c = sym.cse(sym.project("vol", self.radiation_tag)(self.c))
+        rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
+        rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
 
         rad_bc = sym.cse(join_fields(rad_c,
-                0.5*(rad_u - sym.interp("vol", self.radiation_tag)(self.sign)
+                0.5*(rad_u - sym.project("vol", self.radiation_tag)(self.sign)
                     * np.dot(rad_normal, rad_v)),
                 0.5*rad_normal*(np.dot(rad_normal, rad_v)
-                    - sym.interp("vol", self.radiation_tag)(self.sign)*rad_u)
+                    - sym.project("vol", self.radiation_tag)(self.sign)*rad_u)
                 ), "rad_bc")
 
         # entire operator -----------------------------------------------------
         def flux(pair):
-            return sym.interp(pair.dd, "all_faces")(self.flux(pair))
+            return sym.project(pair.dd, "all_faces")(self.flux(pair))
 
         result = sym.InverseMassOperator()(
                 join_fields(
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 7b36ee9f73a7fdf0143619a06617789f7a8a7862..abd1d78b6bf3bf27b7b65cb3e7c14a256ba3a162 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -95,6 +95,7 @@ class LoopyKernelInstruction(Instruction):
     scope_indicator = ""
 
     def __init__(self, kernel_descriptor):
+        super(LoopyKernelInstruction, self).__init__()
         self.kernel_descriptor = kernel_descriptor
 
     @memoize_method
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index c8def89882ea06b1002135983b9a1cdf570ea5db..e8a63c1c697ae964689459ce1a056ef7fc0c8d0e 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -131,7 +131,7 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
 
     map_elementwise_linear = _map_op_base
 
-    map_interpolation = _map_op_base
+    map_projection = _map_op_base
 
     map_elementwise_sum = _map_op_base
     map_elementwise_min = _map_op_base
@@ -184,7 +184,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
         # it's a leaf--no changing children
         return expr
 
-    map_interpolation = map_elementwise_linear
+    map_projection = map_elementwise_linear
 
     map_elementwise_sum = map_elementwise_linear
     map_elementwise_min = map_elementwise_linear
@@ -358,9 +358,9 @@ class DOFDescReplacer(IdentityMapper):
             field = self.rec(expr.field)
             return op.OppositePartitionFaceSwap(dd_in=self.new_dd,
                                                 dd_out=self.new_dd)(field)
-        elif (isinstance(expr.op, op.InterpolationOperator)
+        elif (isinstance(expr.op, op.ProjectionOperator)
                     and expr.op.dd_out == self.prev_dd):
-            return op.InterpolationOperator(dd_in=expr.op.dd_in,
+            return op.ProjectionOperator(dd_in=expr.op.dd_in,
                                             dd_out=self.new_dd)(expr.field)
         elif (isinstance(expr.op, op.RefDiffOperatorBase)
                     and expr.op.dd_out == self.prev_dd
@@ -420,14 +420,14 @@ class DistributedMapper(CSECachingMapperMixin, IdentityMapper):
         from meshmode.mesh import BTAG_PARTITION
         from meshmode.discretization.connection import (FACE_RESTR_ALL,
                                                         FACE_RESTR_INTERIOR)
-        if (isinstance(expr.op, op.InterpolationOperator)
+        if (isinstance(expr.op, op.ProjectionOperator)
                 and expr.op.dd_in.domain_tag is FACE_RESTR_INTERIOR
                 and expr.op.dd_out.domain_tag is FACE_RESTR_ALL):
             distributed_work = 0
             for i_remote_part in self.connected_parts:
                 mapped_field = RankGeometryChanger(i_remote_part)(expr.field)
                 btag_part = BTAG_PARTITION(i_remote_part)
-                distributed_work += op.InterpolationOperator(dd_in=btag_part,
+                distributed_work += op.ProjectionOperator(dd_in=btag_part,
                                              dd_out=expr.op.dd_out)(mapped_field)
             return expr + distributed_work
         else:
@@ -458,9 +458,9 @@ class RankGeometryChanger(CSECachingMapperMixin, IdentityMapper):
                     dd_in=self.new_dd,
                     dd_out=self.new_dd,
                     unique_id=expr.op.unique_id)(field)
-        elif (isinstance(expr.op, op.InterpolationOperator)
+        elif (isinstance(expr.op, op.ProjectionOperator)
                     and expr.op.dd_out == self.prev_dd):
-            return op.InterpolationOperator(dd_in=expr.op.dd_in,
+            return op.ProjectionOperator(dd_in=expr.op.dd_in,
                                             dd_out=self.new_dd)(expr.field)
         elif (isinstance(expr.op, op.RefDiffOperator)
                     and expr.op.dd_out == self.prev_dd
@@ -850,8 +850,8 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
     def map_function_symbol(self, expr, enclosing_prec):
         return expr.name
 
-    def map_interpolation(self, expr, enclosing_prec):
-        return "Interp" + self._format_op_dd(expr)
+    def map_projection(self, expr, enclosing_prec):
+        return "Project" + self._format_op_dd(expr)
 
 
 class PrettyStringifyMapper(
@@ -908,14 +908,14 @@ class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
         import grudge.symbolic.operators as op
         # changed
 
-        if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator):
+        if dd_in == dd_out and isinstance(expr.op, op.ProjectionOperator):
             # This used to be to-quad interpolation and has become a no-op.
             # Remove it.
             return self.rec(expr.field)
 
         if isinstance(expr.op, op.StiffnessTOperator):
             new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out)
-        elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)):
+        elif isinstance(expr.op, (op.FaceMassOperator, op.ProjectionOperator)):
             new_op = type(expr.op)(dd_in, dd_out)
         else:
             raise NotImplementedError("do not know how to modify dd_in and dd_out "
@@ -978,7 +978,7 @@ class EmptyFluxKiller(CSECachingMapperMixin, IdentityMapper):
 
     def map_operator_binding(self, expr):
         from meshmode.mesh import is_boundary_tag_empty
-        if (isinstance(expr.op, sym.InterpolationOperator)
+        if (isinstance(expr.op, sym.ProjectionOperator)
                 and expr.op.dd_out.is_boundary_or_partition_interface()):
             domain_tag = expr.op.dd_out.domain_tag
             assert isinstance(domain_tag, sym.DTAG_BOUNDARY)
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 8f7609c016bc88ca1f2e065dd7dc97021a5e44f8..0324592367aa48da0dd47aa7a90dea58aee9fad9 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -36,9 +36,9 @@ Basic Operators
 
 .. autoclass:: Operator
 .. autoclass:: ElementwiseLinearOperator
-.. autoclass:: InterpolationOperator
+.. autoclass:: ProjectionOperator
 
-.. data:: interp
+.. data:: project
 
 Reductions
 ^^^^^^^^^^
@@ -146,17 +146,17 @@ class ElementwiseLinearOperator(Operator):
     mapper_method = intern("map_elementwise_linear")
 
 
-class InterpolationOperator(Operator):
+class ProjectionOperator(Operator):
     def __init__(self, dd_in, dd_out):
-        super(InterpolationOperator, self).__init__(dd_in, dd_out)
+        super(ProjectionOperator, self).__init__(dd_in, dd_out)
 
     def __call__(self, expr):
         from pytools.obj_array import with_object_array_or_scalar
 
-        def interp_one(subexpr):
+        def project_one(subexpr):
             from pymbolic.primitives import is_constant
             if self.dd_in == self.dd_out:
-                # no-op interpolation, go away
+                # no-op projection, go away
                 return subexpr
             elif is_constant(subexpr):
                 return subexpr
@@ -164,12 +164,30 @@ class InterpolationOperator(Operator):
                 from grudge.symbolic.primitives import OperatorBinding
                 return OperatorBinding(self, subexpr)
 
-        return with_object_array_or_scalar(interp_one, expr)
+        return with_object_array_or_scalar(project_one, expr)
+
+    mapper_method = intern("map_projection")
+
+
+project = ProjectionOperator
+
+
+class InterpolationOperator(ProjectionOperator):
+    def __init__(self, dd_in, dd_out):
+        from warnings import warn
+        warn("'InterpolationOperator' is deprecated, "
+                "use 'ProjectionOperator' instead.",
+                DeprecationWarning, stacklevel=2)
+
+        super(InterpolationOperator, self).__init__(dd_in, dd_out)
 
-    mapper_method = intern("map_interpolation")
 
+def interp(dd_in, dd_out):
+    from warnings import warn
+    warn("using 'interp' is deprecated, use 'project' instead.",
+            DeprecationWarning, stacklevel=2)
 
-interp = InterpolationOperator
+    return ProjectionOperator(dd_in, dd_out)
 
 
 # {{{ element reduction: sum, min, max
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 7e26223ba401aee32d60a0071e33ac1b0c782b68..5d146c323595b55d56d0e42115d85d1d88ce38e5 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -528,8 +528,8 @@ def forward_metric_nth_derivative(xyz_axis, ref_axes, dd=None):
         result = RefDiffOperator(rst_axis, inner_dd)(result)
 
     if dd.uses_quadrature():
-        from grudge.symbolic.operators import interp
-        result = interp(inner_dd, dd)(result)
+        from grudge.symbolic.operators import project
+        result = project(inner_dd, dd)(result)
 
     prefix = "dx%d_%s" % (
             xyz_axis,
@@ -773,10 +773,10 @@ def mv_normal(dd, ambient_dim, dim=None):
     # correct exterior face normal vector.
     assert dim == ambient_dim - 2
 
-    from grudge.symbolic.operators import interp
+    from grudge.symbolic.operators import project
     volm_normal = (
             surface_normal(ambient_dim, dim=dim + 1, dd=DD_VOLUME)
-            .map(interp(DD_VOLUME, dd)))
+            .map(project(DD_VOLUME, dd)))
     pder = pseudoscalar(ambient_dim, dim, dd=dd)
 
     mv = cse(-(volm_normal ^ pder) << volm_normal.I.inv(),
@@ -868,7 +868,7 @@ class TracePair:
 
 
 def int_tpair(expression, qtag=None, from_dd=None):
-    from grudge.symbolic.operators import interp, OppositeInteriorFaceSwap
+    from grudge.symbolic.operators import project, OppositeInteriorFaceSwap
 
     if from_dd is None:
         from_dd = DD_VOLUME
@@ -878,12 +878,12 @@ def int_tpair(expression, qtag=None, from_dd=None):
     if from_dd.domain_tag == trace_dd.domain_tag:
         i = expression
     else:
-        i = interp(from_dd, trace_dd.with_qtag(None))(expression)
+        i = project(from_dd, trace_dd.with_qtag(None))(expression)
     e = cse(OppositeInteriorFaceSwap()(i))
 
     if trace_dd.uses_quadrature():
-        i = cse(interp(trace_dd.with_qtag(None), trace_dd)(i))
-        e = cse(interp(trace_dd.with_qtag(None), trace_dd)(e))
+        i = cse(project(trace_dd.with_qtag(None), trace_dd)(i))
+        e = cse(project(trace_dd.with_qtag(None), trace_dd)(e))
 
     return TracePair(trace_dd, i, e)
 
@@ -909,8 +909,8 @@ def bv_tpair(dd, interior, exterior):
         representing the exterior value to be used
         for the flux.
     """
-    from grudge.symbolic.operators import interp
-    interior = cse(interp("vol", dd)(interior))
+    from grudge.symbolic.operators import project
+    interior = cse(project("vol", dd)(interior))
     return TracePair(dd, interior, exterior)
 
 # }}}
diff --git a/setup.py b/setup.py
index 4c2b718dc4cf873e06b2b9c47845c7de5dd50d08..59dcb5c85c39006f335642f400b9f991778f985f 100644
--- a/setup.py
+++ b/setup.py
@@ -42,6 +42,7 @@ def main():
 
           packages=find_packages(),
 
+          python_requires="~=3.6",
           install_requires=[
               "pytest>=2.3",
               "pytools>=2018.5.2",
diff --git a/test/test_grudge.py b/test/test_grudge.py
index e27fe03bfe6c2c039302e909e18009ed3f6b1ef1..371b0f19b1c539950bd3ebef13371c2c4a4278f6 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -376,7 +376,7 @@ def test_face_normal_surface(ctx_factory, mesh_name):
     ambient_dim = mesh.ambient_dim
     dim = mesh.dim
 
-    sym_surf_normal = sym.interp(dv, df)(
+    sym_surf_normal = sym.project(dv, df)(
             sym.surface_normal(ambient_dim, dim=dim, dd=dv).as_vector()
             )
     sym_surf_normal = sym_surf_normal / sym.sqrt(sum(sym_surf_normal**2))
@@ -528,7 +528,7 @@ def test_2d_gauss_theorem(ctx_factory):
                 ).sum())
             -  # noqa: W504
             sym.integral(
-                sym.interp("vol", sym.BTAG_ALL)(f(sym.nodes(2)))
+                sym.project("vol", sym.BTAG_ALL)(f(sym.nodes(2)))
                 .dot(sym.normal(sym.BTAG_ALL, 2)),
                 dd=sym.BTAG_ALL)
             )(queue)
@@ -636,7 +636,7 @@ def test_surface_divergence_theorem(ctx_factory, mesh_name, visualize=False):
         sym_normal = sym.surface_normal(ambient_dim, dim=dim, dd=dq).as_vector()
 
         sym_face_normal = sym.normal(df, ambient_dim, dim=dim - 1)
-        sym_face_f = sym.interp(dd, df)(sym_f)
+        sym_face_f = sym.project(dd, df)(sym_f)
 
         # operators
         sym_stiff = sum(
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index f9b12d3b4ec8922204cb7201c715adbb51f037b4..0b5ae492798b2783564f4afc75ec3051a43c3461 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -70,15 +70,15 @@ def simple_mpi_communication_entrypoint():
     myfunc = bind(vol_discr, myfunc_symb)(queue)
 
     sym_all_faces_func = sym.cse(
-        sym.interp("vol", "all_faces")(sym.var("myfunc")))
+        sym.project("vol", "all_faces")(sym.var("myfunc")))
     sym_int_faces_func = sym.cse(
-        sym.interp("vol", "int_faces")(sym.var("myfunc")))
+        sym.project("vol", "int_faces")(sym.var("myfunc")))
     sym_bdry_faces_func = sym.cse(
-        sym.interp(sym.BTAG_ALL, "all_faces")(
-            sym.interp("vol", sym.BTAG_ALL)(sym.var("myfunc"))))
+        sym.project(sym.BTAG_ALL, "all_faces")(
+            sym.project("vol", sym.BTAG_ALL)(sym.var("myfunc"))))
 
     bound_face_swap = bind(vol_discr,
-        sym.interp("int_faces", "all_faces")(
+        sym.project("int_faces", "all_faces")(
             sym.OppositeInteriorFaceSwap("int_faces")(
                 sym_int_faces_func)
             ) - (sym_all_faces_func - sym_bdry_faces_func)