diff --git a/README.md b/README.md
index 549245db..3dc05d8f 100644
--- a/README.md
+++ b/README.md
@@ -222,6 +222,11 @@ with respect to documented and/or tested features.
 - Changed: Default tags ('left', 'right', 'top', ...) are no more
   added automatically during mesh initialization, as a workaround you
   can add them explicitly by calling `mesh = mesh.with_defaults()`
+- Changed: All indices within the library are now using `np.int32` for
+  around 10% boost in performance and the corresponding reduction in
+  memory usage for larger meshes - theoretically, the largest possible
+  tetrahedral tensor product mesh is roughly 550 ** 3 = 166 M vertices
+  and 993 M elements, without the facet indexing wrapping over INT_MAX
 
 ### [9.1.1] - 2024-04-23
 
diff --git a/docs/advanced.rst b/docs/advanced.rst
index 8731a040..d6986de3 100644
--- a/docs/advanced.rst
+++ b/docs/advanced.rst
@@ -185,7 +185,7 @@ The DOFs corresponding to the nodes (or vertices) of the mesh are
 .. doctest::
 
    >>> basis.nodal_dofs
-   array([[0, 1, 2, 3, 4, 5, 6, 7]])
+   array([[0, 1, 2, 3, 4, 5, 6, 7]], dtype=int32)
 
 This means that the first (zeroth) entry in the DOF array corresponds to the
 first node/vertex in the finite element mesh (see ``m.p`` for a list of
@@ -207,9 +207,9 @@ edges) and the facets (``m.facets`` for a list of facets) of the mesh are
 .. doctest::
 
    >>> basis.edge_dofs
-   array([[ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]])
+   array([[ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]], dtype=int32)
    >>> basis.facet_dofs
-   array([[20, 21, 22, 23, 24, 25]])
+   array([[20, 21, 22, 23, 24, 25]], dtype=int32)
 
 .. plot::
 
@@ -238,7 +238,7 @@ The remaining DOFs are internal to the element and not shared:
 .. doctest::
 
    >>> basis.interior_dofs
-   array([[26]])
+   array([[26]], dtype=int32)
    
 Each DOF is associated either with a node (``nodal_dofs``), a facet
 (``facet_dofs``), an edge (``edge_dofs``), or an element (``interior_dofs``).
diff --git a/docs/examples/ex41.py b/docs/examples/ex41.py
index 385a4cf8..db255681 100644
--- a/docs/examples/ex41.py
+++ b/docs/examples/ex41.py
@@ -36,7 +36,7 @@
 A = asm(laplace, basis)
 f = asm(unit_load, basis)
 
-y = solve(*condense(A, f, D=out[0]['boundary']['line'].astype(np.int64)))
+y = solve(*condense(A, f, D=out[0]['boundary']['line'].astype(np.int32)))
 
 def visualize():
     from skfem.visuals.matplotlib import plot, draw
diff --git a/docs/howto.rst b/docs/howto.rst
index 3780887e..440561af 100644
--- a/docs/howto.rst
+++ b/docs/howto.rst
@@ -37,9 +37,9 @@ DOFs on the matching facets:
 
    >>> dofs = basis.get_dofs(lambda x: x[0] == 0.)
    >>> dofs.nodal
-   {'u': array([ 0,  2,  5, 10, 14])}
+   {'u': array([ 0,  2,  5, 10, 14], dtype=int32)}
    >>> dofs.facet
-   {'u': array([26, 30, 39, 40])}
+   {'u': array([26, 30, 39, 40], dtype=int32)}
 
 This element has one DOF per node and one DOF per facet.  The facets have their
 own numbering scheme starting from zero, however, the corresponding DOFs are
@@ -48,7 +48,7 @@ offset by the total number of nodal DOFs:
 .. doctest::
 
    >>> dofs.facet['u']
-   array([26, 30, 39, 40])
+   array([26, 30, 39, 40], dtype=int32)
 
 .. plot::
 
@@ -92,9 +92,9 @@ An array of all DOFs with the key ``u`` can be obtained as follows:
 .. doctest::
 
    >>> dofs.all(['u'])
-   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40])
+   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40], dtype=int32)
    >>> dofs.flatten()  # all DOFs, no matter which key
-   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40])
+   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40], dtype=int32)
 
 If a set of facets is tagged, the name of the tag can be passed
 to :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs`:
@@ -103,7 +103,7 @@ to :meth:`~skfem.assembly.basis.AbstractBasis.get_dofs`:
 
    >>> dofs = basis.get_dofs('left')
    >>> dofs.flatten()
-   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40])
+   array([ 0,  2,  5, 10, 14, 26, 30, 39, 40], dtype=int32)
    
 Many DOF types have a well-defined location.  These DOF locations can be found
 as follows:
diff --git a/skfem/assembly/basis/abstract_basis.py b/skfem/assembly/basis/abstract_basis.py
index 579ade03..3be4345f 100644
--- a/skfem/assembly/basis/abstract_basis.py
+++ b/skfem/assembly/basis/abstract_basis.py
@@ -137,7 +137,7 @@ def get_dofs(self,
         >>> m = MeshTri().refined()
         >>> basis = Basis(m, ElementTriP1())
         >>> basis.get_dofs().flatten()
-        array([0, 1, 2, 3, 4, 5, 7, 8])
+        array([0, 1, 2, 3, 4, 5, 7, 8], dtype=int32)
 
         Get DOFs via a function query:
 
@@ -146,7 +146,7 @@ def get_dofs(self,
         >>> m = MeshTri().refined()
         >>> basis = Basis(m, ElementTriP1())
         >>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).flatten()
-        array([0, 2, 5])
+        array([0, 2, 5], dtype=int32)
 
         Get DOFs via named boundaries:
 
@@ -157,7 +157,7 @@ def get_dofs(self,
         ...      .with_boundaries({'left': lambda x: np.isclose(x[0], 0)}))
         >>> basis = Basis(m, ElementTriP1())
         >>> basis.get_dofs('left').flatten()
-        array([0, 2, 5])
+        array([0, 2, 5], dtype=int32)
 
         Get DOFs via named subdomains:
 
@@ -167,7 +167,7 @@ def get_dofs(self,
         ...      .with_subdomains({'left': lambda x: x[0] < .5}))
         >>> basis = Basis(m, ElementTriP1())
         >>> basis.get_dofs(elements='left').flatten()
-        array([0, 2, 4, 5, 6, 8])
+        array([0, 2, 4, 5, 6, 8], dtype=int32)
 
         Further reduce the set of DOFs:
 
@@ -178,7 +178,7 @@ def get_dofs(self,
         >>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).nodal.keys()
         dict_keys(['u', 'u_x', 'u_y', 'u_xx', 'u_xy', 'u_yy'])
         >>> basis.get_dofs(lambda x: np.isclose(x[0], 0)).all(['u', 'u_x'])
-        array([ 0,  1, 12, 13, 30, 31])
+        array([ 0,  1, 12, 13, 30, 31], dtype=int32)
 
         Skip some DOF names altogether:
 
@@ -199,7 +199,7 @@ def get_dofs(self,
         ...                        'right': lambda x: np.isclose(x[0], 1)}))
         >>> basis = Basis(m, ElementTriP1())
         >>> basis.get_dofs({'left', 'right'}).flatten()
-        array([0, 1, 2, 3])
+        array([0, 1, 2, 3], dtype=int32)
 
         Parameters
         ----------
@@ -324,7 +324,7 @@ def split_indices(self) -> List[ndarray]:
         output: List[ndarray] = []
         if isinstance(self.elem, ElementComposite):
             nelems = len(self.elem.elems)
-            o = np.zeros(4, dtype=np.int64)
+            o = np.zeros(4, dtype=np.int32)
             for k in range(nelems):
                 e = self.elem.elems[k]
                 output.append(np.concatenate((
@@ -333,7 +333,7 @@ def split_indices(self) -> List[ndarray]:
                     self.facet_dofs[o[2]:(o[2] + e.facet_dofs)].flatten('F'),
                     (self.interior_dofs[o[3]:(o[3] + e.interior_dofs)]
                      .flatten('F'))
-                )).astype(np.int64))
+                )).astype(np.int32))
                 o += np.array([e.nodal_dofs,
                                e.edge_dofs,
                                e.facet_dofs,
@@ -348,7 +348,7 @@ def split_indices(self) -> List[ndarray]:
                     self.edge_dofs[k::ndims].flatten('F'),
                     self.facet_dofs[k::ndims].flatten('F'),
                     self.interior_dofs[k::ndims].flatten('F'),
-                )).astype(np.int64))
+                )).astype(np.int32))
             return output
         return [np.unique(self.dofs.element_dofs)]
 
diff --git a/skfem/assembly/basis/facet_basis.py b/skfem/assembly/basis/facet_basis.py
index d6776791..43eb1877 100644
--- a/skfem/assembly/basis/facet_basis.py
+++ b/skfem/assembly/basis/facet_basis.py
@@ -75,7 +75,7 @@ def __init__(self,
 
         # by default use boundary facets
         if facets is None:
-            self.find = np.nonzero(self.mesh.f2t[1] == -1)[0]
+            self.find = np.nonzero(self.mesh.f2t[1] == -1)[0].astype(np.int32)
         else:
             self.find = mesh.normalize_facets(facets)
 
diff --git a/skfem/assembly/basis/interior_facet_basis.py b/skfem/assembly/basis/interior_facet_basis.py
index 633481d8..17ae04fc 100644
--- a/skfem/assembly/basis/interior_facet_basis.py
+++ b/skfem/assembly/basis/interior_facet_basis.py
@@ -30,7 +30,7 @@ def __init__(self,
         """Precomputed global basis on interior facets."""
 
         if facets is None:
-            facets = np.nonzero(mesh.f2t[1] != -1)[0]
+            facets = np.nonzero(mesh.f2t[1] != -1)[0].astype(np.int32)
 
         facets = mesh.normalize_facets(facets)
 
diff --git a/skfem/assembly/dofs.py b/skfem/assembly/dofs.py
index 241f56fa..153a6acb 100644
--- a/skfem/assembly/dofs.py
+++ b/skfem/assembly/dofs.py
@@ -260,7 +260,7 @@ def __init__(self, topo, element, offset=0):
         self.element = element
 
         self.nodal_dofs = np.reshape(
-            np.arange(element.nodal_dofs * topo.nvertices, dtype=np.int64),
+            np.arange(element.nodal_dofs * topo.nvertices, dtype=np.int32),
             (element.nodal_dofs, topo.nvertices),
             order='F') + offset
         offset += element.nodal_dofs * topo.nvertices
@@ -269,32 +269,32 @@ def __init__(self, topo, element, offset=0):
         if element.dim == 3 and element.edge_dofs > 0:
             self.edge_dofs = np.reshape(
                 np.arange(element.edge_dofs * topo.nedges,
-                          dtype=np.int64),
+                          dtype=np.int32),
                 (element.edge_dofs, topo.nedges),
                 order='F') + offset
             offset += element.edge_dofs * topo.nedges
         else:
-            self.edge_dofs = np.empty((0, 0), dtype=np.int64)
+            self.edge_dofs = np.empty((0, 0), dtype=np.int32)
 
         # facet dofs
         if element.facet_dofs > 0:
             self.facet_dofs = np.reshape(
                 np.arange(element.facet_dofs * topo.nfacets,
-                          dtype=np.int64),
+                          dtype=np.int32),
                 (element.facet_dofs, topo.nfacets),
                 order='F') + offset
             offset += element.facet_dofs * topo.nfacets
         else:
-            self.facet_dofs = np.empty((0, 0), dtype=np.int64)
+            self.facet_dofs = np.empty((0, 0), dtype=np.int32)
 
         # interior dofs
         self.interior_dofs = np.reshape(
-            np.arange(element.interior_dofs * topo.nelements, dtype=np.int64),
+            np.arange(element.interior_dofs * topo.nelements, dtype=np.int32),
             (element.interior_dofs, topo.nelements),
             order='F') + offset
 
         # global numbering
-        self.element_dofs = np.zeros((0, topo.nelements), dtype=np.int64)
+        self.element_dofs = np.zeros((0, topo.nelements), dtype=np.int32)
 
         # nodal dofs
         for itr in range(topo.t.shape[0]):
@@ -350,9 +350,9 @@ def get_vertex_dofs(
         return DofsView(
             self,
             nodes,
-            np.empty((0,), dtype=np.int64),
-            np.empty((0,), dtype=np.int64),
-            np.empty((0,), dtype=np.int64),
+            np.empty((0,), dtype=np.int32),
+            np.empty((0,), dtype=np.int32),
+            np.empty((0,), dtype=np.int32),
             r1,
             r2,
             r3,
@@ -376,13 +376,13 @@ def get_element_dofs(
             An array of dofnames to skip.
 
         """
-        nodal_ix = (np.empty((0,), dtype=np.int64)
+        nodal_ix = (np.empty((0,), dtype=np.int32)
                     if self.element.nodal_dofs == 0
                     else np.unique(self.topo.t[:, elements]))
-        edge_ix = (np.empty((0,), dtype=np.int64)
+        edge_ix = (np.empty((0,), dtype=np.int32)
                    if self.element.edge_dofs == 0
                    else np.unique(self.topo.t2e[:, elements]))
-        facet_ix = (np.empty((0,), dtype=np.int64)
+        facet_ix = (np.empty((0,), dtype=np.int32)
                     if self.element.facet_dofs == 0
                     else np.unique(self.topo.t2f[:, elements]))
         interior_ix = elements
@@ -424,13 +424,13 @@ def get_facet_dofs(
         if self.element.nodal_dofs > 0 or self.element.edge_dofs > 0:
             nodal_ix, edge_ix = self.topo._expand_facets(facets)
 
-        nodal_ix = (np.empty((0,), dtype=np.int64)
+        nodal_ix = (np.empty((0,), dtype=np.int32)
                     if self.element.nodal_dofs == 0
                     else nodal_ix)
-        edge_ix = (np.empty((0,), dtype=np.int64)
+        edge_ix = (np.empty((0,), dtype=np.int32)
                    if self.element.edge_dofs == 0
                    else edge_ix)
-        facet_ix = (np.empty((0,), dtype=np.int64)
+        facet_ix = (np.empty((0,), dtype=np.int32)
                     if self.element.facet_dofs == 0
                     else facets)
 
@@ -444,7 +444,7 @@ def get_facet_dofs(
             nodal_ix,
             facet_ix,
             edge_ix,
-            np.empty((0,), dtype=np.int64),
+            np.empty((0,), dtype=np.int32),
             r1,
             r2,
             r3,
@@ -466,7 +466,7 @@ def _by_name(self,
 
         ents = {
             self.element.dofnames[rows[i] + off]: np.zeros((0, n_ents),
-                                                           dtype=np.int64)
+                                                           dtype=np.int32)
             for i in range(n_dofs)
         }
         for i in range(n_dofs):
diff --git a/skfem/assembly/form/bilinear_form.py b/skfem/assembly/form/bilinear_form.py
index ab7cc4f7..0324c9dd 100644
--- a/skfem/assembly/form/bilinear_form.py
+++ b/skfem/assembly/form/bilinear_form.py
@@ -79,8 +79,8 @@ def _assemble(self,
         # initialize COO data structures
         sz = ubasis.Nbfun * vbasis.Nbfun * nt
         data = np.zeros((ubasis.Nbfun, vbasis.Nbfun, nt), dtype=self.dtype)
-        rows = np.zeros(sz, dtype=np.int64)
-        cols = np.zeros(sz, dtype=np.int64)
+        rows = np.zeros(sz, dtype=np.int32)
+        cols = np.zeros(sz, dtype=np.int32)
 
         # loop over the indices of local stiffness matrix
         for j in range(ubasis.Nbfun):
diff --git a/skfem/assembly/form/linear_form.py b/skfem/assembly/form/linear_form.py
index cd7a73da..67b807da 100644
--- a/skfem/assembly/form/linear_form.py
+++ b/skfem/assembly/form/linear_form.py
@@ -36,7 +36,7 @@ def _assemble(self,
         # initialize COO data structures
         sz = vbasis.Nbfun * nt
         data = np.zeros(sz, dtype=self.dtype)
-        rows = np.zeros(sz, dtype=np.int64)
+        rows = np.zeros(sz, dtype=np.int32)
 
         for i in range(vbasis.Nbfun):
             ixs = slice(nt * i, nt * (i + 1))
diff --git a/skfem/assembly/form/trilinear_form.py b/skfem/assembly/form/trilinear_form.py
index 2dba2f8d..de2dc74f 100644
--- a/skfem/assembly/form/trilinear_form.py
+++ b/skfem/assembly/form/trilinear_form.py
@@ -33,9 +33,9 @@ def _assemble(self,
         # initialize COO data structures
         sz = (ubasis.Nbfun, vbasis.Nbfun, wbasis.Nbfun, nt)
         data = np.zeros(sz, dtype=self.dtype)
-        rows = np.zeros(sz, dtype=np.int64)
-        cols = np.zeros(sz, dtype=np.int64)
-        mats = np.zeros(sz, dtype=np.int64)
+        rows = np.zeros(sz, dtype=np.int32)
+        cols = np.zeros(sz, dtype=np.int32)
+        mats = np.zeros(sz, dtype=np.int32)
 
         # loop over the indices of local stiffness matrix
         for k in range(ubasis.Nbfun):
diff --git a/skfem/element/element_hcurl.py b/skfem/element/element_hcurl.py
index 4642ab0f..a8ade43f 100644
--- a/skfem/element/element_hcurl.py
+++ b/skfem/element/element_hcurl.py
@@ -16,7 +16,7 @@ def orient(self, mapping, i, tind=None):
         ix = int(i / divide_by)
         if mapping.mesh.dim() == 2 and ix >= self.refdom.nfacets:
             # no orientation required for interior DOFs => return 1
-            ori = np.ones(mapping.mesh.t.shape[1], dtype=np.int64)
+            ori = np.ones(mapping.mesh.t.shape[1], dtype=np.int32)
             return ori[tind]
         if mapping.mesh.dim() == 3:
             t1, t2 = mapping.mesh.refdom.edges[ix]
diff --git a/skfem/element/element_hdiv.py b/skfem/element/element_hdiv.py
index 11f2f947..340be373 100644
--- a/skfem/element/element_hdiv.py
+++ b/skfem/element/element_hdiv.py
@@ -16,7 +16,7 @@ def orient(self, mapping, i, tind=None):
             # TODO support edge DOFs
             # TODO can you just skip np.arange here? len(tind)?
             return np.ones(len(np.arange(mapping.mesh.t.shape[1])[tind]),
-                           dtype=np.int64)
+                           dtype=np.int32)
         ori = -1 + 2 * (mapping.mesh.f2t[0, mapping.mesh.t2f[ix]]
                         == np.arange(mapping.mesh.t.shape[1]))
         return ori[tind]
diff --git a/skfem/experimental/autodiff/__init__.py b/skfem/experimental/autodiff/__init__.py
index a0cd3f89..d4bc8c7e 100644
--- a/skfem/experimental/autodiff/__init__.py
+++ b/skfem/experimental/autodiff/__init__.py
@@ -131,12 +131,12 @@ def assemble(self, basis, x=None, **kwargs):
         # initialize COO data structures
         sz = basis.Nbfun * basis.Nbfun * nt
         data = np.zeros((basis.Nbfun, basis.Nbfun, nt), dtype=self.dtype)
-        rows = np.zeros(sz, dtype=np.int64)
-        cols = np.zeros(sz, dtype=np.int64)
+        rows = np.zeros(sz, dtype=np.int32)
+        cols = np.zeros(sz, dtype=np.int32)
 
         sz1 = basis.Nbfun * nt
         data1 = np.zeros(sz1, dtype=self.dtype)
-        rows1 = np.zeros(sz1, dtype=np.int64)
+        rows1 = np.zeros(sz1, dtype=np.int32)
 
         # # autograd version
         # def _make_jacobian(V):
diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
index 930d2aa4..f9913531 100644
--- a/skfem/io/meshio.py
+++ b/skfem/io/meshio.py
@@ -118,7 +118,7 @@ def from_meshio(m,
 
     # parse any subdomains from cell_sets
     if m.cell_sets:
-        subdomains = {k: v[meshio_type].astype(np.int64)
+        subdomains = {k: v[meshio_type].astype(np.int32)
                       for k, v in m.cell_sets_dict.items()
                       if meshio_type in v}
 
@@ -136,7 +136,7 @@ def from_meshio(m,
                 ind = p2f[:, sorted_facets[0]]
                 for itr in range(sorted_facets.shape[0] - 1):
                     ind = ind.multiply(p2f[:, sorted_facets[itr + 1]])
-                boundaries[k] = np.nonzero(ind)[0]
+                boundaries[k] = np.nonzero(ind)[0].astype(np.int32)
 
                 if not ignore_orientation:
                     try:
@@ -179,7 +179,7 @@ def find_tagname(tag):
                 return None
 
             for tag in tags:
-                t_set = np.nonzero(tag == elements_tag)[0]
+                t_set = np.nonzero(tag == elements_tag)[0].astype(np.int32)
                 subdomains[find_tagname(tag)] = t_set
 
             # find tagged boundaries
@@ -201,7 +201,7 @@ def find_tagname(tag):
             tags = index[:, 0]
             boundaries = {}
             for tag in np.unique(tags):
-                tagindex = np.nonzero(tags == tag)[0]
+                tagindex = np.nonzero(tags == tag)[0].astype(np.int32)
                 boundaries[find_tagname(tag)] = index[tagindex, 1]
 
         except Exception:
diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py
index 41255983..150b3b1a 100644
--- a/skfem/mapping/mapping_affine.py
+++ b/skfem/mapping/mapping_affine.py
@@ -272,7 +272,7 @@ def normals(self, X, tind, find, t2f):
         N = np.empty((self.dim, len(find)))
 
         for itr in range(Nref.shape[0]):
-            ix = np.nonzero(t2f[itr, tind] == find)[0]
+            ix = np.nonzero(t2f[itr, tind] == find)[0].astype(np.int32)
             for jtr in range(Nref.shape[1]):
                 N[jtr, ix] = Nref[itr, jtr]
 
diff --git a/skfem/mapping/mapping_isoparametric.py b/skfem/mapping/mapping_isoparametric.py
index 6efc8eaf..d2c7ad03 100644
--- a/skfem/mapping/mapping_isoparametric.py
+++ b/skfem/mapping/mapping_isoparametric.py
@@ -232,7 +232,7 @@ def normals(self, X, tind, find, t2f):
         N = np.empty((self.dim, len(find)))
 
         for itr in range(Nref.shape[0]):
-            ix = np.nonzero(t2f[itr, tind] == find)[0]
+            ix = np.nonzero(t2f[itr, tind] == find)[0].astype(np.int32)
             for jtr in range(Nref.shape[1]):
                 N[jtr, ix] = Nref[itr, jtr]
 
diff --git a/skfem/mesh/__init__.py b/skfem/mesh/__init__.py
index dd6ec989..6685e6df 100644
--- a/skfem/mesh/__init__.py
+++ b/skfem/mesh/__init__.py
@@ -58,7 +58,7 @@ def __call__(self, p=None, t=None, **kwargs):
             p = np.atleast_2d(p)
 
         if p is not None and t is None:
-            tmp = np.arange(p.shape[1] - 1, dtype=np.int64)
+            tmp = np.arange(p.shape[1] - 1, dtype=np.int32)
             t = np.vstack((tmp, tmp + 1))
             return MeshLine1(p, t, **kwargs)
 
diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py
index 17293974..15a6fc23 100644
--- a/skfem/mesh/mesh.py
+++ b/skfem/mesh/mesh.py
@@ -196,7 +196,7 @@ def dim(self):
 
     def boundary_facets(self) -> ndarray:
         """Return an array of boundary facet indices."""
-        return np.nonzero(self.f2t[1] == -1)[0]
+        return np.nonzero(self.f2t[1] == -1)[0].astype(np.int32)
 
     def boundary_edges(self) -> ndarray:
         """Return an array of boundary edge indices."""
@@ -355,7 +355,7 @@ def encode_boundary(boundary: Union[ndarray, OrientedBoundary]
             **{
                 f"skfem:s:{name}": [
                     np.isin(np.arange(self.t.shape[1],
-                                      dtype=np.int64), subdomain).astype(int)
+                                      dtype=np.int32), subdomain).astype(int)
                 ]
                 for name, subdomain in subdomains.items()
             },
@@ -379,7 +379,7 @@ def _decode_cell_data(self, cell_data: Dict[str, List[ndarray]]):
             elif subnames[1] == "b":
                 mask = (
                     (1 << np.arange(self.refdom.nfacets))[:, None]
-                    & data[0].astype(np.int64)
+                    & data[0].astype(np.int32)
                 ).astype(bool)
                 facets = np.sort(self.t2f[mask])
                 cells = mask.nonzero()[1]
@@ -413,7 +413,7 @@ def nodes_satisfying(self,
             If ``True``, include only boundary facets.
 
         """
-        nodes = np.nonzero(test(self.p))[0]
+        nodes = np.nonzero(test(self.p))[0].astype(np.int32)
         if boundaries_only:
             nodes = np.intersect1d(nodes, self.boundary_nodes())
         return nodes
@@ -436,7 +436,7 @@ def facets_satisfying(self,
 
         """
         midp = self.p[:, self.facets].mean(axis=1)
-        facets = np.nonzero(test(midp))[0]
+        facets = np.nonzero(test(midp))[0].astype(np.int32)
         if boundaries_only:
             facets = np.intersect1d(facets, self.boundary_facets())
         if normal is not None:
@@ -466,9 +466,11 @@ def facets_around(self, elements, flip=False) -> OrientedBoundary:
         facets, counts = np.unique(self.t2f[:, elements], return_counts=True)
         facets = facets[counts == 1]
         if flip:
-            ori = np.nonzero(~np.isin(self.f2t[:, facets], elements).T)[1]
+            ori = (np.nonzero(~np.isin(self.f2t[:, facets], elements).T)[1]
+                   .astype(np.int32))
         else:
-            ori = np.nonzero(np.isin(self.f2t[:, facets], elements).T)[1]
+            ori = (np.nonzero(np.isin(self.f2t[:, facets], elements).T)[1]
+                   .astype(np.int32))
         return OrientedBoundary(facets, ori)
 
     def elements_satisfying(self,
@@ -483,7 +485,7 @@ def elements_satisfying(self,
 
         """
         midp = self.p[:, self.t].mean(axis=1)
-        return np.nonzero(test(midp))[0].astype(np.int64)
+        return np.nonzero(test(midp))[0].astype(np.int32)
 
     def _expand_facets(self, ix: ndarray) -> Tuple[ndarray, ndarray]:
         """Return vertices and edges corresponding to given facet indices.
@@ -499,7 +501,7 @@ def _expand_facets(self, ix: ndarray) -> Tuple[ndarray, ndarray]:
         if self.dim() == 3 and self.bndelem is not None:
             edges = np.unique(self.f2e[:, ix])
         else:
-            edges = np.array([], dtype=np.int64)
+            edges = np.array([], dtype=np.int32)
 
         return vertices, edges
 
@@ -549,7 +551,7 @@ def __post_init__(self):
             self.t = np.sort(self.t, axis=0)
 
         self.doflocs = np.asarray(self.doflocs, dtype=np.float64, order="K")
-        self.t = np.asarray(self.t, dtype=np.int64, order="K")
+        self.t = np.asarray(self.t, dtype=np.int32, order="K")
 
         M = self.elem.refdom.nnodes
 
@@ -560,7 +562,7 @@ def __post_init__(self):
             p, t = self.doflocs, self.t
             t_nodes = t[:M]
             uniq, ix = np.unique(t_nodes, return_inverse=True)
-            self.t = (np.arange(len(uniq), dtype=np.int64)[ix]
+            self.t = (np.arange(len(uniq), dtype=np.int32)[ix]
                       .reshape(t_nodes.shape))
             doflocs = np.hstack((
                 p[:, uniq],
@@ -583,6 +585,18 @@ def __post_init__(self):
                                "to C_CONTIGUOUS.")
             self.t = np.ascontiguousarray(self.t)
 
+        # normalize data types
+        # if self._boundaries is not None:
+        #     self._boundaries = {
+        #         k: v.astype(np.int32)
+        #         for k, v in self._boundaries.items()
+        #     }
+        # if self._subdomains is not None:
+        #     self._subdomains = {
+        #         k: v.astype(np.int32)
+        #         for k, v in self._subdomains.items()
+        #     }
+
         # run validation
         if self.validate and logger.getEffectiveLevel() <= logging.DEBUG:
             self.is_valid()
@@ -867,8 +881,8 @@ def refined(self, times_or_ix: Union[int, ndarray] = 1):
                 # fix subdomains for remaining mesh types
                 if m._subdomains is not None and mtmp._subdomains is None:
                     N = int(mtmp.t.shape[1] / m.t.shape[1])
-                    new_t = np.zeros((N, m.t.shape[1]), dtype=np.int64)
-                    new_t[0] = np.arange(m.t.shape[1], dtype=np.int64)
+                    new_t = np.zeros((N, m.t.shape[1]), dtype=np.int32)
+                    new_t[0] = np.arange(m.t.shape[1], dtype=np.int32)
                     for itr in range(N - 1):
                         new_t[itr + 1] = new_t[itr] + m.t.shape[1]
                     mtmp = replace(
@@ -1069,7 +1083,7 @@ def build_inverse(t, mapping):
         e_last, ix_last = np.unique(e[::-1], return_index=True)
         ix_last = e.shape[0] - ix_last - 1
 
-        inverse = np.zeros((2, np.max(mapping) + 1), dtype=np.int64)
+        inverse = np.zeros((2, np.max(mapping) + 1), dtype=np.int32)
         inverse[0, e_first] = tix[ix_first]
         inverse[1, e_last] = tix[ix_last]
 
@@ -1089,8 +1103,8 @@ def param(self) -> float:
     def _reix(self, ix: ndarray) -> Tuple[ndarray, ndarray, ndarray]:
         """Connect ``self.p`` based on the indices ``ix``."""
         ixuniq = np.unique(ix)
-        t = np.zeros(np.max(ix) + 1, dtype=np.int64)
-        t[ixuniq] = np.arange(len(ixuniq), dtype=np.int64)
+        t = np.zeros(np.max(ix) + 1, dtype=np.int32)
+        t[ixuniq] = np.arange(len(ixuniq), dtype=np.int32)
         return (
             np.ascontiguousarray(self.p[:, ixuniq]),
             np.ascontiguousarray(t[ix]),
@@ -1147,21 +1161,21 @@ def restrict(self,
         new_subdomains = None
         if not skip_subdomains and self.subdomains is not None:
             # map from old to new element index
-            newt = np.zeros(self.t.shape[1], dtype=np.int64) - 1
-            newt[elements] = np.arange(len(elements), dtype=np.int64)
+            newt = np.zeros(self.t.shape[1], dtype=np.int32) - 1
+            newt[elements] = np.arange(len(elements), dtype=np.int32)
             # remove 'elements' from each subdomain and remap
             new_subdomains = {
                 k: newt[np.intersect1d(self.subdomains[k],
-                                       elements).astype(np.int64)]
+                                       elements).astype(np.int32)]
                 for k in self.subdomains
             }
 
         new_boundaries = None
         if not skip_boundaries and self.boundaries is not None:
             # map from old to new facet index
-            newf = np.zeros(self.facets.shape[1], dtype=np.int64) - 1
+            newf = np.zeros(self.facets.shape[1], dtype=np.int32) - 1
             facets = np.unique(self.t2f[:, elements])
-            newf[facets] = np.arange(len(facets), dtype=np.int64)
+            newf[facets] = np.arange(len(facets), dtype=np.int32)
             new_boundaries = {k: newf[self.boundaries[k]]
                               for k in self.boundaries}
             # filter facets not existing in the new mesh, value is -1
@@ -1192,7 +1206,7 @@ def remove_elements(self, elements):
         """
         elements = self.normalize_elements(elements)
         return self.restrict(np.setdiff1d(np.arange(self.t.shape[1],
-                                                    dtype=np.int64),
+                                                    dtype=np.int32),
                                           elements))
 
     def remove_unused_nodes(self):
@@ -1317,7 +1331,7 @@ def normalize_elements(self, elements) -> ndarray:
 
         """
         if isinstance(elements, bool) and elements:
-            return np.arange(self.nelements, dtype=np.int64)
+            return np.arange(self.nelements, dtype=np.int32)
         if isinstance(elements, int):
             # Make  normalize_elements([1,2,3]) have the same behavior as
             # normalize_elements(np.array([1,2,3]))
diff --git a/skfem/mesh/mesh_3d.py b/skfem/mesh/mesh_3d.py
index e8a64302..0fd9afcf 100644
--- a/skfem/mesh/mesh_3d.py
+++ b/skfem/mesh/mesh_3d.py
@@ -29,7 +29,8 @@ def edges_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray:
             belonging to the output set.
 
         """
-        return np.nonzero(test(self.p[:, self.edges].mean(1)))[0]
+        return (np.nonzero(test(self.p[:, self.edges].mean(1)))[0]
+                .astype(np.int32))
 
     def boundary_edges(self) -> ndarray:
         """Return an array of boundary edge indices."""
@@ -51,5 +52,5 @@ def boundary_edges(self) -> ndarray:
 
     def interior_edges(self) -> ndarray:
         """Return an array of interior edge indices."""
-        return np.setdiff1d(np.arange(self.edges.shape[1], dtype=np.int64),
+        return np.setdiff1d(np.arange(self.edges.shape[1], dtype=np.int32),
                             self.boundary_edges())
diff --git a/skfem/mesh/mesh_dg.py b/skfem/mesh/mesh_dg.py
index d6f9c401..bafdee3f 100644
--- a/skfem/mesh/mesh_dg.py
+++ b/skfem/mesh/mesh_dg.py
@@ -10,8 +10,8 @@ def init_tensor(cls,
                     *args,
                     periodic=[]):
         mesh = cls.__bases__[-1].init_tensor(*args)
-        ix = np.empty((0,), dtype=np.int64)
-        ix0 = np.empty((0,), dtype=np.int64)
+        ix = np.empty((0,), dtype=np.int32)
+        ix0 = np.empty((0,), dtype=np.int32)
         for dim in periodic:
             argmin = args[dim].min()
             argmax = args[dim].max()
@@ -59,12 +59,12 @@ def periodic(cls, mesh, ix, ix0):
                              "creation of a periodic mesh should be equal.")
 
         # reorder vertices: eliminated nodes must have highest index values
-        remap = np.empty(mesh.nvertices, dtype=np.int64)
+        remap = np.empty(mesh.nvertices, dtype=np.int32)
         remap[ix] = np.arange(mesh.nvertices - len(ix),
                               mesh.nvertices,
-                              dtype=np.int64)
-        oix = np.setdiff1d(np.arange(mesh.nvertices, dtype=np.int64), ix)
-        remap[oix] = np.arange(mesh.nvertices - len(ix), dtype=np.int64)
+                              dtype=np.int32)
+        oix = np.setdiff1d(np.arange(mesh.nvertices, dtype=np.int32), ix)
+        remap[oix] = np.arange(mesh.nvertices - len(ix), dtype=np.int32)
 
         doflocs = np.hstack((mesh.doflocs[:, oix], mesh.doflocs[:, ix]))
         t = remap[mesh.t]
@@ -77,7 +77,7 @@ def periodic(cls, mesh, ix, ix0):
         )
 
         # make periodic
-        reremap = np.arange(mesh.nvertices, dtype=np.int64)
+        reremap = np.arange(mesh.nvertices, dtype=np.int32)
         ix1 = remap[ix]
         ix2 = remap[ix0]
         reremap[ix1] = ix2
diff --git a/skfem/mesh/mesh_hex_1.py b/skfem/mesh/mesh_hex_1.py
index f76719a7..945ca605 100644
--- a/skfem/mesh/mesh_hex_1.py
+++ b/skfem/mesh/mesh_hex_1.py
@@ -41,7 +41,7 @@ class MeshHex1(Mesh3D):
     )
     t: ndarray = field(
         default_factory=lambda: np.array(
-            [[0, 1, 2, 3, 4, 5, 6, 7]], dtype=np.int64
+            [[0, 1, 2, 3, 4, 5, 6, 7]], dtype=np.int32
         ).T
     )
     elem: Type[Element] = ElementHex1
@@ -60,7 +60,7 @@ def _uniform(self):
         sz = p.shape[1]
         t2e = self.t2e.copy() + sz
         t2f = self.t2f.copy() + np.max(t2e) + 1
-        mid = np.arange(self.t.shape[1], dtype=np.int64) + np.max(t2f) + 1
+        mid = np.arange(self.t.shape[1], dtype=np.int32) + np.max(t2f) + 1
 
         doflocs = np.hstack((
             p,
@@ -152,7 +152,7 @@ def init_tensor(cls: Type,
                 .reshape(ne, 1, order='F')
                 .copy()
                 .flatten())
-        return cls(p, t.astype(np.int64))
+        return cls(p, t.astype(np.int32))
 
     def to_meshtet(self):
         """Split each hexahedron into six tetrahedra."""
diff --git a/skfem/mesh/mesh_line_1.py b/skfem/mesh/mesh_line_1.py
index 2ebffbb9..9b5aedac 100644
--- a/skfem/mesh/mesh_line_1.py
+++ b/skfem/mesh/mesh_line_1.py
@@ -18,7 +18,7 @@ class MeshLine1(MeshSimplex, Mesh):
     doflocs: ndarray = field(
         default_factory=lambda: np.array([[0., 1.]], dtype=np.float64))
     t: ndarray = field(default_factory=lambda: np.array(
-        [[0], [1]], dtype=np.int64))
+        [[0], [1]], dtype=np.int32))
     elem: Type[Element] = ElementLineP1
     affine: bool = True
 
@@ -84,7 +84,7 @@ def finder(x):
             xin = x.copy()  # bring endpoint inside for np.digitize
             xin[x == self.p[0, ix[-1]]] = self.p[0, ix[-2:]].mean()
             elems = np.nonzero(ix[np.digitize(xin, self.p[0, ix])][:, None]
-                               == maxt)[1]
+                               == maxt)[1].astype(np.int32)
             if len(elems) < len(x):
                 raise ValueError("Point is outside of the mesh.")
             return elems
diff --git a/skfem/mesh/mesh_quad_1.py b/skfem/mesh/mesh_quad_1.py
index 4526aafd..de1b1d9a 100644
--- a/skfem/mesh/mesh_quad_1.py
+++ b/skfem/mesh/mesh_quad_1.py
@@ -29,7 +29,7 @@ class MeshQuad1(Mesh2D):
         ).T
     )
     t: ndarray = field(
-        default_factory=lambda: np.array([[0, 1, 2, 3]], dtype=np.int64).T
+        default_factory=lambda: np.array([[0, 1, 2, 3]], dtype=np.int32).T
     )
     elem: Type[Element] = ElementQuad1
 
@@ -39,7 +39,7 @@ def _uniform(self):
         t = self.t
         sz = p.shape[1]
         t2f = self.t2f.copy()
-        mid = np.arange(t.shape[1], dtype=np.int64) + np.max(t2f) + sz + 1
+        mid = np.arange(t.shape[1], dtype=np.int32) + np.max(t2f) + sz + 1
 
         m = replace(
             self,
@@ -60,8 +60,8 @@ def _uniform(self):
 
         if self._boundaries is not None:
             # mapping of indices between old and new facets
-            new_facets = np.zeros((2, self.facets.shape[1]), dtype=np.int64)
-            ix0 = np.arange(t.shape[1], dtype=np.int64)
+            new_facets = np.zeros((2, self.facets.shape[1]), dtype=np.int32)
+            ix0 = np.arange(t.shape[1], dtype=np.int32)
             ix1 = ix0 + t.shape[1]
             ix2 = ix0 + 2 * t.shape[1]
             ix3 = ix0 + 3 * t.shape[1]
@@ -130,7 +130,7 @@ def init_tensor(cls: Type,
         t[3] = (ix[0:(npy-1), 1:npx].reshape(nt, 1, order='F')
                 .copy()
                 .flatten())
-        return cls(p, t.astype(np.int64))
+        return cls(p, t.astype(np.int32))
 
     def to_meshtri(self,
                    x: Optional[ndarray] = None,
@@ -150,7 +150,7 @@ def to_meshtri(self,
         if style == 'x':
             tnew = np.arange(np.max(self.t) + 1,
                              np.max(self.t) + 1 + self.t.shape[1],
-                             dtype=np.int64)
+                             dtype=np.int32)
             t = np.hstack((
                 np.vstack((self.t[[0, 1]], tnew)),
                 np.vstack((self.t[[1, 2]], tnew)),
diff --git a/skfem/mesh/mesh_simplex.py b/skfem/mesh/mesh_simplex.py
index 5f4e71a8..ba09d19a 100644
--- a/skfem/mesh/mesh_simplex.py
+++ b/skfem/mesh/mesh_simplex.py
@@ -11,7 +11,7 @@ def orientation(self):
         mapping = self._mapping()
         return (np.sign(mapping.detDF(np.zeros(self.p.shape[0])[:, None]))
                 .flatten()
-                .astype(np.int64))
+                .astype(np.int32))
 
     def oriented(self):
         """Return a oriented mesh with positive Jacobian determinant.
@@ -19,7 +19,7 @@ def oriented(self):
         For triangular meshes this corresponds to CCW orientation.
 
         """
-        flip = np.nonzero(self.orientation() == -1)[0]
+        flip = np.nonzero(self.orientation() == -1)[0].astype(np.int32)
         t = self.t.copy()
         t0 = t[0, flip]
         t1 = t[1, flip]
diff --git a/skfem/mesh/mesh_tet_1.py b/skfem/mesh/mesh_tet_1.py
index a220b845..3117c217 100644
--- a/skfem/mesh/mesh_tet_1.py
+++ b/skfem/mesh/mesh_tet_1.py
@@ -37,7 +37,7 @@ class MeshTet1(MeshSimplex, Mesh3D):
                 [2, 3, 1, 7],
                 [1, 2, 4, 7],
             ],
-            dtype=np.int64,
+            dtype=np.int32,
         ).T
     )
     elem: Type[Element] = ElementTetP1
@@ -63,7 +63,7 @@ def finder(x, y, z, _search_all=False):
                 _, ix_ind = np.unique(ix, return_index=True)
                 ix = ix[np.sort(ix_ind)]
             else:
-                ix = np.arange(nelems, dtype=np.int64)
+                ix = np.arange(nelems, dtype=np.int32)
 
             X = mapping.invF(np.array([x, y, z])[:, None], ix)
             eps = np.finfo(X.dtype).eps
@@ -127,8 +127,8 @@ def _uniform(self):
 
         if self._subdomains is not None:
             nt = t.shape[1]
-            new_t = np.zeros((8, nt), dtype=np.int64)
-            new_t[0] = np.arange(nt, dtype=np.int64)
+            new_t = np.zeros((8, nt), dtype=np.int32)
+            new_t[0] = np.arange(nt, dtype=np.int32)
             new_t[1] = new_t[0] + nt
             new_t[2] = new_t[1] + nt
             new_t[3] = new_t[2] + nt
@@ -136,20 +136,20 @@ def _uniform(self):
             n2 = np.sum(c2)
             n3 = np.sum(c3)
             if n1 > 0:
-                new_t[4, c1] = np.arange(n1, dtype=np.int64) + 4 * nt
-                new_t[5, c1] = np.arange(n1, dtype=np.int64) + 5 * nt
-                new_t[6, c1] = np.arange(n1, dtype=np.int64) + 6 * nt
-                new_t[7, c1] = np.arange(n1, dtype=np.int64) + 7 * nt
+                new_t[4, c1] = np.arange(n1, dtype=np.int32) + 4 * nt
+                new_t[5, c1] = np.arange(n1, dtype=np.int32) + 5 * nt
+                new_t[6, c1] = np.arange(n1, dtype=np.int32) + 6 * nt
+                new_t[7, c1] = np.arange(n1, dtype=np.int32) + 7 * nt
             if n2 > 0:
-                new_t[4, c2] = np.arange(n2, dtype=np.int64) + 4 * nt + n1
-                new_t[5, c2] = np.arange(n2, dtype=np.int64) + 5 * nt + n1
-                new_t[6, c2] = np.arange(n2, dtype=np.int64) + 6 * nt + n1
-                new_t[7, c2] = np.arange(n2, dtype=np.int64) + 7 * nt + n1
+                new_t[4, c2] = np.arange(n2, dtype=np.int32) + 4 * nt + n1
+                new_t[5, c2] = np.arange(n2, dtype=np.int32) + 5 * nt + n1
+                new_t[6, c2] = np.arange(n2, dtype=np.int32) + 6 * nt + n1
+                new_t[7, c2] = np.arange(n2, dtype=np.int32) + 7 * nt + n1
             if n3 > 0:
-                new_t[4, c3] = np.arange(n3, dtype=np.int64) + 4 * nt + n1 + n2
-                new_t[5, c3] = np.arange(n3, dtype=np.int64) + 5 * nt + n1 + n2
-                new_t[6, c3] = np.arange(n3, dtype=np.int64) + 6 * nt + n1 + n2
-                new_t[7, c3] = np.arange(n3, dtype=np.int64) + 7 * nt + n1 + n2
+                new_t[4, c3] = np.arange(n3, dtype=np.int32) + 4 * nt + n1 + n2
+                new_t[5, c3] = np.arange(n3, dtype=np.int32) + 5 * nt + n1 + n2
+                new_t[6, c3] = np.arange(n3, dtype=np.int32) + 6 * nt + n1 + n2
+                new_t[7, c3] = np.arange(n3, dtype=np.int32) + 7 * nt + n1 + n2
             subdomains = {
                 name: np.sort(new_t[:, ixs].flatten())
                 for name, ixs in self._subdomains.items()
@@ -237,29 +237,29 @@ def _find_nz(self, rows, cols, shape, transform=None):
     def _adaptive(self, marked):
         """Longest edge bisection."""
         if isinstance(marked, list):
-            marked = np.array(marked, dtype=np.int64)
+            marked = np.array(marked, dtype=np.int32)
         marked = np.unique(marked)
         nt = self.t.shape[1]
         nv = self.p.shape[1]
         p = np.zeros((3, 9 * nv), dtype=np.float64)
-        t = np.zeros((4, 8 * nt), dtype=np.int64)
+        t = np.zeros((4, 8 * nt), dtype=np.int32)
         p[:, :nv] = self.p.copy()
         t[:, :nt] = self.t.copy()
 
         nonconf = np.ones(8 * nv, dtype=np.int8)
-        split_edge = np.zeros((3, 8 * nv), dtype=np.int64)
+        split_edge = np.zeros((3, 8 * nv), dtype=np.int32)
         ns = 0
 
         while len(marked) > 0:
             nm = len(marked)
-            tnew = np.zeros(nm, dtype=np.int64) - 1
+            tnew = np.zeros(nm, dtype=np.int32) - 1
             t = self._adaptive_sort_mesh(p, t, marked)
             t0, t1, t2, t3 = t[:, marked]
 
             if ns == 0:
-                ix = np.arange(nm, dtype=np.int64)
+                ix = np.arange(nm, dtype=np.int32)
             else:
-                nonconf_edge = np.nonzero(nonconf[:ns])[0]
+                nonconf_edge = np.nonzero(nonconf[:ns])[0].astype(np.int32)
                 i, j = self._find_nz(
                     np.hstack((split_edge[0, nonconf_edge],
                                split_edge[1, nonconf_edge])),
@@ -268,7 +268,7 @@ def _adaptive(self, marked):
                     lambda I: I[t0].multiply(I[t1])
                 )
                 tnew[i] = j
-                ix = np.nonzero(tnew == -1)[0]
+                ix = np.nonzero(tnew == -1)[0].astype(np.int32)
 
             if len(ix) > 0:
                 i, j = self._find_nz(
@@ -280,7 +280,7 @@ def _adaptive(self, marked):
 
                 split_edge[0, nix] = i
                 split_edge[1, nix] = j
-                split_edge[2, nix] = np.arange(nv, nv + nn, dtype=np.int64)
+                split_edge[2, nix] = np.arange(nv, nv + nn, dtype=np.int32)
 
                 # add new points
                 p[:, nv:(nv + nn)] = .5 * (p[:, i] + p[:, j])
@@ -301,10 +301,11 @@ def _adaptive(self, marked):
             t[:, nt:(nt + nm)] = np.vstack((t2, t1, t3, tnew))
             nt += nm
 
-            check = np.nonzero(nonconf[:ns])[0]
-            check_node = np.zeros(nv, dtype=np.int64)
+            check = np.nonzero(nonconf[:ns])[0].astype(np.int32)
+            check_node = np.zeros(nv, dtype=np.int32)
             check_node[split_edge[:2, check]] = 1
-            check_elem = np.nonzero(check_node[t[:, :nt]].sum(axis=0))[0]
+            check_elem = (np.nonzero(check_node[t[:, :nt]].sum(axis=0))[0]
+                          .astype(np.int32))
 
             i, j = self._find_nz(
                 t[:, check_elem],
@@ -389,7 +390,7 @@ def init_tensor(cls: Type,
         T[:, (4 * ne):(5 * ne)] = t[[0, 2, 6, 7]]
         T[:, (5 * ne):] = t[[0, 3, 6, 7]]
 
-        return cls(p, T.astype(np.int64))
+        return cls(p, T.astype(np.int32))
 
     @classmethod
     def init_ball(cls: Type,
@@ -416,7 +417,7 @@ def init_ball(cls: Type,
                       [0, 2, 3, 4],
                       [0, 4, 5, 3],
                       [0, 4, 6, 2],
-                      [0, 5, 6, 1]], dtype=np.int64).T
+                      [0, 5, 6, 1]], dtype=np.int32).T
         m = cls(p, t)
         for _ in range(nrefs):
             m = m.refined()
diff --git a/skfem/mesh/mesh_tri_1.py b/skfem/mesh/mesh_tri_1.py
index 8c227dd8..f5458f75 100644
--- a/skfem/mesh/mesh_tri_1.py
+++ b/skfem/mesh/mesh_tri_1.py
@@ -20,7 +20,7 @@ class MeshTri1(MeshSimplex, Mesh2D):
     )
     t: ndarray = field(
         default_factory=lambda: np.array(
-            [[0, 1, 2], [1, 3, 2]], dtype=np.int64
+            [[0, 1, 2], [1, 3, 2]], dtype=np.int32
         ).T
     )
     elem: Type[Element] = ElementTriP1
@@ -76,7 +76,7 @@ def init_tensor(cls: Type, x: ndarray, y: ndarray):
                      .copy()
                      .flatten())
 
-        return cls(p, t.astype(np.int64))
+        return cls(p, t.astype(np.int32))
 
     @classmethod
     def init_symmetric(cls: Type) -> Mesh2D:
@@ -100,7 +100,7 @@ def init_symmetric(cls: Type) -> Mesh2D:
         t = np.array([[0, 1, 4],
                       [1, 2, 4],
                       [2, 3, 4],
-                      [0, 3, 4]], dtype=np.int64).T
+                      [0, 3, 4]], dtype=np.int32).T
         return cls(p, t)
 
     @classmethod
@@ -129,7 +129,7 @@ def init_sqsymmetric(cls: Type) -> Mesh2D:
                       [3, 4, 6],
                       [4, 6, 7],
                       [4, 7, 8],
-                      [4, 5, 8]], dtype=np.int64).T
+                      [4, 5, 8]], dtype=np.int32).T
         return cls(p, t)
 
     @classmethod
@@ -156,7 +156,7 @@ def init_lshaped(cls: Type) -> Mesh2D:
                       [0, 6, 3],
                       [0, 7, 4],
                       [0, 4, 5],
-                      [0, 3, 5]], dtype=np.int64).T
+                      [0, 3, 5]], dtype=np.int32).T
         return cls(p, t)
 
     @classmethod
@@ -194,7 +194,7 @@ def init_circle(cls: Type,
         t = np.array([[0, 1, 2],
                       [0, 1, 4],
                       [0, 2, 3],
-                      [0, 3, 4]], dtype=np.int64).T
+                      [0, 3, 4]], dtype=np.int32).T
         m = cls(p, t)
         for _ in range(nrefs):
             m = m.refined()
@@ -228,8 +228,8 @@ def _uniform(self):
 
         if self._boundaries is not None:
             # mapping of indices between old and new facets
-            new_facets = np.zeros((2, self.facets.shape[1]), dtype=np.int64)
-            ix0 = np.arange(t.shape[1], dtype=np.int64)
+            new_facets = np.zeros((2, self.facets.shape[1]), dtype=np.int32)
+            ix0 = np.arange(t.shape[1], dtype=np.int32)
             ix1 = ix0 + t.shape[1]
             ix2 = ix0 + 2 * t.shape[1]
 
@@ -277,7 +277,7 @@ def _adaptive_sort_mesh(p, t):
     @staticmethod
     def _adaptive_find_facets(m, marked_elems):
         """Find the facets to split."""
-        facets = np.zeros(m.facets.shape[1], dtype=np.int64)
+        facets = np.zeros(m.facets.shape[1], dtype=np.int32)
         facets[m.t2f[:, marked_elems].flatten('F')] = 1
         prev_nnz = -1e10
 
@@ -292,7 +292,7 @@ def _adaptive_find_facets(m, marked_elems):
     @staticmethod
     def _adaptive_split_elements(m, facets, subdomains):
         """Define new elements."""
-        ix = (-1) * np.ones(m.facets.shape[1], dtype=np.int64)
+        ix = (-1) * np.ones(m.facets.shape[1], dtype=np.int32)
         ix[facets == 1] = (np.arange(np.count_nonzero(facets))
                            + m.p.shape[1])
         ix = ix[m.t2f]
@@ -335,28 +335,28 @@ def _adaptive_split_elements(m, facets, subdomains):
                   m.p[:, m.facets[1, facets == 1]])
 
         if subdomains is not None:
-            new_t = np.zeros((4, m.t.shape[1]), dtype=np.int64) - 1
+            new_t = np.zeros((4, m.t.shape[1]), dtype=np.int32) - 1
             nred = np.sum(red)
             nblue1 = np.sum(blue1)
             nblue2 = np.sum(blue2)
             ngreen = np.sum(green)
             offset = np.sum(rest)
-            new_t[0, rest] = np.arange(offset, dtype=np.int64)
+            new_t[0, rest] = np.arange(offset, dtype=np.int32)
             new_t[:, red] = np.arange(offset,
                                       offset + 4 * nred,
-                                      dtype=np.int64).reshape(4, -1)
+                                      dtype=np.int32).reshape(4, -1)
             offset += 4 * nred
             new_t[:3, blue1] = np.arange(offset,
                                          offset + 3 * nblue1,
-                                         dtype=np.int64).reshape(3, -1)
+                                         dtype=np.int32).reshape(3, -1)
             offset += 3 * nblue1
             new_t[:3, blue2] = np.arange(offset,
                                          offset + 3 * nblue2,
-                                         dtype=np.int64).reshape(3, -1)
+                                         dtype=np.int32).reshape(3, -1)
             offset += 3 * nblue2
             new_t[:2, green] = np.arange(offset,
                                          offset + 2 * ngreen,
-                                         dtype=np.int64).reshape(2, -1)
+                                         dtype=np.int32).reshape(2, -1)
             subdomains = {
                 name: np.setdiff1d(np.unique(new_t[:, ixs]), [-1])
                 for name, ixs in subdomains.items()
@@ -397,7 +397,7 @@ def __mul__(self, other):
 
         if isinstance(other, MeshLine1):
             points = np.zeros((3, 0), dtype=np.float64)
-            wedges = np.zeros((6, 0), dtype=np.int64)
+            wedges = np.zeros((6, 0), dtype=np.int32)
             diff = 0
             for i, p in enumerate(np.sort(other.p[0])):
                 points = np.hstack((
@@ -438,7 +438,7 @@ def finder(x, y, _search_all=False):
                 _, ix_ind = np.unique(ix, return_index=True)
                 ix = ix[np.sort(ix_ind)]
             else:
-                ix = np.arange(nelems, dtype=np.int64)
+                ix = np.arange(nelems, dtype=np.int32)
 
             X = mapping.invF(np.array([x, y])[:, None], ix)
             eps = np.finfo(X.dtype).eps
diff --git a/skfem/mesh/mesh_tri_1_dg.py b/skfem/mesh/mesh_tri_1_dg.py
index 06c2c53a..78aca088 100644
--- a/skfem/mesh/mesh_tri_1_dg.py
+++ b/skfem/mesh/mesh_tri_1_dg.py
@@ -32,7 +32,7 @@ class MeshTri1DG(MeshDG, MeshTri1):
                 [0, 1, 3],
                 [0, 2, 3],
             ],
-            dtype=np.int64,
+            dtype=np.int32,
         ).T
     )
     elem: Type[Element] = ElementTriP1DG
diff --git a/skfem/mesh/mesh_wedge_1.py b/skfem/mesh/mesh_wedge_1.py
index 64711d6b..f275246b 100644
--- a/skfem/mesh/mesh_wedge_1.py
+++ b/skfem/mesh/mesh_wedge_1.py
@@ -29,7 +29,7 @@ class MeshWedge1(Mesh3D):
     )
     t: ndarray = field(
         default_factory=lambda: np.array(
-            [[0, 2, 3, 1, 4, 5], [2, 3, 6, 4, 5, 7]], dtype=np.int64
+            [[0, 2, 3, 1, 4, 5], [2, 3, 6, 4, 5, 7]], dtype=np.int32
         ).T
     )
     elem: Type[Element] = ElementWedge1
diff --git a/skfem/refdom.py b/skfem/refdom.py
index 33025d30..ace1fc02 100644
--- a/skfem/refdom.py
+++ b/skfem/refdom.py
@@ -34,14 +34,14 @@ def on_facet(cls, i, X):
 class RefPoint(Refdom):
 
     p = np.array([[0.]], dtype=np.float64)
-    t = np.array([[0]], dtype=np.int64)
+    t = np.array([[0]], dtype=np.int32)
     name = "Zero-dimensional"
 
 
 class RefLine(Refdom):
 
     p = np.array([[0., 1.]], dtype=np.float64)
-    t = np.array([[0], [1]], dtype=np.int64)
+    t = np.array([[0], [1]], dtype=np.int32)
     normals = np.array([[-1.],
                         [1.]])
     facets = [[0],
@@ -56,7 +56,7 @@ class RefTri(Refdom):
 
     p = np.array([[0., 1., 0.],
                   [0., 0., 1.]], dtype=np.float64)
-    t = np.array([[0], [1], [2]], dtype=np.int64)
+    t = np.array([[0], [1], [2]], dtype=np.int32)
     normals = np.array([[0., -1.],
                         [1., 1.],
                         [-1., 0.]])
@@ -93,7 +93,7 @@ class RefTet(Refdom):
     p = np.array([[0., 1., 0., 0.],
                   [0., 0., 1., 0.],
                   [0., 0., 0., 1.]], dtype=np.float64)
-    t = np.array([[0], [1], [2], [3]], dtype=np.int64)
+    t = np.array([[0], [1], [2], [3]], dtype=np.int32)
     facets = [[0, 1, 2],
               [0, 1, 3],
               [0, 2, 3],
@@ -153,7 +153,7 @@ class RefQuad(Refdom):
 
     p = np.array([[0., 1., 1., 0.],
                   [0., 0., 1., 1.]], dtype=np.float64)
-    t = np.array([[0], [1], [2], [3]], dtype=np.int64)
+    t = np.array([[0], [1], [2], [3]], dtype=np.int32)
     normals = np.array([[0., -1.],
                         [1., 0.],
                         [0., 1.],
@@ -178,7 +178,7 @@ class RefHex(Refdom):
                   [0., 1., 0.],
                   [0., 0., 1.],
                   [0., 0., 0.]], dtype=np.float64).T
-    t = np.array([[0], [1], [2], [3], [4], [5], [6], [7]], dtype=np.int64)
+    t = np.array([[0], [1], [2], [3], [4], [5], [6], [7]], dtype=np.int32)
     normals = np.array([[1., 0., 0.],
                         [0., 0., 1.],
                         [0., 1., 0.],
@@ -264,7 +264,7 @@ class RefWedge(Refdom):
                   [0., 0., 1.],
                   [1., 0., 1.],
                   [0., 1., 1.]], dtype=np.float64).T
-    t = np.array([[0], [1], [2], [3], [4], [5]], dtype=np.int64)
+    t = np.array([[0], [1], [2], [3], [4], [5]], dtype=np.int32)
     normals = np.array([[0., -1., 0.],
                         [1., 1., 0.],
                         [-1., 0., 0.],
diff --git a/skfem/supermeshing.py b/skfem/supermeshing.py
index 322aa619..02324ce6 100644
--- a/skfem/supermeshing.py
+++ b/skfem/supermeshing.py
@@ -62,8 +62,8 @@ def _intersect2d(p1, t1, p2, t2):
                 ix2.append(jtr)
     return (
         MeshTri(p, t),
-        np.array(ix1, dtype=np.int64),
-        np.array(ix2, dtype=np.int64),
+        np.array(ix1, dtype=np.int32),
+        np.array(ix2, dtype=np.int32),
     )
 
 
diff --git a/skfem/utils.py b/skfem/utils.py
index e4c18ef2..92358d4d 100644
--- a/skfem/utils.py
+++ b/skfem/utils.py
@@ -307,9 +307,9 @@ def _init_bc(A: spmatrix,
     if I is None and D is None:
         raise Exception("Either I or D must be given!")
     elif I is None and D is not None:
-        I = np.setdiff1d(np.arange(A.shape[0], dtype=np.int64), D)
+        I = np.setdiff1d(np.arange(A.shape[0], dtype=np.int32), D)
     elif D is None and I is not None:
-        D = np.setdiff1d(np.arange(A.shape[0], dtype=np.int64), I)
+        D = np.setdiff1d(np.arange(A.shape[0], dtype=np.int32), I)
     else:
         raise Exception("Give only I or only D!")
 
@@ -377,7 +377,7 @@ def enforce(A: spmatrix,
     start = Aout.indptr[D]
     stop = Aout.indptr[D + 1]
     count = stop - start
-    idx = np.ones(count.sum(), dtype=np.int64)
+    idx = np.ones(count.sum(), dtype=np.int32)
     idx[np.cumsum(count)[:-1]] -= count[:-1]
     idx = np.repeat(start, count) + np.cumsum(idx) - 1
     Aout.data[idx] = 0.
@@ -624,11 +624,11 @@ def mpc(A: spmatrix,
 
     """
     if M is None:
-        M = np.array([], dtype=np.int64)
+        M = np.array([], dtype=np.int32)
     if S is None:
-        S = np.array([], dtype=np.int64)
+        S = np.array([], dtype=np.int32)
 
-    U = np.setdiff1d(np.arange(A.shape[0], dtype=np.int64),
+    U = np.setdiff1d(np.arange(A.shape[0], dtype=np.int32),
                      np.concatenate((M, S)))
 
     if T is None:
@@ -705,9 +705,9 @@ def rcm(A: spmatrix,
 def adaptive_theta(est, theta=0.5, max=None):
     """For choosing which elements to refine in an adaptive strategy."""
     if max is None:
-        return np.nonzero(theta * np.max(est) < est)[0]
+        return np.nonzero(theta * np.max(est) < est)[0].astype(np.int32)
     else:
-        return np.nonzero(theta * max < est)[0]
+        return np.nonzero(theta * max < est)[0].astype(np.int32)
 
 
 @deprecated("Basis.project")
diff --git a/skfem/visuals/glvis.py b/skfem/visuals/glvis.py
index fd3ecab0..c6b23b76 100644
--- a/skfem/visuals/glvis.py
+++ b/skfem/visuals/glvis.py
@@ -80,15 +80,15 @@ def plot(basis, x, keys: Optional[str] = None):
         m.dim(),
         m.nelements,
         _to_int_string(np.hstack((
-            np.ones(m.nelements, dtype=np.int64)[:, None],
-            (np.zeros(m.nelements, dtype=np.int64)[:, None]
+            np.ones(m.nelements, dtype=np.int32)[:, None],
+            (np.zeros(m.nelements, dtype=np.int32)[:, None]
              + MESH_TYPE_MAPPING[type(m)]),
             m.t.T,
         ))),
         nbfacets,
         _to_int_string(np.hstack((
-            np.ones(nbfacets, dtype=np.int64)[:, None],
-            (np.zeros(nbfacets, dtype=np.int64)[:, None]
+            np.ones(nbfacets, dtype=np.int32)[:, None],
+            (np.zeros(nbfacets, dtype=np.int32)[:, None]
              + BOUNDARY_TYPE_MAPPING[MESH_TYPE_MAPPING[type(m)]]),
             m.facets[:, bfacets].T,
         ))),
diff --git a/tests/test_dofs.py b/tests/test_dofs.py
index b7f72a02..3cc8abab 100644
--- a/tests/test_dofs.py
+++ b/tests/test_dofs.py
@@ -28,7 +28,7 @@ def runTest(self):
         assert_allclose(all_dofs,
                         all_dofs.drop('does_not_exist'))
 
-        assert_allclose(np.empty((0,), dtype=np.int64),
+        assert_allclose(np.empty((0,), dtype=np.int32),
                         all_dofs.keep('does_not_exist'))
 
 
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 87a355ef..8d7c1ef3 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -242,7 +242,7 @@ def test_adaptive_splitting_3d_1():
 def test_adaptive_splitting_3d_2():
     m = MeshTet()
     for itr in range(5):
-        m = m.refined(np.arange(m.nelements, dtype=np.int64))
+        m = m.refined(np.arange(m.nelements, dtype=np.int32))
         assert m.is_valid()
 
 def test_adaptive_splitting_3d_3():