diff --git a/README.md b/README.md
index ad604d7b..05d8e302 100644
--- a/README.md
+++ b/README.md
@@ -212,6 +212,11 @@ with respect to documented and/or tested features.
 ### Unreleased
 
 - Removed: Python 3.7 support
+- Added: `Mesh.load` supports new keyword arguments
+  `ignore_orientation=True` and `ignore_interior_facets=True` which
+  will both speed up the loading of larger three-dimensional meshes by
+  ignoring facet orientation and all tags not on the boundary,
+  respectively.
 - Fixed: `MeshTet` uniform refine was reindexing subdomains incorrectly
 
 ## [8.1.0] - 2023-06-16
diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
index a058cfeb..b3ee2a78 100644
--- a/skfem/io/meshio.py
+++ b/skfem/io/meshio.py
@@ -50,7 +50,9 @@
 def from_meshio(m,
                 out=None,
                 int_data_to_sets=False,
-                force_meshio_type=None):
+                force_meshio_type=None,
+                ignore_orientation=False,
+                ignore_interior_facets=False):
 
     cells = m.cells_dict
     meshio_type = None
@@ -126,37 +128,49 @@ def from_meshio(m,
         }
         sorted_facets = {k: [tuple(np.sort(f)) for f in v]
                          for k, v in oriented_facets.items()}
-
         for k, v in oriented_facets.items():
-            indices = []
-            oris = []
-            for i, f in enumerate(map(tuple, mtmp.facets.T)):
-                if f in sorted_facets[k]:
-                    indices.append(i)
-                    ix = sorted_facets[k].index(f)
-                    facet = v[ix]
-                    t1, t2 = mtmp.f2t[:, i]
-                    if t2 == -1:
-                        # orientation on boundary is 0
-                        oris.append(0)
-                        continue
-                    if len(f) == 2:
-                        # rotate tangent to find normal
-                        tangent = mtmp.p[:, facet[1]] - mtmp.p[:, facet[0]]
-                        normal = np.array([-tangent[1], tangent[0]])
-                    elif len(f) == 3:
-                        # cross product to find normal
-                        tangent1 = mtmp.p[:, facet[1]] - mtmp.p[:, facet[0]]
-                        tangent2 = mtmp.p[:, facet[2]] - mtmp.p[:, facet[0]]
-                        normal = -np.cross(tangent1, tangent2)
-                    else:
-                        raise NotImplementedError
-                    # find another vector using node outside of boundary
-                    third = np.setdiff1d(mtmp.t[:, t1],
-                                         np.array(f))[0]
-                    outplane = mtmp.p[:, f[0]] - mtmp.p[:, third]
-                    oris.append(1 if np.dot(normal, outplane) > 0 else 0)
-            boundaries[k] = OrientedBoundary(indices, oris)
+            if ignore_orientation or ignore_interior_facets:
+                a = np.array(sorted_facets[k])
+                if ignore_interior_facets:
+                    b = mtmp.facets[:, mtmp.boundary_facets()].T
+                else:
+                    b = mtmp.facets.T
+                boundaries[k] = np.nonzero((a == b[:, None])
+                                           .all(axis=2)
+                                           .any(axis=1))[0]
+            else:
+                indices = []
+                oris = []
+                for i, f in enumerate(map(tuple, mtmp.facets.T)):
+                    if f in sorted_facets[k]:
+                        indices.append(i)
+                        ix = sorted_facets[k].index(f)
+                        facet = v[ix]
+                        t1, t2 = mtmp.f2t[:, i]
+                        if t2 == -1:
+                            # orientation on boundary is 0
+                            oris.append(0)
+                            continue
+                        if len(f) == 2:
+                            # rotate tangent to find normal
+                            tangent = (mtmp.p[:, facet[1]]
+                                       - mtmp.p[:, facet[0]])
+                            normal = np.array([-tangent[1], tangent[0]])
+                        elif len(f) == 3:
+                            # cross product to find normal
+                            tangent1 = (mtmp.p[:, facet[1]]
+                                        - mtmp.p[:, facet[0]])
+                            tangent2 = (mtmp.p[:, facet[2]]
+                                        - mtmp.p[:, facet[0]])
+                            normal = -np.cross(tangent1, tangent2)
+                        else:
+                            raise NotImplementedError
+                        # find another vector using node outside of boundary
+                        third = np.setdiff1d(mtmp.t[:, t1],
+                                             np.array(f))[0]
+                        outplane = mtmp.p[:, f[0]] - mtmp.p[:, third]
+                        oris.append(1 if np.dot(normal, outplane) > 0 else 0)
+                boundaries[k] = OrientedBoundary(indices, oris)
 
     # MSH 2.2 tag parsing
     if len(boundaries) == 0 and m.cell_data and m.field_data:
diff --git a/skfem/visuals/matplotlib.py b/skfem/visuals/matplotlib.py
index 1421a09e..4c7b49d5 100644
--- a/skfem/visuals/matplotlib.py
+++ b/skfem/visuals/matplotlib.py
@@ -1,3 +1,4 @@
+# type: ignore
 """Drawing meshes and solutions using matplotlib."""
 
 from functools import singledispatch
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 2fd3aeda..5f862a9a 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -571,7 +571,21 @@ def test_saveload_cycle_vtk(m):
         MeshTet(),
     ]
 )
-def test_saveload_cycle_tags(fmt, kwargs, m):
+@pytest.mark.parametrize(
+    "ignore_orientation",
+    [
+        True,
+        False,
+    ]
+)
+@pytest.mark.parametrize(
+    "ignore_interior_facets",
+    [
+        True,
+        False,
+    ]
+)
+def test_saveload_cycle_tags(fmt, kwargs, m, ignore_orientation, ignore_interior_facets):
 
     m = (m
          .refined(2)
@@ -582,7 +596,10 @@ def test_saveload_cycle_tags(fmt, kwargs, m):
     with NamedTemporaryFile(suffix=fmt) as f:
         m.save(f.name, point_data={'foo': m.p[0]}, **kwargs)
         out = ['point_data', 'cells_dict']
-        m2 = Mesh.load(f.name, out=out)
+        m2 = Mesh.load(f.name,
+                       out=out,
+                       ignore_orientation=ignore_orientation,
+                       ignore_interior_facets=ignore_interior_facets)
 
 
         assert_array_equal(m.p, m2.p)