use int32 everywhere

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
 

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``).

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

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:

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)]
 

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)
 

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)
 

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):

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):

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))

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):

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]

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]

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):