diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index ddbd7d65..c446fceb 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -1230,7 +1230,12 @@ def show_scalar_in_matplotlib_3d(self, field, **kwargs):
             while len(nodes) < 3:
                 nodes.append(0*nodes[0])
 
-            from matplotlib.tri.triangulation import Triangulation
+            try:
+                from matplotlib.tri import Triangulation
+            except ImportError:
+                # NOTE: deprecated starting with v3.7
+                from matplotlib.tri.triangulation import Triangulation
+
             tri, _, kwargs = \
                 Triangulation.get_from_args_and_kwargs(
                         *nodes,
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 3f52ec51..3c2d16df 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1092,6 +1092,7 @@ class Mesh:
     .. autoattribute:: _nodal_adjacency
     .. autoattribute:: _facial_adjacency_groups
 
+    .. automethod:: copy
     .. automethod:: __eq__
     """
 
@@ -1242,12 +1243,15 @@ def __init__(
                 node_vertex_consistency_tolerance=node_vertex_consistency_tolerance,
                 skip_element_orientation_test=skip_element_orientation_test)
 
-    def copy(self, **kwargs: Any) -> "Mesh":
-        warn(f"'{type(self).__name__}.copy' is deprecated and will be removed in "
-             f"2025. '{type(self).__name__}' is a dataclass and can use the "
-             "standard 'replace' function.",
-             DeprecationWarning, stacklevel=2)
-
+    def copy(self, *,
+             skip_tests: bool = False,
+             node_vertex_consistency_tolerance:
+                 Optional[Union[Literal[False], bool]] = None,
+             skip_element_orientation_test: bool = False,
+             # NOTE: this is set to *True* to avoid the meaningless warning in
+             # `__init__` when calling `Mesh.copy`
+             factory_constructed: bool = True,
+             **kwargs: Any) -> "Mesh":
         if "nodal_adjacency" in kwargs:
             kwargs["_nodal_adjacency"] = kwargs.pop("nodal_adjacency")
 
@@ -1255,7 +1259,14 @@ def copy(self, **kwargs: Any) -> "Mesh":
             kwargs["_facial_adjacency_groups"] = (
                 kwargs.pop("facial_adjacency_groups"))
 
-        return replace(self, **kwargs)
+        mesh = replace(self, factory_constructed=factory_constructed, **kwargs)
+        if __debug__ and not skip_tests:
+            check_mesh_consistency(
+                mesh,
+                node_vertex_consistency_tolerance=node_vertex_consistency_tolerance,
+                skip_element_orientation_test=skip_element_orientation_test)
+
+        return mesh
 
     @property
     def ambient_dim(self) -> int:
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index a61349da..1b2a64dd 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -1359,7 +1359,10 @@ def glue_mesh_boundaries(
 
     return mesh.copy(
         nodal_adjacency=False,
-        facial_adjacency_groups=facial_adjacency_groups)
+        _facial_adjacency_groups=tuple([
+            tuple(fagrps) for fagrps in facial_adjacency_groups
+            ]),
+        )
 
 # }}}
 
@@ -1403,7 +1406,8 @@ def map_mesh(mesh: Mesh, f: Callable[[np.ndarray], np.ndarray]) -> Mesh:
     # }}}
 
     return mesh.copy(
-            vertices=vertices, groups=new_groups,
+            vertices=vertices,
+            groups=tuple(new_groups),
             is_conforming=mesh.is_conforming)
 
 # }}}
@@ -1499,8 +1503,11 @@ def compute_new_map(old_map: AffineMap) -> AffineMap:
     # }}}
 
     return mesh.copy(
-            vertices=vertices, groups=new_groups,
-            facial_adjacency_groups=facial_adjacency_groups,
+            vertices=vertices,
+            groups=tuple(new_groups),
+            _facial_adjacency_groups=tuple([
+                tuple(fagrps) for fagrps in facial_adjacency_groups
+                ]) if facial_adjacency_groups is not None else None,
             is_conforming=mesh.is_conforming)
 
 
diff --git a/test/test_array.py b/test/test_array.py
index f01e83ae..a2908a61 100644
--- a/test/test_array.py
+++ b/test/test_array.py
@@ -30,6 +30,7 @@
     dataclass_array_container, pytest_generate_tests_for_array_contexts,
     with_container_arithmetic)
 from pytools.obj_array import make_obj_array
+from pytools.tag import Tag
 
 from meshmode import _acf  # noqa: F401
 from meshmode.array_context import (
@@ -48,7 +49,9 @@
 
 # {{{ test_flatten_unflatten
 
-@with_container_arithmetic(bcast_obj_array=False, rel_comparison=True)
+@with_container_arithmetic(bcast_obj_array=False,
+                           rel_comparison=True,
+                           _cls_has_array_context_attr=True)
 @dataclass_array_container
 @dataclass(frozen=True)
 class MyContainer:
@@ -182,9 +185,6 @@ def test_dof_array_pickling(actx_factory):
     assert actx.to_numpy(flat_norm(dc_of_dofs - dc2_of_dofs, np.inf)) == 0
 
 
-from pytools.tag import Tag
-
-
 class FooTag(Tag):
     pass
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 2257f679..5be44eec 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -562,7 +562,7 @@ def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls,
     # {{{ volume calculation check
 
     if isinstance(mg, SimplexElementGroup):
-        from pytools import factorial
+        from math import factorial
         true_vol = 1/factorial(dim) * 2**dim
     elif isinstance(mg, TensorProductElementGroup):
         true_vol = 2**dim