diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 4a0f84e1d..60851b45b 100644 --- a/.github/workflows/test.yaml +++ b/.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 diff --git a/docs/sphinx_mods.py b/docs/sphinx_mods.py index 3996c325d..e1e17e2e2 100644 --- a/docs/sphinx_mods.py +++ b/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: diff --git a/nutils/__init__.py b/nutils/__init__.py index aabc1fb0a..a0d899be6 100644 --- a/nutils/__init__.py +++ b/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) diff --git a/nutils/function.py b/nutils/function.py index a853167ff..0e279369c 100644 --- a/nutils/function.py +++ b/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: diff --git a/nutils/matrix/_mkl.py b/nutils/matrix/_mkl.py index 45a2f2acf..5c116b9e9 100644 --- a/nutils/matrix/_mkl.py +++ b/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 = { diff --git a/nutils/topology.py b/nutils/topology.py index 31b363cec..8212a0fc2 100644 --- a/nutils/topology.py +++ b/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,15 +2363,10 @@ 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' @@ -2366,6 +2374,17 @@ def __init__(self, space: str, simplices: numpy.ndarray, transforms: transformse 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) diff --git a/pyproject.toml b/pyproject.toml index ac1d91257..0ce5a3667 100644 --- a/pyproject.toml +++ b/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"] diff --git a/tests/test_function.py b/tests/test_function.py index 7f7480b0d..5c26206bf 100644 --- a/tests/test_function.py +++ b/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' diff --git a/tests/test_topology.py b/tests/test_topology.py index 5d1308edf..ab4e0920a 100644 --- a/tests/test_topology.py +++ b/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)