Skip to content

Commit

Permalink
use int32 everywhere (#1135)
Browse files Browse the repository at this point in the history
  • Loading branch information
kinnala authored Jun 19, 2024
1 parent 61e802e commit 54ffefb
Showing 35 changed files with 196 additions and 175 deletions.
5 changes: 5 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
@@ -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

8 changes: 4 additions & 4 deletions docs/advanced.rst
Original file line number Diff line number Diff line change
@@ -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``).
2 changes: 1 addition & 1 deletion docs/examples/ex41.py
Original file line number Diff line number Diff line change
@@ -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
12 changes: 6 additions & 6 deletions docs/howto.rst
Original file line number Diff line number Diff line change
@@ -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:
18 changes: 9 additions & 9 deletions skfem/assembly/basis/abstract_basis.py
Original file line number Diff line number Diff line change
@@ -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)]

2 changes: 1 addition & 1 deletion skfem/assembly/basis/facet_basis.py
Original file line number Diff line number Diff line change
@@ -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)

2 changes: 1 addition & 1 deletion skfem/assembly/basis/interior_facet_basis.py
Original file line number Diff line number Diff line change
@@ -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)

36 changes: 18 additions & 18 deletions skfem/assembly/dofs.py
Original file line number Diff line number Diff line change
@@ -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):
4 changes: 2 additions & 2 deletions skfem/assembly/form/bilinear_form.py
Original file line number Diff line number Diff line change
@@ -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):
2 changes: 1 addition & 1 deletion skfem/assembly/form/linear_form.py
Original file line number Diff line number Diff line change
@@ -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))
6 changes: 3 additions & 3 deletions skfem/assembly/form/trilinear_form.py
Original file line number Diff line number Diff line change
@@ -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):
2 changes: 1 addition & 1 deletion skfem/element/element_hcurl.py
Original file line number Diff line number Diff line change
@@ -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]
2 changes: 1 addition & 1 deletion skfem/element/element_hdiv.py
Original file line number Diff line number Diff line change
@@ -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]
6 changes: 3 additions & 3 deletions skfem/experimental/autodiff/__init__.py
Original file line number Diff line number Diff line change
@@ -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):
Loading

0 comments on commit 54ffefb

Please sign in to comment.