diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 71c8dfd2e58c7936bb08955be47a3757020a4f38..5dd3b22ec7d654cea909ea00a893dc5c7810fe53 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -728,6 +728,9 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     el_vertices = []
 
+    vertex_indices_left = {}
+    vertex_indices_right = {}
+
     if dim == 1:
         for i in range(shape[0]-1):
             # a--b
@@ -738,6 +741,11 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
             el_vertices.append((a, b,))
 
     elif dim == 2:
+        # For custom boundary tagging.
+        vertex_indices_left["x"] = vertex_indices[0, :].flatten()
+        vertex_indices_right["x"] = vertex_indices[-1, :].flatten()
+        vertex_indices_left["y"] = vertex_indices[:, 0].flatten()
+        vertex_indices_right["y"] = vertex_indices[:, -1].flatten()
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
 
@@ -757,6 +765,13 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
                     el_vertices.append((d, c, b))
 
     elif dim == 3:
+        # For custom boundary tagging.
+        vertex_indices_left["x"] = vertex_indices[0, :, :].flatten()
+        vertex_indices_right["x"] = vertex_indices[-1, :, :].flatten()
+        vertex_indices_left["y"] = vertex_indices[:, 0, :].flatten()
+        vertex_indices_right["y"] = vertex_indices[:, -1, :].flatten()
+        vertex_indices_left["z"] = vertex_indices[:, :, 0].flatten()
+        vertex_indices_right["z"] = vertex_indices[:, :, -1].flatten()
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
                 for k in range(shape[2]-1):
@@ -800,7 +815,6 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
     face_vertex_indices_to_tags = {}
     boundary_tags = list(boundary_tag_to_face.keys())
     axes = ["x", "y", "z", "w"]
-    face_ids = [1, 0, 0, 0]
     from meshmode.mesh import _compute_facial_adjacency_from_vertices
     if boundary_tags:
         for tag_idx, tag in enumerate(boundary_tags):
@@ -815,6 +829,18 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
                     raise ValueError("Face ", face, " given for boundary tag ",
                             tag, " is not formatted properly. Use +x, -x, +y, "
                             "+z, -z, +w, or -w")
+                if face.endswith("x"):
+                    pass
+                elif face.endswith("y"):
+                    pass
+                elif face.endswith("z"):
+                    pass
+                elif face.endswith("w"):
+                    pass
+                else:
+                    raise ValueError("Face ", face, " given for boundary tag ",
+                            tag, " is not formatted properly. Use +x, -x, +y, "
+                            "+z, -z, +w, or -w")
                 for i_ax, axis in enumerate(axes):
                     if face == "-" + axis:
                         if dim < i_ax + 1:
@@ -822,26 +848,25 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
                             "mismatch: facial dimension ", face, " for tag ",
                             tag, " is higher than the dimension of the "
                             "problem")
-                        face_id = face_ids[i_ax]
-                        dim_crit = i_ax
-                        node_crit = axis_coords[i_ax][0]
+                        vert_crit = vertex_indices_left[axis]
                     elif face == "+" + axis:
                         if dim < i_ax + 1:
                             raise ValueError("Boundary condition dimension "
                             "mismatch: facial dimension ", face, " for tag ",
                             tag, " is higher than the dimension of the "
                             "problem")
-                        face_id = face_ids[i_ax]
-                        dim_crit = i_ax
-                        node_crit = axis_coords[i_ax][-1]
+                        vert_crit = vertex_indices_right[axis]
                 for ielem in range(0, grp.nelements):
-                    if grp.nodes[dim_crit][ielem][0] == node_crit:
+                    # Loop through potential boundary-facing faces.
+                    for i in range(0, dim+1):
+                        # Element vertices:
                         fvi = grp.vertex_indices[ielem,
                                                  grp.face_vertex_indices()[
-                                                     face_id]]
-                        if frozenset(fvi) not in face_vertex_indices_to_tags:
-                            face_vertex_indices_to_tags[frozenset(fvi)] = []
-                        face_vertex_indices_to_tags[frozenset(fvi)].append(tag)
+                                                     i]]
+                        if all(ind in vert_crit for ind in fvi):
+                            if frozenset(fvi) not in face_vertex_indices_to_tags:
+                                face_vertex_indices_to_tags[frozenset(fvi)] = []
+                            face_vertex_indices_to_tags[frozenset(fvi)].append(tag)
         # Need to add BTAG ALL and BTAG_REALLY_ALL to list of tags
         from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
         boundary_tags.append(BTAG_ALL)
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 609e426b77558f063a0cd7f1d2708172c60e6386..a0248712a7f052d697cbc99b6d09caa81fea3013 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -126,18 +126,29 @@ def test_boundary_tags():
 # {{{ test custom boundary tags on box mesh
 @pytest.mark.parametrize(("dim", "nelem"), [
     (2, 20),
+    (3, 20),
     ])
 def test_box_boundary_tags(dim, nelem):
     from meshmode.mesh.generation import generate_regular_rect_mesh
     from meshmode.mesh import is_boundary_tag_empty
     from meshmode.mesh import check_bc_coverage
-    mesh = generate_regular_rect_mesh(a=(0, -1), b=(1, 1),
-                                      n=(nelem, nelem), order=3,
-                                      boundary_tag_to_face={
-                                      "btag_test_1": ["+x", "-y"],
-                                      "btag_test_2": ["+y", "-x"]})
+    if dim == 2:
+        a = (0, -1)
+        b = (1, 1)
+        n = (nelem, nelem)
+        btag_to_face = {"btag_test_1": ["+x", "-y"],
+                        "btag_test_2": ["+y", "-x"]}
+    elif dim == 3:
+        a = (0, -1, -1)
+        b = (1, 1, 1)
+        n = (nelem, nelem, nelem)
+        btag_to_face = {"btag_test_1": ["+x", "-y", "-z"],
+                        "btag_test_2": ["+y", "-x", "+z"]}
+    mesh = generate_regular_rect_mesh(a=a, b=b,
+                                      n=n, order=3,
+                                      boundary_tag_to_face=btag_to_face)
     # correct answer
-    num_on_bdy = dim*(nelem-1)
+    num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1)
 
     assert not is_boundary_tag_empty(mesh, "btag_test_1")
     assert not is_boundary_tag_empty(mesh, "btag_test_2")
@@ -159,6 +170,7 @@ def test_box_boundary_tags(dim, nelem):
                 num_marked_bdy_1 += 1
             if (-nbrs) & btag_2_bit:
                 num_marked_bdy_2 += 1
+
     # raise errors if wrong number of elements marked
     if num_marked_bdy_1 != num_on_bdy:
         raise ValueError("%i marked on custom boundary 1, should be %i" %