backport of #837, #838, #834, distutils fix (#835)

.github/workflows/test.yaml

@@ -100,7 +100,7 @@ jobs:
       - name: Configure MKL
         if: ${{ matrix.matrix-backend == 'mkl' }}
         run: |
-          python -um pip install --upgrade --upgrade-strategy eager mkl
+          python -um pip install --upgrade --upgrade-strategy eager 'mkl<2024'
           python -um devtools.gha.configure_mkl
       - name: Test
         run: python -um coverage run -m unittest discover -b -q -t . -s tests

docs/sphinx_mods.py

@@ -78,7 +78,7 @@ def generate_api(app):
     dst_root.mkdir(parents=True, exist_ok=True)
 
     srcs = tuple(f for f in sorted(nutils.glob('**/*.py')) if f != nutils/'__init__.py' and (f.name == '__init__.py' or not f.name.startswith('_')))
-    for src in sphinx.util.status_iterator(srcs, 'generating api... ', 'purple', len(srcs), app.verbosity):
+    for src in sphinx.util.display.status_iterator(srcs, 'generating api... ', 'purple', len(srcs), app.verbosity):
         module = '.'.join((src.parent if src.name == '__init__.py' else src.with_suffix('')).relative_to(nutils).parts)
         dst = dst_root/(module+'.rst')
         with dst.open('w', encoding='utf-8') as f:

nutils/__init__.py

@@ -1,13 +1,6 @@
 'Numerical Utilities for Finite Element Analysis'
 
-import sys
-import numpy
-from distutils.version import LooseVersion
-
-assert sys.version_info >= (3, 5)
-assert LooseVersion(numpy.version.version) >= LooseVersion('1.16'), 'nutils requires numpy 1.16 or higher, got {}'.format(numpy.version.version)
-
-__version__ = version = '8.3'
+__version__ = version = '8.4'
 version_name = 'idiyappam'
 long_version = ('{} "{}"' if version_name else '{}').format(version, version_name)
 

nutils/function.py

@@ -2733,12 +2733,6 @@ def get(__array: IntoArray, __axis: int, __index: IntoArray) -> Array:
     return numpy.take(Array.cast(__array), Array.cast(__index, dtype=int, ndim=0), __axis)
 
 
-def _range(__length: int, __offset: int) -> Array:
-    length = Array.cast(__length, dtype=int, ndim=0)
-    offset = Array.cast(__offset, dtype=int, ndim=0)
-    return _Wrapper(lambda l, o: evaluable.Range(l) + o, _WithoutPoints(length), _WithoutPoints(offset), shape=(__length,), dtype=int)
-
-
 def _takeslice(__array: IntoArray, __s: slice, __axis: int) -> Array:
     array = Array.cast(__array)
     s = __s
@@ -2749,7 +2743,8 @@ def _takeslice(__array: IntoArray, __s: slice, __axis: int) -> Array:
         stop = n if s.stop is None else s.stop if s.stop >= 0 else s.stop + n
         if start == 0 and stop == n:
             return array
-        index = _range(stop-start, start)
+        length = stop - start
+        index = _Wrapper(evaluable.Range, _WithoutPoints(_Constant(length)), shape=(length,), dtype=int) + start
     elif isinstance(n, numbers.Integral):
         index = Array.cast(numpy.arange(*s.indices(int(n))))
     else:

nutils/matrix/_mkl.py

@@ -30,8 +30,8 @@ def assemble(data, index, shape):
 class Pardiso:
     '''Wrapper for libmkl.pardiso.
 
-    https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/
-      sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface.html
+    https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/
+      onemkl-pardiso-parallel-direct-sparse-solver-iface.html
     '''
 
     _errorcodes = {

nutils/topology.py

@@ -1866,6 +1866,19 @@ def __or__(self, other):
     def __rsub__(self, other):
         return other
 
+    @property
+    def connectivity(self):
+        return ()
+
+    @property
+    def boundary(self):
+        if self.ndims:
+            return EmptyTopology(self.space, self.transforms.todims, self.ndims - 1)
+        else:
+            raise ValueError('A 0D topology has no boundary.')
+
+    interfaces = boundary
+
 
 def StructuredLine(space, root: transform.TransformItem, i: int, j: int, periodic: bool = False, bnames: Optional[Tuple[str, str]] = None):
     assert isinstance(i, int), f'i={i!r}'
@@ -2350,22 +2363,28 @@ def __init__(self, space: str, references: References, transforms: transformseq.
 class SimplexTopology(TransformChainsTopology):
     'simpex topology'
 
-    def _renumber(simplices):
-        simplices = numpy.asarray(simplices)
-        keep = numpy.zeros(simplices.max()+1, dtype=bool)
-        keep[simplices.flat] = True
-        return types.arraydata(simplices if keep.all() else (numpy.cumsum(keep)-1)[simplices])
-
     def __init__(self, space: str, simplices: numpy.ndarray, transforms: transformseq.Transforms, opposites: transformseq.Transforms):
         assert isinstance(space, str), f'space={space!r}'
         assert isinstance(simplices, numpy.ndarray), f'simplices={simplices!r}'
+        assert len(simplices), f'simplices is empty'
         assert simplices.shape == (len(transforms), transforms.fromdims+1)
         self.simplices = numpy.asarray(simplices)
         assert numpy.greater(self.simplices[:, 1:], self.simplices[:, :-1]).all(), 'nodes should be sorted'
         assert not numpy.equal(self.simplices[:, 1:], self.simplices[:, :-1]).all(), 'duplicate nodes'
         references = References.uniform(element.getsimplex(transforms.fromdims), len(transforms))
         super().__init__(space, references, transforms, opposites)
 
+    @cached_property
+    def contiguous_simplices(self):
+        simplices = numpy.asarray(self.simplices)
+        keep = numpy.zeros(simplices.max()+1, dtype=bool)
+        keep[simplices.flat] = True
+        return simplices if keep.all() else (numpy.cumsum(keep)-1)[simplices]
+
+    @cached_property
+    def nverts(self):
+        return self.contiguous_simplices.max() + 1
+
     def take_unchecked(self, indices):
         space, = self.spaces
         return SimplexTopology(space, self.simplices[indices], self.transforms[indices], self.opposites[indices])
@@ -2374,6 +2393,8 @@ def take_unchecked(self, indices):
     def boundary(self):
         space, = self.spaces
         ielem, iedge = (self.connectivity == -1).nonzero()
+        if len(ielem) == 0:
+            return self.empty_like().boundary
         nd = self.ndims
         edges = numpy.arange(nd+1).repeat(nd).reshape(nd,nd+1).T[::-1]
         simplices = self.simplices[ielem, edges[iedge].T].T
@@ -2398,7 +2419,7 @@ def connectivity(self):
     def basis_std(self, degree):
         if degree == 1:
             coeffs = element.getsimplex(self.ndims).get_poly_coeffs('bernstein', degree=1)
-            return function.PlainBasis([coeffs] * len(self), self.simplices, self.simplices.max()+1, self.f_index, self.f_coords)
+            return function.PlainBasis([coeffs] * len(self), self.contiguous_simplices, self.nverts, self.f_index, self.f_coords)
         return super().basis_std(degree)
 
     def basis_bubble(self):
@@ -2410,9 +2431,8 @@ def basis_bubble(self):
         coeffs[:-1] = poly.change_degree(bernstein, self.ndims, 1 + self.ndims) - bubble[None] / (self.ndims+1)
         coeffs[-1] = bubble
         coeffs = types.frozenarray(coeffs, copy=False)
-        nverts = self.simplices.max() + 1
-        ndofs = nverts + len(self)
-        nmap = [types.frozenarray(numpy.hstack([idofs, nverts+ielem]), copy=False) for ielem, idofs in enumerate(self.simplices)]
+        ndofs = self.nverts + len(self)
+        nmap = [types.frozenarray(numpy.hstack([idofs, self.nverts+ielem]), copy=False) for ielem, idofs in enumerate(self.contiguous_simplices)]
         return function.PlainBasis([coeffs] * len(self), nmap, ndofs, self.f_index, self.f_coords)
 
 

pyproject.toml

@@ -24,7 +24,7 @@ classifiers = [
 [project.optional-dependencies]
 docs = ["Sphinx>=1.8"]
 export_mpl = ["matplotlib>=1.3", "pillow>2.6"]
-matrix_mkl = ["mkl"]
+matrix_mkl = ["mkl<2024"]
 matrix_scipy = ["scipy>=0.13"]
 import_gmsh = ["meshio"]
 

tests/test_function.py

@@ -6,6 +6,7 @@
 import warnings as _builtin_warnings
 import functools
 import fractions
+import pickle
 
 
 class Array(TestCase):
@@ -171,6 +172,12 @@ def test_lower_eval(self):
         desired = self.n_op(*self.args)
         self.assertArrayAlmostEqual(actual, desired, decimal=15)
 
+    def test_pickle(self):
+        f = self.op(*self.args)
+        s = pickle.dumps(f)
+        f_ = pickle.loads(s)
+        self.assertEqual(f.as_evaluable_array, f_.as_evaluable_array)
+
 
 def generate(*shape, real, imag, zero, negative):
     'generate array values that cover certain numerical classes'

tests/test_topology.py

@@ -1352,9 +1352,30 @@ def setUp(self):
         self.desired_references = [element.TriangleReference()]*4
         self.desired_vertices = coords[simplices].tolist()
 
+    def test_contiguous_simplices(self):
+        # Same structure as the simplices in `setUp`, but renumbered.
+        simplices = numpy.array([[10,13,25],[10,18,25],[13,21,25],[18,21,25]])
+        transforms = transformseq.IndexTransforms(2, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(self.topo.contiguous_simplices.tolist(), self.topo.simplices.tolist())
+        self.assertEqual(self.topo.nverts, 5)
+
     def test_boundary(self):
         self.assertIsInstance(self.topo.boundary, topology.SimplexTopology)
 
+    def test_empty_boundary(self):
+        simplices = numpy.array([[0, 1, 2, 3]])
+        transforms = transformseq.IndexTransforms(3, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(len(topo.boundary), 4)
+        self.assertEqual(len(topo.boundary.boundary), 0)
+
+    def test_empty_interfaces(self):
+        simplices = numpy.array([[0, 1, 2, 3]])
+        transforms = transformseq.IndexTransforms(3, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(len(topo.interfaces), 0)
+
     def test_getitem(self):
         self.assertIsInstance(self.topo[numpy.arange(4) < 2], topology.SimplexTopology)
         self.assertIsInstance(self.topo[numpy.arange(2)], topology.SimplexTopology)