diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 61ad5852385460f05dbd19f5b4befe464b455839..1d2fbfcfa02885681ebc5e189e25dc1dfe91e694 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -225,6 +225,7 @@ class MeshElementGroup(Record):
     def nunit_nodes(self):
         return self.unit_nodes.shape[-1]
 
+    @property
     @memoize_method
     def is_affine(self):
         return is_affine_group(self)
@@ -1360,19 +1361,36 @@ def is_boundary_tag_empty(mesh, boundary_tag):
 
 # {{{
 
-def is_affine_group(group):
+def is_affine_group(group, abs_tol=None):
+    if abs_tol is None:
+        abs_tol = 1.0e-13
+
+    # get matrices
     basis = mp.simplex_best_available_basis(group.dim, group.order)
     grad_basis = mp.grad_simplex_best_available_basis(group.dim, group.order)
 
-    diff_matrix = mp.differentiation_matrices(basis, grad_basis, group.unit_nodes)
     vinv = la.inv(mp.vandermonde(basis, group.unit_nodes))
+    diff = mp.differentiation_matrices(basis, grad_basis, group.unit_nodes)
+    if not isinstance(diff, tuple):
+        diff = (diff,)
+
+    # construct all second derivative matrices (including cross terms)
+    from itertools import product
+    mats = []
+    for n in product(range(group.dim), repeat=2):
+        if n[0] > n[1]:
+            continue
+        mats.append(vinv.dot(diff[n[0]].dot(diff[n[1]])))
 
-    if not isinstance(diff_matrix, tuple):
-        diff_matrix = (diff_matrix,)
-    mat = tuple(vinv.dot(d.dot(d)) for d in diff_matrix)
+    # check if any element has a non-affine parametrization
+    ddx_coeffs = np.einsum("aij,bj->abi", mats, group.nodes[:, 0, :])
+    norm_inf = np.max(np.abs(ddx_coeffs))
+    if norm_inf > abs_tol:
+        return False
 
-    ddx_coeffs = np.einsum("aij,bj->abi", mat, group.nodes[:, 0, :])
-    return np.max(np.abs(ddx_coeffs)) < 1.0e-13
+    ddx_coeffs = np.einsum("aij,bcj->abci", mats, group.nodes)
+    norm_inf = np.max(np.abs(ddx_coeffs))
+    return norm_inf < abs_tol
 
 # }}}
 
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 92c49b7d0e6183da5375b79a53daab261dc2f540..31221b573c26a7405eed60867e4b8daec6b9b590 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -843,6 +843,7 @@ def generate_warped_rect_mesh(dim, order, n):
     from meshmode.mesh.processing import map_mesh
     return map_mesh(mesh, m)
 
+
 # }}}
 
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 1ae0b63442d9e26d577848e98055ad4c61a853eb..79c462b9ae6fbd434673c93239fa991fc556bfe0 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -1205,8 +1205,25 @@ def test_open_curved_mesh(curve_name):
             closed=closed)
 
 
+def _generate_cross_warped_rect_mesh(dim, order, n):
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(
+            a=(0,)*dim, b=(1,)*dim,
+            n=(n,)*dim, order=order)
+
+    def m(x):
+        results = np.empty_like(x)
+        results[0] = 1 + 1.5 * (x[0] + 0.25) * (x[1] + 0.3)
+        results[1] = x[1]
+        return results
+
+    from meshmode.mesh.processing import map_mesh
+    return map_mesh(mesh, m)
+
+
 @pytest.mark.parametrize("mesh_name", [
-    "box", "warped_box",
+    "box2d", "box3d",
+    "warped_box2d", "warped_box3d", "cross_warped_box",
     "circle", "ellipse",
     "sphere", "torus"
     ])
@@ -1219,16 +1236,20 @@ def test_is_affine_group_check(mesh_name):
     order = 4
     nelements = 32
 
-    if mesh_name == "box":
-        dim = 2
+    if mesh_name.startswith("box"):
+        dim = int(mesh_name[-2])
         is_affine = True
         mesh = generate_regular_rect_mesh(
                 a=(-0.5,)*dim, b=(0.5,)*dim,
                 n=(nelements,)*dim, order=order)
-    elif mesh_name == "warped_box":
+    elif mesh_name.startswith("warped_box"):
+        dim = int(mesh_name[-2])
         is_affine = False
-        dim = 2
         mesh = generate_warped_rect_mesh(dim, order, nelements)
+    elif mesh_name == "cross_warped_box":
+        dim = 2
+        is_affine = False
+        mesh = _generate_cross_warped_rect_mesh(dim, order, nelements)
     elif mesh_name == "circle":
         is_affine = False
         mesh = make_curve_mesh(
@@ -1248,7 +1269,7 @@ def test_is_affine_group_check(mesh_name):
     else:
         raise ValueError("unknown mesh name: {}".format(mesh_name))
 
-    assert all(grp.is_affine() for grp in mesh.groups) == is_affine
+    assert all(grp.is_affine for grp in mesh.groups) == is_affine
 
 
 if __name__ == "__main__":