From 4d4d544fc819e5b576baf275a75b0dd351882c98 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 28 Jan 2025 10:04:13 +0100 Subject: [PATCH 01/10] Fix merging of nested Assemble operations This patch fixes Assemble._optimized_for_numpy, which was incorrect in certain scenarios involving higher dimensional indices. A unit tests triggering this error follows in the text commit as tests.test_evaluable.check:inflate-twice. --- nutils/evaluable.py | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index b33de5f3c..06faa7a7b 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -3422,22 +3422,27 @@ def _compile_with_out(self, builder, out, out_block_id, mode): def _optimized_for_numpy(self): if isinstance(self.func, Assemble): - nontrivial_indices = [] - for f in self.func, self: # order is important to maintain insertion order of ndim=0 indices - offset = 0 - for index in f.indices: - if index != Range(f.shape[offset]): - nontrivial_indices.append((offset, index)) - offset += index.ndim - assert offset == f.func.ndim - nontrivial_indices.sort(key=lambda item: item[0]) # order of equal offsets is maintained - indices = [] - for i, index in nontrivial_indices: - if i < len(indices): - return # inflations overlap - while i > len(indices): - indices.append(Range(self.shape[len(indices)])) - indices.append(index) + indices = list(self.indices) + # We aim to merge the indices from the nested Assemble operations + # if they are separable, i.e. preceded by or following on full + # slices, by replacing Range instances in indices by the + # corresponding index from self.func. + for i, index in enumerate(self.func.indices): + if index != Range(self.func.shape[i]): # non-trivial index of self.func + # we need to account for the different axis numberings + # between self and self.func to find the right insertion + # point. + ax1 = 0 # axis of self.func + ax2 = 0 # axis of self + while ax1 < i: # find ax1, ax2 corresponding to i + ax1 += indices[ax2].ndim + ax2 += 1 + if ax1 != i or ax2 >= self.ndim or indices[ax2] != Range(self.shape[ax2]): + # Any nontrivial nesting scenario would have been + # handled by Inflate if possible, so we simply bail out + # at the first sign of difficulty. + return + indices[ax2] = index # merge! return Assemble(self.func.func, tuple(indices), self.shape) From 8efe8f93c532733150be6576c3c8684f60923c43 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 27 Jan 2025 15:41:28 +0100 Subject: [PATCH 02/10] Fix Assemble._compile_with_out This patch fixes the compiler of function.Assemble, which failed to account for all of the function's axes in case any of the indices was multi-dimensional. A unit tests triggering this error is added as tests.test_evaluable.check:inflate-twice. --- nutils/evaluable.py | 31 ++++++++++++++++++++++--------- tests/test_evaluable.py | 1 + 2 files changed, 23 insertions(+), 9 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 06faa7a7b..3ec07028c 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -3394,27 +3394,40 @@ def evalf(func, *args): return numeric.accumulate(func, reshaped_indices, shape) def _compile_with_out(self, builder, out, out_block_id, mode): + # Compiles to an assignment (or in place addition) of the form: + # + # out[index1, index2, ...] = func.transpose(trans) + # + # The order of the indices is unchanged, but ranges are converted to + # slices and any other indices (dubbed 'advanced') reshaped in similar + # fashion to numpy._ix to index the cross product. The right hand side + # is transposed if a slice operation separates two advanced indices to + # match Numpy's advanced indexing rules. assert mode in ('iadd', 'assign') if mode == 'assign': builder.get_block_for_evaluable(self, block_id=out_block_id, comment='zero').array_fill_zeros(out) - compiled_indices = [] advanced_ndim = builtins.sum(index.ndim for index in self.indices if not isinstance(index, Range)) + compiled_indices = [] + trans = [] # axes of func corresponding to advanced indices i = 0 for index in self.indices: + j = i + index.ndim if isinstance(index, Range): n = builder.compile(index.shape[0]) - compiled_indices.append(_pyast.Variable('slice').call(n)) + compiled_index = _pyast.Variable('slice').call(n) else: - j = i + index.ndim + prefix = len(trans) + trans.extend(range(i, j)) + suffix = advanced_ndim - len(trans) compiled_index = builder.compile(index) if index.ndim < advanced_ndim: - compiled_index = compiled_index.get_item(_pyast.Raw(','.join(['None'] * i + [':'] * index.ndim + ['None'] * (advanced_ndim-j)))) - compiled_indices.append(compiled_index) - i = j - assert i == advanced_ndim + compiled_index = compiled_index.get_item(_pyast.Raw(','.join(['None'] * prefix + [':'] * index.ndim + ['None'] * suffix))) + compiled_indices.append(compiled_index) + i = j + assert i == self.func.ndim + assert len(trans) == advanced_ndim compiled_func = builder.compile(self.func) - trans = [i for i, index in enumerate(self.indices) if not isinstance(index, Range)] - if advanced_ndim and any(numpy.diff(trans) > 1): + if advanced_ndim > 1 and trans[-1] - trans[0] != advanced_ndim - 1: # trans is noncontiguous # see https://numpy.org/doc/stable/user/basics.indexing.html#combining-advanced-and-basic-indexing trans.extend(i for i, index in enumerate(self.indices) if isinstance(index, Range)) compiled_func = compiled_func.get_attr('transpose').call(*[_pyast.LiteralInt(i) for i in trans]) diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py index 338bb7784..7486e78a9 100644 --- a/tests/test_evaluable.py +++ b/tests/test_evaluable.py @@ -559,6 +559,7 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2): _check('inflate-diagonal', lambda f: evaluable.Inflate(evaluable.Inflate(f, evaluable.constant(1), evaluable.constant(3)), evaluable.constant(1), evaluable.constant(3)), lambda a: numpy.diag(numpy.array([0, a, 0])), numpy.array(.5)) _check('inflate-one', lambda f: evaluable.Inflate(f, evaluable.constant(0), evaluable.constant(1)), lambda a: numpy.array([a]), numpy.array(.5)) _check('inflate-range', lambda f: evaluable.Inflate(f, evaluable.Range(evaluable.constant(3)), evaluable.constant(3)), lambda a: a, ANY(3)) +_check('inflate-twice', lambda f: evaluable._inflate(evaluable.Inflate(f, evaluable.constant([[0,1],[2,0]]), evaluable.constant(3)), evaluable.constant([0, 1, 0]), evaluable.constant(2), axis=0), lambda a: numpy.einsum('ij,jklm,lmn->ikn', [[1,0,1],[0,1,0]], a, [[[1,0,0],[0,1,0]],[[0,0,1],[1,0,0]]]), ANY(3,5,2,2)) _check('take', lambda f: evaluable.Take(f, evaluable.constant([0, 3, 2])), lambda a: a[:, [0, 3, 2]], ANY(2, 4)) _check('take-duplicate', lambda f: evaluable.Take(f, evaluable.constant([0, 3, 0])), lambda a: a[:, [0, 3, 0]], ANY(2, 4)) _check('choose', lambda a, b, c: evaluable.Choose(a % 2, evaluable.stack([b, c], -1)), lambda a, b, c: numpy.choose(a % 2, [b, c]), INT(3, 3), ANY(3, 3), ANY(3, 3)) From b8aa58c7256b01c68c4f706bfcf155619c5714d2 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 23 Jan 2025 15:59:09 +0100 Subject: [PATCH 03/10] Enhance linearize and replace_arguments This patch simplifies the function.linearize and function.replace_arguments functions by offloading argument processing to the helper function _argument_to_array. In the new form, linearize gains the ability to specity arguments directly as Argument objects, while replace_arguments gainst the ability to specify argument replacements as 'old:new' strings. --- nutils/function.py | 80 +++++++++++++++++++++++++++++++--------------- 1 file changed, 54 insertions(+), 26 deletions(-) diff --git a/nutils/function.py b/nutils/function.py index 3f74942a0..df334b6ea 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -911,24 +911,12 @@ class _Replace(Array): def __init__(self, arg: Array, replacements: Dict[str, Array]) -> None: self._arg = arg self._replacements = {} + for old, new in _argument_to_array(replacements, arg): + if new.spaces: + raise ValueError(f'replacement functions cannot contain spaces, but replacement for Argument {old.name!r} contains space {", ".join(new.spaces)}.') + self._replacements[old.name] = new # Build arguments map with replacements. - unreplaced = {} - for name, (shape, dtype) in arg.arguments.items(): - if name in replacements: - replacement = replacements[name] - if isinstance(replacement, str): - replacement = Argument(replacement, shape, dtype) - else: - replacement = Array.cast(replacement) - if replacement.shape != shape: - raise ValueError(f'Argument {name!r} has shape {shape} but the replacement has shape {replacement.shape}.') - if replacement.dtype != dtype: - raise ValueError(f'Argument {name!r} has dtype {dtype.__name__} but the replacement has dtype {replacement.dtype.__name__}.') - if replacement.spaces: - raise ValueError(f'replacement functions cannot contain spaces, but replacement for Argument {name!r} contains space {", ".join(replacement.spaces)}.') - self._replacements[name] = replacement - else: - unreplaced[name] = shape, dtype + unreplaced = {name: shape_dtype for name, shape_dtype in arg.arguments.items() if name not in replacements} arguments = _join_arguments([unreplaced] + [replacement.arguments for replacement in self._replacements.values()]) super().__init__(arg.shape, arg.dtype, arg.spaces, arguments) @@ -1687,6 +1675,53 @@ def kronecker(__array: IntoArray, axis: int, length: int, pos: IntoArray) -> Arr return _Transpose.from_end(scatter(__array, length, pos), axis) +def _argument_to_array(d: Any, array: Array) -> Iterable[Tuple[Argument, Array]]: + '''Helper function for argument replacement. + + Given dictionary-like input, along with an array to look up argument names, + yield ``(arg, new)`` tuples where ``arg`` is an ``Argument`` and ``new`` an + ``Array`` of equal shape and dtype. A ``ValueError`` is raised if any key + cannot be made a valid argument for the provided array or any value cannot + be made an array of the same shape and dtype. + + Dictionary-like input is either an actual dictionary, a sequence of (key, + value) tuples or strings, or a string, where any string is replaced by the + relevant Argument object. The following inputs are all equivalent: + + d = {'u': 'v', 'p': 'q'} + d = 'u:v,p:q' + d = ('u:v', 'p:q') + d = [('u', 'v'), ('p', 'q')] + d = {'u': Argument('v', ...), 'p': Argument('q', ...)} + d = [(Argument('u', ...), Argument('v', ...)), 'p:q'] + ''' + + for item in d.split(',') if isinstance(d, str) else d.items() if isinstance(d, dict) else d: + arg, new = item.split(':', 1) if isinstance(item, str) else item + + if isinstance(arg, str): + if arg not in array.arguments: + continue + arg = Argument(arg, *array.arguments[arg]) + elif not isinstance(arg, Argument): + raise ValueError('Key must be string or argument') + elif arg.name not in arguments: + continue + elif array.arguments[arg.name] != (arg.shape, arg.dtype): + raise ValueError(f'Argument {arg.name!r} has wrong shape or dtype') + + if isinstance(new, str): + new = Argument(new, arg.shape, arg.dtype) + else: + new = Array.cast(new) + if new.shape != arg.shape: + raise ValueError(f'Argument {arg.name!r} has shape {arg.shape} but the replacement has shape {new.shape}.') + elif new.dtype != arg.dtype: + raise ValueError(f'Argument {arg.name!r} has dtype {arg.dtype.__name__} but the replacement has dtype {new.dtype.__name__}.') + + yield arg, new + + @nutils_dispatch def replace_arguments(__array: IntoArray, __arguments: Mapping[str, Union[IntoArray, str]]) -> Array: '''Replace arguments with :class:`Array` objects. @@ -1735,15 +1770,8 @@ def linearize(__array: IntoArray, __arguments: Union[str, Dict[str, str], Iterab ''' array = Array.cast(__array) - args = __arguments.split(',') if isinstance(__arguments, str) \ - else __arguments.items() if isinstance(__arguments, dict) \ - else __arguments - parts = [] - for kv in args: - k, v = kv.split(':', 1) if isinstance(kv, str) else kv - f = derivative(array, k) - parts.append(numpy.sum(f * Argument(v, *array.arguments[k]), tuple(range(array.ndim, f.ndim)))) - return util.sum(parts) + return util.sum(numpy.sum(derivative(array, arg) * lin, array.ndim + numpy.arange(arg.ndim)) + for arg, lin in _argument_to_array(__arguments, array)) def broadcast_arrays(*arrays: IntoArray) -> Tuple[Array, ...]: From 171f205ced2261e70e36bc3343313d05262b8170 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 3 Feb 2025 14:36:09 +0100 Subject: [PATCH 04/10] Rephrase bound replacement error message This patch changes the error message that is raised when replacing an argument with a bound array from "replacement functions cannot contain spaces, but replacement for Argument .. contains space ..'" to "replacement functions cannot be bound to a space, but replacement for Argument .. is bound to ..". The new wording is thought to better reflect the situation with spaces. --- nutils/function.py | 2 +- tests/test_function.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/nutils/function.py b/nutils/function.py index df334b6ea..7d6c202c5 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -913,7 +913,7 @@ def __init__(self, arg: Array, replacements: Dict[str, Array]) -> None: self._replacements = {} for old, new in _argument_to_array(replacements, arg): if new.spaces: - raise ValueError(f'replacement functions cannot contain spaces, but replacement for Argument {old.name!r} contains space {", ".join(new.spaces)}.') + raise ValueError(f'replacement functions cannot be bound to a space, but replacement for Argument {old.name!r} is bound to {", ".join(new.spaces)}.') self._replacements[old.name] = new # Build arguments map with replacements. unreplaced = {name: shape_dtype for name, shape_dtype in arg.arguments.items() if name not in replacements} diff --git a/tests/test_function.py b/tests/test_function.py index 4fe6bb486..47846b813 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -790,7 +790,7 @@ def test_different_dtype(self): def test_nonempty_spaces(self): topo, geom = mesh.unitsquare(1, 'square') - with self.assertRaisesRegex(ValueError, "replacement functions cannot contain spaces, but replacement for Argument 'foo' contains space X."): + with self.assertRaisesRegex(ValueError, "replacement functions cannot be bound to a space, but replacement for Argument 'foo' is bound to X."): function.replace_arguments(function.Argument('foo', (2,), dtype=float), dict(foo=geom)) From 64b8a16e851858b4b3c8d55e63b029956d65c23e Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 24 Jan 2025 09:53:37 +0100 Subject: [PATCH 05/10] Allow iterator argument in function.vectorize This patch generalizes function.vectorize to accept any iterable argument, which makes it possible to use generator comprehensions without double bracketing. --- nutils/function.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nutils/function.py b/nutils/function.py index 7d6c202c5..90647c5f8 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -2346,6 +2346,7 @@ def vectorize(args: Sequence[IntoArray]) -> Array: :class:`Array` ''' + args = tuple(args) return numpy.concatenate([kronecker(arg, axis=-1, length=len(args), pos=iarg) for iarg, arg in enumerate(args)]) From 30148b3cc9597e8a1d1356e749e20f56e8d87dbf Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 24 Jan 2025 09:03:22 +0100 Subject: [PATCH 06/10] Rename dotarg to field This patch introduces a new function function.field to replace function.dotarg. The latter is kept as an alias. All use of dotarg is replaced by field. --- examples/platewithhole.py | 2 +- examples/poisson.py | 2 +- nutils/SI.py | 4 ++-- nutils/expression_v2.py | 8 +++---- nutils/function.py | 12 +++++++--- tests/test_function.py | 24 ++++++++++---------- tests/test_solver.py | 46 +++++++++++++++++++-------------------- 7 files changed, 52 insertions(+), 46 deletions(-) diff --git a/examples/platewithhole.py b/examples/platewithhole.py index 5f8845858..cda9bc8bf 100644 --- a/examples/platewithhole.py +++ b/examples/platewithhole.py @@ -75,7 +75,7 @@ def generate(self, radius): if self.nrefine: topo = topo.refine(self.nrefine) bsplinebasis = topo.basis('spline', degree=2) - sqr = topo.integral((function.dotarg('w', bsplinebasis) - weightfunc)**2, degree=9) + sqr = topo.integral((function.field('w', bsplinebasis) - weightfunc)**2, degree=9) controlweights = System(sqr, trial='w').solve()['w'] nurbsbasis = bsplinebasis * controlweights / weightfunc return topo.withboundary(hole='left', sym='top,bottom', far='right'), geom, nurbsbasis, 5 diff --git a/examples/poisson.py b/examples/poisson.py index 10081a1d0..be294675c 100644 --- a/examples/poisson.py +++ b/examples/poisson.py @@ -20,7 +20,7 @@ def main(nelems: int = 32): ''' topo, x = mesh.unitsquare(nelems, etype='square') - u = function.dotarg('u', topo.basis('std', degree=1)) + u = function.field('u', topo.basis('std', degree=1)) g = u.grad(x) J = function.J(x) diff --git a/nutils/SI.py b/nutils/SI.py index 284f7bdf7..245a67499 100644 --- a/nutils/SI.py +++ b/nutils/SI.py @@ -406,8 +406,8 @@ def __evaluate(op, *args, **kwargs): dims, args = zip(*Quantity.__unpack(*args)) return tuple(dim.wrap(ret) for (dim, ret) in zip(dims, op(*args, **kwargs))) - @register('dotarg') - def __dotarg(op, *args, **kwargs): + @register('field') + def __field(op, *args, **kwargs): dims, args = zip(*Quantity.__unpack(*args)) # we abuse the fact that unpack str returns dimensionless return functools.reduce(operator.mul, dims).wrap(op(*args, **kwargs)) diff --git a/nutils/expression_v2.py b/nutils/expression_v2.py index ade879738..be90be1d8 100644 --- a/nutils/expression_v2.py +++ b/nutils/expression_v2.py @@ -695,8 +695,8 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio >>> topo, ns.x = mesh.rectilinear([2, 2]) >>> ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) >>> ns.basis = topo.basis('spline', degree=1) - >>> ns.u = function.dotarg('u', ns.basis) - >>> ns.v = function.dotarg('v', ns.basis) + >>> ns.u = function.field('u', ns.basis) + >>> ns.v = function.field('v', ns.basis) >>> res = topo.integral('-∇_i(v) ∇_i(u) dV' @ ns, degree=2) >>> res += topo.boundary.integral('∇_i(v) u n_i dS' @ ns, degree=2) ''' @@ -718,7 +718,7 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio setattr(self, jacobian, function.jacobian(geom, numpy.size(geom) - i)) def add_field(self, __names: Union[str, Sequence[str]], *__bases, shape: Tuple[int, ...] = (), dtype: function.DType = float): - '''Add field(s) of the form ns.u = function.dotarg('u', ...) + '''Add field(s) of the form ns.u = function.field('u', ...) Parameters ---------- @@ -733,7 +733,7 @@ def add_field(self, __names: Union[str, Sequence[str]], *__bases, shape: Tuple[i ''' for name in (__names,) if isinstance(__names, str) else __names: - setattr(self, name, function.dotarg(name, *__bases, shape=shape, dtype=dtype)) + setattr(self, name, function.field(name, *__bases, shape=shape, dtype=dtype)) def copy_(self, **replacements: Mapping[str, function.Array]) -> 'Namespace': '''Return a copy of this namespace. diff --git a/nutils/function.py b/nutils/function.py index 90647c5f8..356386535 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -2368,8 +2368,14 @@ def rotmat(__arg: IntoArray) -> Array: return Array.cast(numpy.stack([trignormal(__arg), trigtangent(__arg)], 0)) +def dotarg(*args, **kwargs): + '''Alias for :func:`field`.''' + + return field(*args, **kwargs) + + @nutils_dispatch -def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtype: DType = float) -> Array: +def field(name: str, /, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtype: DType = float) -> Array: '''Return the inner product of the first axes of the given arrays with an argument with the given name. An argument with shape ``(arrays[0].shape[0], ..., arrays[-1].shape[0]) + @@ -2379,7 +2385,7 @@ def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtyp Parameters ---------- - argname : :class:`str` + name : :class:`str` The name of the argument. *arrays : :class:`Array` or something that can be :meth:`~Array.cast` into one The arrays to take inner products with. @@ -2394,7 +2400,7 @@ def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtyp The inner product with shape ``shape + arrays[0].shape[1:] + ... + arrays[-1].shape[1:]``. ''' - result = Argument(__argname, tuple(array.shape[0] for array in arrays) + tuple(shape), dtype=dtype) + result = Argument(name, tuple(array.shape[0] for array in arrays) + tuple(shape), dtype=dtype) for array in arrays: result = numpy.sum(_append_axes(result.transpose((*range(1, result.ndim), 0)), array.shape[1:]) * array, result.ndim-1) return result diff --git a/tests/test_function.py b/tests/test_function.py index 47846b813..d3faa8750 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -794,7 +794,7 @@ def test_nonempty_spaces(self): function.replace_arguments(function.Argument('foo', (2,), dtype=float), dict(foo=geom)) -class dotarg(TestCase): +class field(TestCase): def assertEvalAlmostEqual(self, f_actual, desired, **arguments): self.assertEqual(f_actual.shape, desired.shape) @@ -810,13 +810,13 @@ def test(self): a243 = numpy.arange(24, dtype=float).reshape(2, 4, 3) a4 = numpy.arange(4, dtype=float) a45 = numpy.arange(20, dtype=float).reshape(4, 5) - self.assertEvalAlmostEqual(function.dotarg('arg'), a, arg=a) - self.assertEvalAlmostEqual(function.dotarg('arg', shape=(2, 3)), a23, arg=a23) - self.assertEvalAlmostEqual(function.dotarg('arg', a4), numpy.einsum('i,i->', a4, a4), arg=a4) - self.assertEvalAlmostEqual(function.dotarg('arg', a45), numpy.einsum('i,ij->j', a4, a45), arg=a4) - self.assertEvalAlmostEqual(function.dotarg('arg', a24, shape=(3,)), numpy.einsum('ij,ik->jk', a23, a24), arg=a23) - self.assertEvalAlmostEqual(function.dotarg('arg', a2, a45), numpy.einsum('ij,i,jk->k', a24, a2, a45), arg=a24) - self.assertEvalAlmostEqual(function.dotarg('arg', a2, a45, shape=(3,)), numpy.einsum('ijk,i,jl->kl', a243, a2, a45), arg=a243) + self.assertEvalAlmostEqual(function.field('arg'), a, arg=a) + self.assertEvalAlmostEqual(function.field('arg', shape=(2, 3)), a23, arg=a23) + self.assertEvalAlmostEqual(function.field('arg', a4), numpy.einsum('i,i->', a4, a4), arg=a4) + self.assertEvalAlmostEqual(function.field('arg', a45), numpy.einsum('i,ij->j', a4, a45), arg=a4) + self.assertEvalAlmostEqual(function.field('arg', a24, shape=(3,)), numpy.einsum('ij,ik->jk', a23, a24), arg=a23) + self.assertEvalAlmostEqual(function.field('arg', a2, a45), numpy.einsum('ij,i,jk->k', a24, a2, a45), arg=a24) + self.assertEvalAlmostEqual(function.field('arg', a2, a45, shape=(3,)), numpy.einsum('ijk,i,jl->kl', a243, a2, a45), arg=a243) class jacobian(TestCase): @@ -1388,13 +1388,13 @@ def test_stokes(self): class Eval(TestCase): def test_single(self): - f = function.dotarg('v', numpy.array([1, 2, 3])) + f = function.field('v', numpy.array([1, 2, 3])) retval = function.eval(f, v=numpy.array([4, 5, 6])) self.assertEqual(retval, 4+10+18) def test_multiple(self): - f = function.dotarg('v', numpy.array([1, 2, 3])) - g = function.dotarg('v', numpy.array([3, 2, 1])) + f = function.field('v', numpy.array([1, 2, 3])) + g = function.field('v', numpy.array([3, 2, 1])) retvals = function.eval([f, g], v=numpy.array([4, 5, 6])) self.assertEqual(retvals, (4+10+18, 12+10+6)) @@ -1433,7 +1433,7 @@ class factor(TestCase): def test_lower_with_points(self): topo, geom = mesh.rectilinear([3]) - f = function.dotarg('dof') + f = function.field('dof') v = topo.sample('uniform', 1).eval(function.factor(f**2) * geom[0], dof=2.) self.assertAllAlmostEqual(v, [2, 6, 10]) diff --git a/tests/test_solver.py b/tests/test_solver.py index 13f7b5822..78fe72630 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -101,8 +101,8 @@ def setUp(self): super().setUp() domain, geom = mesh.rectilinear([8, 8]) basis = domain.basis('std', degree=1) - u = function.dotarg('u', basis) - v = function.dotarg('v', basis) + u = function.field('u', basis) + v = function.field('v', basis) sqr = domain.boundary['left'].integral(u**2, degree=2) self.cons = solver.optimize('u,', sqr) self.residual = domain.integral((v.grad(geom) @ u.grad(geom))*function.J(geom), degree=2) \ @@ -205,9 +205,9 @@ def setUp(self): domain, geom = mesh.rectilinear([numpy.linspace(0, 1, 9)] * 2) ubasis = domain.basis('std', degree=2) if self.vector: - u = function.dotarg('dofs', ubasis.vector(2)) + u = function.field('dofs', ubasis.vector(2)) else: - u = function.dotarg('dofs', ubasis, shape=(2,)) + u = function.field('dofs', ubasis, shape=(2,)) Geom = geom * [1.1, 1] + u self.cons = solver.optimize('dofs', domain.boundary['left,right'].integral((u**2).sum(0), degree=4), droptol=1e-15) self.boolcons = ~numpy.isnan(self.cons) @@ -258,7 +258,7 @@ def setUp(self): self.domain, geom = mesh.rectilinear([2, 2]) self.ubasis = self.domain.basis('std', degree=1) self.ns = Namespace() - self.ns.u = function.dotarg('dofs', self.ubasis) + self.ns.u = function.field('dofs', self.ubasis) def test_linear(self): err = self.domain.boundary['bottom'].integral('(u - 1)^2' @ self.ns, degree=2) @@ -337,8 +337,8 @@ class System(TestCase): def test_constant(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) - v = function.dotarg('v', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', 1)) + v = function.field('v', domain.basis('std', 1)) f = domain.integral(v * (1 + u) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') self.assertTrue(sys.is_linear) @@ -358,7 +358,7 @@ def test_constant(self): def test_constant_symmetric(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', 1)) f = domain.integral((1 + u - u**2) * function.J(geom), degree=2) sys = solver.System(f, trial='u') self.assertTrue(sys.is_linear) @@ -379,9 +379,9 @@ def test_constant_symmetric(self): def test_constant_matrix(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) - v = function.dotarg('v', domain.basis('std', 1)) - t = function.dotarg('t') + u = function.field('u', domain.basis('std', 1)) + v = function.field('v', domain.basis('std', 1)) + t = function.field('t') f = domain.integral(v * (t + u) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') self.assertTrue(sys.is_linear) @@ -401,9 +401,9 @@ def test_constant_matrix(self): def test_linear(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) - v = function.dotarg('v', domain.basis('std', 1)) - t = function.dotarg('t') + u = function.field('u', domain.basis('std', 1)) + v = function.field('v', domain.basis('std', 1)) + t = function.field('t') f = domain.integral(v * (u * t + 1) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') self.assertTrue(sys.is_linear) @@ -419,9 +419,9 @@ def test_linear(self): def test_nonlinear(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) - v = function.dotarg('v', domain.basis('std', 1)) - t = function.dotarg('t') + u = function.field('u', domain.basis('std', 1)) + v = function.field('v', domain.basis('std', 1)) + t = function.field('t') f = domain.integral(v * (u**2 + t) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') self.assertFalse(sys.is_linear) @@ -437,7 +437,7 @@ def test_nonlinear(self): def test_nonlinear_symmetric(self): domain, geom = mesh.rectilinear([9]) - u = function.dotarg('u', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', 1)) f = domain.integral(numpy.exp(u) * function.J(geom), degree=2) sys = solver.System(f, trial='u') self.assertFalse(sys.is_linear) @@ -457,7 +457,7 @@ class system_finitestrain(TestCase): def setUp(self): super().setUp() domain, geom = mesh.rectilinear([numpy.linspace(0, 1, 9)] * 2) - u = function.dotarg('u', domain.basis('std', degree=2), shape=(2,)) + u = function.field('u', domain.basis('std', degree=2), shape=(2,)) Geom = geom * [1.1, 1] + u self.cons = solver.optimize('u,', domain.boundary['left,right'].integral((u**2).sum(0), degree=4), droptol=1e-15) strain = .5 * (function.outer(Geom.grad(geom), axis=1).sum(0) - function.eye(2)) @@ -497,10 +497,10 @@ def setUp(self): dx = function.J(geom) ubasis = domain.basis('std', degree=2) pbasis = domain.basis('std', degree=1) - u = function.dotarg('u', ubasis, shape=(2,)) - v = function.dotarg('v', ubasis, shape=(2,)) - p = function.dotarg('p', pbasis) - q = function.dotarg('q', pbasis) + u = function.field('u', ubasis, shape=(2,)) + v = function.field('v', ubasis, shape=(2,)) + p = function.field('p', pbasis) + q = function.field('q', pbasis) self.cons = solver.optimize('u,', domain.boundary['top,bottom'].integral((u**2).sum(), degree=4), droptol=1e-10) self.cons = solver.optimize('u,', domain.boundary['left'].integral((u[0]-uin)**2 + u[1]**2, degree=4), droptol=1e-10, constrain=self.cons) res = gauss.integral((viscosity * numpy.einsum('ij,ij', v.grad(geom), u.grad(geom) + u.grad(geom).T) - v.div(geom) * p + q * u.div(geom)) * dx) From 97198c56e91283f18defafd8c73459fce4db54fd Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 24 Jan 2025 14:44:50 +0100 Subject: [PATCH 07/10] Reorganize basis specialization This patch changes the definition of Topology.basis, making it the sole entry point for basis creation and dispatching to the _tensorial_bases generator. Additionally, the degree argument is made keyword-only, with a deprecation warning on postional use, and 'name' argument is renamed (via deprecation) to 'btype' to avoid future clashes with the field name argument, while still allowing positional use. --- nutils/mesh.py | 2 +- nutils/topology.py | 99 ++++++++++++++++++++++++++++-------------- tests/test_function.py | 2 +- tests/test_solver.py | 20 ++++----- tests/test_topology.py | 6 +-- 5 files changed, 81 insertions(+), 48 deletions(-) diff --git a/nutils/mesh.py b/nutils/mesh.py index cb85a8e55..78470dac7 100644 --- a/nutils/mesh.py +++ b/nutils/mesh.py @@ -835,7 +835,7 @@ def unitcircle(nelems: int, variant: str) -> Tuple[Topology, function.Array]: # denominator onto a second order basis for efficient evaluation and # correct boundary gradients at the patch interfaces. - basis = topo.basis('spline', 2) + basis = topo.basis('spline', degree=2) # Minimize e = (f - B c)^2 = f^2 - 2 f B c + cT BT B c -> BT B c = BT f BB, BX, BY, BW, XX, YY, WW = topo.integrate([basis[:,None] * basis[None,:], basis * X, basis * Y, basis * W, X**2, Y**2, W**2], degree=4) diff --git a/nutils/topology.py b/nutils/topology.py index e6c7276c8..abbb2b93c 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -341,19 +341,56 @@ def f_coords(self) -> function.Array: raise NotImplementedError - def basis(self, name: str, *args, **kwargs) -> function.Basis: - ''' - Create a basis. + def _tensorial_bases(self, btype: str, **kwargs) -> function.Basis: + '''Generate bases for different tensorial directions as appropriate.''' + + # This default implementation generates a single basis through dispatch + # to the basis_ method, after stripping any hierarchical prefix + # to allow the same basis type to be used on the base topology as on + # its refinements. + + if not self.ndims: + return + + if btype.startswith('h-'): + btype = btype[2:] + elif btype.startswith('th-'): + btype = btype[3:] + kwargs.pop('truncation_tolerance', None) + + f = getattr(self, 'basis_' + btype) + yield f(**kwargs) + + @log.withcontext + def basis(self, btype=None, __degree=None, **kwargs): + '''Create a basis. + + A basis is a set of basis functions organized in a one dimensional + array. The available types varies depending on the topology. ''' + + # The basis method is not redefined in Topology subclasses; + # specialization happens via the _tensorial_bases method. + + if btype is None: + if 'name' not in kwargs: + raise TypeError("basis() missing 1 required positional argument: 'btype'") + btype = kwargs.pop('name') + warnings.deprecation("the basis type parameter has been renamed from 'name' to 'btype', please update your code") + + if __degree is not None: + if 'degree' in kwargs: + raise TypeError("basis() got multiple values for argument 'degree'") + kwargs['degree'] = __degree + warnings.deprecation('the degree parameter must be specified as a keyword argument, please update your code') + if self.ndims == 0: return function.PlainBasis([[[1]]], [[0]], 1, self.f_index, self.f_coords) - split = name.split('-', 1) - if len(split) == 2 and split[0] in ('h', 'th'): - name = split[1] # default to non-hierarchical bases - if split[0] == 'th': - kwargs.pop('truncation_tolerance', None) - f = getattr(self, 'basis_' + name) - return f(*args, **kwargs) + + basis, *other_bases = self._tensorial_bases(btype, **kwargs) + for other_basis in other_bases: + basis = numpy.ravel(basis[..., numpy.newaxis] * other_basis) + return basis def sample(self, ischeme: str, degree: int) -> Sample: 'Create sample.' @@ -1265,8 +1302,7 @@ def interfaces_spaces_unchecked(self, spaces: FrozenSet[str]) -> Topology: interfaces.append(self.topo1.interfaces_spaces(spaces1) * self.topo2) return Topology.disjoint_union(*interfaces) - def basis(self, name: str, degree: Union[int, Sequence[int]], **kwargs) -> function.Basis: - kwargs['degree'] = degree + def _tensorial_bases(self, btype, **kwargs): kwargs1 = {} kwargs2 = {} @@ -1308,10 +1344,8 @@ def basis(self, name: str, degree: Union[int, Sequence[int]], **kwargs) -> funct kwargs1.update(kwargs) kwargs2.update(kwargs) - basis1 = self.topo1.basis(name, **kwargs1) - basis2 = self.topo2.basis(name, **kwargs2) - assert basis1.ndim == basis2.ndim == 1 - return numpy.ravel(basis1[:,None] * basis2[None,:]) + yield from self.topo1._tensorial_bases(btype, **kwargs1) + yield from self.topo2._tensorial_bases(btype, **kwargs2) def sample(self, ischeme: str, degree: int) -> Sample: return self.topo1.sample(ischeme, degree) * self.topo2.sample(ischeme, degree) @@ -1389,8 +1423,8 @@ def f_coords(self) -> function.Array: def connectivity(self) -> Sequence[Sequence[int]]: return self.parent.connectivity - def basis(self, name: str, *args, **kwargs) -> function.Basis: - return self.parent.basis(name, *args, **kwargs) + def _tensorial_bases(self, btype, **kwargs) -> function.Basis: + yield from self.parent._tensorial_bases(btype, **kwargs) def sample(self, ischeme: str, degree: int) -> Sample: return self.parent.sample(ischeme, degree) @@ -1851,8 +1885,8 @@ def points(self): topo = topo.basetopo return UnionTopology(ptopos, pnames) - def basis(self, name, *args, **kwargs): - return self.basetopo.basis(name, *args, **kwargs) + def _tensorial_bases(self, btype, **kwargs): + yield from self.basetopo._tensorial_bases(btype, **kwargs) @cached_property def refined(self): @@ -2685,12 +2719,11 @@ def interfaces(self): irefs[iielem] = ref return SubsetTopology(baseinterfaces, irefs) - @log.withcontext - def basis(self, name, *args, **kwargs): + def _tensorial_bases(self, btype, **kwargs): if isinstance(self.basetopo, HierarchicalTopology): warnings.warn('basis may be linearly dependent; a linearly indepent basis is obtained by trimming first, then creating hierarchical refinements') - basis = self.basetopo.basis(name, *args, **kwargs) - return function.PrunedBasis(basis, self._indices, self.f_index, self.f_coords) + basis = self.basetopo.basis(btype, **kwargs) + yield function.PrunedBasis(basis, self._indices, self.f_index, self.f_coords) def locate(self, geom, coords, *, eps=0, skip_missing=False, **kwargs): sample = self.basetopo.locate(geom, coords, eps=eps, skip_missing=skip_missing, **kwargs) @@ -2903,8 +2936,7 @@ def interfaces(self): hopposites.append(level.interfaces.opposites[selection]) return TransformChainsTopology(self.space, hreferences, transformseq.chain(htransforms, self.transforms.todims, self.ndims-1), transformseq.chain(hopposites, self.transforms.todims, self.ndims-1)) - @log.withcontext - def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs): + def _tensorial_bases(self, btype, truncation_tolerance=1e-15, **kwargs): '''Create hierarchical basis. A hierarchical basis is constructed from bases on different levels of @@ -2945,14 +2977,15 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs): basis : :class:`nutils.function.Array` ''' - if name.startswith('h-'): + if btype.startswith('h-'): truncated = False - name = name[2:] - elif name.startswith('th-'): + btype = btype[2:] + elif btype.startswith('th-'): truncated = True - name = name[3:] + btype = btype[3:] else: - return super().basis(name, *args, **kwargs) + yield from super()._tensorial_bases(btype, **kwargs) + return # 1. identify active (supported) and passive (unsupported) basis functions ubases = [] @@ -2971,7 +3004,7 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs): prev_ielems = ielems_i = numpy.unique(numpy.concatenate([numpy.asarray(touchielems_i, dtype=int), nontouchielems_i], axis=0)) prev_transforms = topo.transforms - basis_i = topo.basis(name, *args, **kwargs) + basis_i = topo.basis(btype, **kwargs) assert isinstance(basis_i, function.Basis) ubases.insert(0, basis_i) # Basis functions that have at least one touchelem in their support. @@ -3050,7 +3083,7 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs): degree = poly.degree(self.ndims, max(c.shape[-1] for c in trans_coeffs)) hbasis_coeffs.append(numpy.concatenate([poly.change_degree(c, self.ndims, degree) for c in trans_coeffs], axis=0)) - return function.PlainBasis(hbasis_coeffs, hbasis_dofs, ndofs, self.f_index, self.f_coords) + yield function.PlainBasis(hbasis_coeffs, hbasis_dofs, ndofs, self.f_index, self.f_coords) class MultipatchTopology(TransformChainsTopology): diff --git a/tests/test_function.py b/tests/test_function.py index d3faa8750..2fb15eac9 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -1439,7 +1439,7 @@ def test_lower_with_points(self): def test_constant(self): topo, geom = mesh.rectilinear([3]) - basis = topo.basis('std', 1) + basis = topo.basis('std', degree=1) f = function.factor(topo.integral(basis[:, numpy.newaxis] * basis, degree=2)) ((i, j), v, shape), = function.evaluate(f, _post=sparse.extract) self.assertAllEqual(i, [0, 0, 1, 1, 1, 2, 2, 2, 3, 3]) diff --git a/tests/test_solver.py b/tests/test_solver.py index 78fe72630..546c5ac00 100644 --- a/tests/test_solver.py +++ b/tests/test_solver.py @@ -337,8 +337,8 @@ class System(TestCase): def test_constant(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) - v = function.field('v', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) + v = function.field('v', domain.basis('std', degree=1)) f = domain.integral(v * (1 + u) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') self.assertTrue(sys.is_linear) @@ -358,7 +358,7 @@ def test_constant(self): def test_constant_symmetric(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) f = domain.integral((1 + u - u**2) * function.J(geom), degree=2) sys = solver.System(f, trial='u') self.assertTrue(sys.is_linear) @@ -379,8 +379,8 @@ def test_constant_symmetric(self): def test_constant_matrix(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) - v = function.field('v', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) + v = function.field('v', domain.basis('std', degree=1)) t = function.field('t') f = domain.integral(v * (t + u) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') @@ -401,8 +401,8 @@ def test_constant_matrix(self): def test_linear(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) - v = function.field('v', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) + v = function.field('v', domain.basis('std', degree=1)) t = function.field('t') f = domain.integral(v * (u * t + 1) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') @@ -419,8 +419,8 @@ def test_linear(self): def test_nonlinear(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) - v = function.field('v', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) + v = function.field('v', domain.basis('std', degree=1)) t = function.field('t') f = domain.integral(v * (u**2 + t) * function.J(geom), degree=2) sys = solver.System(f, trial='u', test='v') @@ -437,7 +437,7 @@ def test_nonlinear(self): def test_nonlinear_symmetric(self): domain, geom = mesh.rectilinear([9]) - u = function.field('u', domain.basis('std', 1)) + u = function.field('u', domain.basis('std', degree=1)) f = domain.integral(numpy.exp(u) * function.J(geom), degree=2) sys = solver.System(f, trial='u') self.assertFalse(sys.is_linear) diff --git a/tests/test_topology.py b/tests/test_topology.py index c1c9f428a..fb98db908 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1190,7 +1190,7 @@ class multipatch_misc(TestCase): def test_basis_merge_dofs(self): patches = (0, 4, 2, 5), (9, 1, 5, 3), (2, 5, 6, 7), (5, 3, 7, 8) domain, geom = mesh.multipatch(patches=patches, nelems=1) - basis = domain.basis('spline', 1) + basis = domain.basis('spline', degree=1) actual = domain.sample('bezier', 2).eval(basis) desired = numpy.zeros((16, 10)) desired[numpy.arange(16), [0, 1, 2, 3, 4, 5, 3, 6, 2, 3, 7, 8, 3, 6, 8, 9]] = 1 @@ -1207,7 +1207,7 @@ def test_ring(self): domain, geom = mesh.multipatch(patches=patches, patchverts=patchverts, nelems=1) self.assertEqual(domain.connectivity.tolist(), [[1, -1, 2, -1], [-1, 0, 2, -1], [1, -1, -1, 0]]) self.assertEqual(len(domain.interfaces['interpatch']), 3) - self.assertEqual(len(domain.basis('std', 1)), 7) + self.assertEqual(len(domain.basis('std', degree=1)), 7) self.assertAlmostEqual(domain.integrate(function.J(geom), degree=1), 1) def test_partial_ring(self): @@ -1221,7 +1221,7 @@ def test_partial_ring(self): domain, geom = mesh.multipatch(patches=patches, nelems=1) self.assertEqual(len(domain.interfaces['interpatch']), 5) # 5th interface via 15-12-5-1 self.assertEqual(len(domain.boundary), 14) - self.assertEqual(len(domain.basis('std', 1)), 16) + self.assertEqual(len(domain.basis('std', degree=1)), 16) # this test fails if MultipatchTopology.basis_std delegates to # TransformChainsTopology.basis_std rather than MultipatchTopology.basis_spline. From 649029efa7a274ef28b6b3c2383fce52924fc7b9 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 22 Jan 2025 15:10:29 +0100 Subject: [PATCH 08/10] Add field method This patch introduces the Topology.field method, and modifies all examples to avoid creating explicit basis objects where possible. --- examples/adaptivity.py | 16 ++++++++-------- examples/burgers.py | 10 ++++++---- examples/cahnhilliard.py | 13 ++++++------- examples/coil.py | 6 ++++-- examples/cylinderflow.py | 12 ++++++++---- examples/drivencavity.py | 38 +++++++++++++++++++++----------------- examples/elasticity.py | 7 ++++--- examples/finitestrain.py | 2 +- examples/laplace.py | 6 +++--- examples/platewithhole.py | 3 ++- examples/poisson.py | 2 +- examples/torsion.py | 5 ++--- nutils/topology.py | 17 +++++++++++++++++ 13 files changed, 83 insertions(+), 54 deletions(-) diff --git a/examples/adaptivity.py b/examples/adaptivity.py index 71d468e1c..b03d2caa5 100644 --- a/examples/adaptivity.py +++ b/examples/adaptivity.py @@ -53,18 +53,18 @@ def main(etype: str = 'square', if irefine: refdom = domain.refined - refbasis = refdom.basis(btype, degree=degree) - ns.add_field('vref', refbasis) - res = refdom.integral('∇_k(vref) ∇_k(u) dV' @ ns, degree=degree*2) - res -= refdom.boundary.integral('vref ∇_k(u) n_k dS' @ ns, degree=degree*2) - indicator = res.derivative('vref').eval(**args) - supp = refbasis.get_support(indicator**2 > numpy.mean(indicator**2)) - domain = domain.refined_by(refdom.transforms[supp]) + ns.refbasis = refdom.basis(btype=btype, degree=degree) + res = refdom.integral('∇_k(refbasis_n) ∇_k(u) dV' @ ns, degree=degree*2) + res -= refdom.boundary.integral('refbasis_n ∇_k(u) n_k dS' @ ns, degree=degree*2) + indicator = numpy.square(res.eval(**args)) + irefelems = ns.refbasis.get_support(indicator > indicator.mean()) + domain = domain.refined_by(refdom.transforms[irefelems]) ns = Namespace() ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 'v'), domain.basis(btype, degree=degree)) + ns.u = domain.field('u', btype=btype, degree=degree) + ns.v = domain.field('v', btype=btype, degree=degree) ns.uexact = exact ns.du = 'u - uexact' diff --git a/examples/burgers.py b/examples/burgers.py index 7e072efea..8f5839ff2 100644 --- a/examples/burgers.py +++ b/examples/burgers.py @@ -40,14 +40,16 @@ def main(nelems: int = 40, ns = Namespace() ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 'u0', 'v'), domain.basis(btype, degree=degree)) - ns.add_field(('t', 't0')) - ns.dudt = '(u - u0) / (t - t0)' + ns.u = domain.field('u', btype=btype, degree=degree) + ns.du = ns.u - function.replace_arguments(ns.u, 'u:u0') + ns.v = domain.field('v', btype=btype, degree=degree) + ns.t = function.field('t') + ns.dt = ns.t - function.field('t0') ns.f = '.5 u^2' ns.C = 1 ns.uinit = 'exp(-25 x^2)' - res = domain.integral('(v dudt - ∇(v) f) dV' @ ns, degree=degree*2) + res = domain.integral('(v du / dt - ∇(v) f) dV' @ ns, degree=degree*2) res -= domain.interfaces.integral('[v] n ({f} - .5 C [u] n) dS' @ ns, degree=degree*2) sqr = domain.integral('(u - uinit)^2 dV' @ ns, degree=max(degree*2, 5)) diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index 9485acb22..5f218de5b 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -147,20 +147,19 @@ def main(size: Length = parse('10cm'), domain, geom = mesh.unitsquare(nelems, etype) ns = Namespace() - ns.x = size * geom + ns.x = geom * size ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('φ', 'φ0'), domain.basis('std', degree=degree)) - ns.add_field('η', domain.basis('std', degree=degree) * stens / epsilon) # basis scaling to give η the required unit - ns.add_field('dt') - ns.dt *= timestep + ns.φ = domain.field('φ', btype='std', degree=degree) + ns.dφ = ns.φ - function.replace_arguments(ns.φ, 'φ:φ0') + ns.η = domain.field('η', btype='std', degree=degree) * (stens / epsilon) + ns.dt = function.field('dt') * timestep ns.ε = epsilon ns.σ = stens ns.σmean = (wtensp + wtensn) / 2 ns.σdiff = (wtensp - wtensn) / 2 ns.σwall = 'σmean + φ σdiff' - ns.dφ = 'φ - φ0' ns.ψ = '.25 (φ^2 - 1)^2' - ns.δψ = '.25 dφ^2 (1 - φ0^2 / 6 - φ φ0 / 3 - φ^2 / 2)' if stable else '0' + ns.δψ = '.25 dφ^2 (1 - φ^2 + 2 φ dφ / 3 - dφ^2 / 6)' if stable else '0' ns.M = mobility ns.J_i = '-M ∇_i(η)' ns.v_i = 'φ J_i' diff --git a/examples/coil.py b/examples/coil.py index 07ae8ec94..6425c8ff2 100644 --- a/examples/coil.py +++ b/examples/coil.py @@ -105,7 +105,8 @@ def main(nelems: int = 50, X = RZ * REV ns.x = ns.rz @ ns.rot ns.define_for('x', gradient='∇', jacobians=('dV', 'dS'), curl='curl') - ns.add_field(('A', 'Atest'), RZ.basis('spline', degree=degree, removedofs=[[0, -1], [-1]])[:,numpy.newaxis] * ns.eθ, dtype=complex) + ns.A = RZ.field('A', btype='spline', degree=degree, removedofs=[[0, -1], [-1]], dtype=complex) * ns.eθ + ns.Atest = function.replace_arguments(ns.A, 'A:Atest') ns.B_i = 'curl_ij(A_j)' ns.E_i = '-j ω A_i' ns.Jind_i = 'σ E_i' @@ -123,7 +124,8 @@ def main(nelems: int = 50, # on a basis. ns.Borig = ns.B - ns.add_field(('B', 'Btest'), RZ.basis('spline', degree=degree), ns.rot, dtype=complex) + ns.B = function.field('B', RZ.basis('spline', degree=degree), ns.rot, dtype=complex) + ns.Btest = function.replace_arguments(ns.B, 'B:Btest') res = REV.integral(RZ.integral('Btest_i (B_i - Borig_i) dV' @ ns, degree=2*degree), degree=0) args = System(res, trial='B', test='Btest').solve(arguments=args) diff --git a/examples/cylinderflow.py b/examples/cylinderflow.py index 3a1f7ace1..f2fdd0bd7 100644 --- a/examples/cylinderflow.py +++ b/examples/cylinderflow.py @@ -115,11 +115,15 @@ def main(nelems: int = 99, ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) J = ns.x.grad(geom) detJ = numpy.linalg.det(J) - ns.add_field(('u', 'u0', 'v'), function.vectorize([ + ns.u = function.field('u', function.vectorize([ domain.basis('spline', degree=(degree, degree-1), removedofs=((0,), None)), domain.basis('spline', degree=(degree-1, degree))]) @ J.T / detJ) - ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1) / detJ) - ns.add_field(('t', 'dt')) + ns.p = domain.field('p', btype='spline', degree=degree-1) / detJ + ns.v = function.replace_arguments(ns.u, 'u:v') + ns.q = function.replace_arguments(ns.p, 'p:q') + ns.t = function.field('t') + ns.du = ns.u - function.replace_arguments(ns.u, 'u:u0') + ns.dt = function.field('dt') ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij' ns.ω = 'ε_ij ∇_j(u_i)' ns.N = 10 * degree / elemangle # Nitsche constant based on element size = elemangle/2 @@ -133,7 +137,7 @@ def main(nelems: int = 99, sqr = domain.integral('(.5 Σ_i (u_i - uinf_i)^2 - ∇_k(u_k) p) dV' @ ns, degree=degree*2) args = System(sqr, trial='u,p').solve(constrain=cons) # set initial condition to potential flow - res = domain.integral('v_i (u_i - u0_i) dV' @ ns, degree=degree*3) + res = domain.integral('v_i du_i dV' @ ns, degree=degree*3) res += domain.integral('(v_i ∇_j(u_i) u_j + ∇_j(v_i) σ_ij + q ∇_k(u_k)) dt dV' @ ns, degree=degree*3) res += domain.boundary['inner'].integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dt dS' @ ns, degree=degree*2) div = numpy.sqrt(abs(function.factor(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)))) # L2 norm of velocity divergence diff --git a/examples/drivencavity.py b/examples/drivencavity.py index c110924af..d1d28f697 100644 --- a/examples/drivencavity.py +++ b/examples/drivencavity.py @@ -100,15 +100,15 @@ def main(nelems: int = 32, ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) if not compatible: - ns.add_field(('u', 'v'), domain.basis('std', degree=degree), shape=(domain.ndims,)) - ns.add_field(('p', 'q'), domain.basis('std', degree=degree-1)[1:]) - ns.add_field('ψ', domain.basis('std', degree=2)[1:]) + ns.u = domain.field('u', btype='std', degree=degree, shape=[2]) + ns.p = domain.field('p', btype='std', degree=degree-1) + ns.ψ = domain.field('ψ', btype='std', degree=2) else: - ns.add_field(('u', 'v'), function.vectorize([ - domain.basis('spline', degree=(degree, degree-1)), - domain.basis('spline', degree=(degree-1, degree))])) - ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1)[1:]) - ns.add_field('ψ', domain.basis('spline', degree=degree)[1:]) + ns.u = function.field('u', function.vectorize([domain.basis('spline', degree=p) for p in degree - 1 + numpy.eye(2, dtype=int)])) + ns.p = domain.field('p', btype='spline', degree=degree-1) + ns.ψ = domain.field('ψ', btype='spline', degree=degree) + ns.v = function.replace_arguments(ns.u, 'u:v') + ns.q = function.replace_arguments(ns.p, 'p:q') ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij' # weak formulation for Stokes flow, over-integrating for improved @@ -119,6 +119,8 @@ def main(nelems: int = 32, # strong enforcement of non-penetrating boundary conditions sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2) cons = System(sqr, trial='u').solve_constraints(droptol=1e-15) + cons['p'] = numpy.zeros(res.argshapes['p'], dtype=bool) + cons['p'].flat[0] = True # point constraint if strongbc: # strong enforcement of tangential boundary conditions @@ -156,7 +158,9 @@ def postprocess(domain, ns, **arguments): # reconstruct velocity streamlines sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4) - arguments = System(sqr, trial='ψ').solve(arguments=arguments) + consψ = numpy.zeros(sqr.argshapes['ψ'], dtype=bool) + consψ.flat[0] = True # point constraint + arguments = System(sqr, trial='ψ').solve(arguments=arguments, constrain={'ψ': consψ}) bezier = domain.sample('bezier', 9) x, u, p, ψ = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'ψ'] @ ns, **arguments) @@ -197,7 +201,7 @@ def test_baseline(self): YcDGBdLMJR3Y/X+zdkhHHrEHM6lENt25+OU8OUi8PUn+klm87lacqiN4uQrZ4tUCLh3g4AFrV6Q/uctG gQ==''') with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' + self.assertAlmostEqual64(args0['p'][1:], ''' eNoBHgDh/+vVsNEXy6jTbiq1z7Av9C0mLJkw1NDTLEEtEC/xNAwED0s=''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' @@ -206,7 +210,7 @@ def test_baseline(self): nM/5nf5xevqZxDOq5w4bCwLlOoD6XDV/n1t//s5ZvjPzTjmdDjx55+Slky/MGaDgHFB3vz4DgynQfS9O A3WcBIkCAB7aSkk=''') with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' + self.assertAlmostEqual64(args1['p'][1:], ''' eNpz01W9oHVmuU7SJYtzgherdcr0n59dfiZT11yP97yCGQDN0Azu''') def test_mixed(self): @@ -218,7 +222,7 @@ def test_mixed(self): n+Y8k3RG+Mwio99nuoHqg4G48WzCmTignYUXDfXNzoedATuL4bMeA0Op9qczWqfXnTl2ioHhINAdHufv ntx18qoZSH7FSRAJAB13Sc0=''') with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' + self.assertAlmostEqual64(args0['p'][1:], ''' eNp7pKl+nf1KznmxS62ns/W+az/TNTL4ondU1/46t6GKKQDiJg1H''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' @@ -227,7 +231,7 @@ def test_mixed(self): n5xefpbrzEYjZ3MGhiogDr/YYbxbjYHhrH6lYcY55zNgZzGcNWBgUL0Uctr3zLzTt08xMOScZmCYdnbl qQMnpcxB8konQSQACVZG3A==''') with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' + self.assertAlmostEqual64(args1['p'][1:], ''' eNrbqjVZs1/ry/n48z1nSrW9L83RkTmneNZMO/TCOUNbMwDktQ3z''') def test_compatible(self): @@ -237,14 +241,14 @@ def test_compatible(self): eNpjYIAAvwvdBr9O2Zk90E8+rXQ6yxzGZ4CDTfr3z0H45hc2mjSagFgn9f1P15+G6Fc0PHSSgQEAx7kX 6A==''') with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' + self.assertAlmostEqual64(args0['p'][1:], ''' eNoL1u+7NOfUR929ugvORxlU6W7V1TcUuyiif/PKCf1yUwDfRw2t''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' eNpjYICA1HNRRkGnZ5r26m86bX3awvyBftS5C6dOmDHAwWxDmbMzTUEsrfMrTA6YgFjKV53OOJ0FsR7o F561OMnAAAC5tRfX''') with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' + self.assertAlmostEqual64(args1['p'][1:], ''' eNoz1VYx3HT6t16w/uKz73Uv6R7RNzx35swh7XdXrQ0TzADbMQ6l''') def test_strong(self): @@ -255,7 +259,7 @@ def test_strong(self): iNedYmDoPc3AsMGEgQGh77beyzPHzkqfYTpTcObp6Zmnlc7B5A5dOH4+9dyDMx807M50GvCdOnki8wRM DhcAAEYiNtQ=''') with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' + self.assertAlmostEqual64(args0['p'][1:], ''' eNoBHgDh/3fRlNYxy7PR0NVKz1ktVi1E1HowsdGJ07Qt/9PINA6QEUk=''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' @@ -263,7 +267,7 @@ def test_strong(self): /9BPPW1gXH6W2eSpkds5mBz72fvnUs/EnHt+etXpCWdZzgSd3W8Ck1M9L3zGGyi/4Pz80+ZnNp44c9L8 BEwOFwAA4RM3QA==''') with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' + self.assertAlmostEqual64(args1['p'][1:], ''' eNoBHgDh//ot489SzMEsntHezWTPAC+jL+XN/8wF02UxTc1JNhf0ELc=''') diff --git a/examples/elasticity.py b/examples/elasticity.py index 7cc82879a..f08d9dfa2 100644 --- a/examples/elasticity.py +++ b/examples/elasticity.py @@ -41,7 +41,7 @@ def main(nelems: int = 24, ns.δ = function.eye(domain.ndims) ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 't'), domain.basis(btype, degree=degree), shape=(2,)) + ns.u = domain.field('u', btype=btype, degree=degree, shape=[2]) ns.X_i = 'x_i + u_i' ns.λ = 1 ns.μ = .5/poisson - 1 @@ -61,8 +61,9 @@ def main(nelems: int = 24, if direct: ns.t_i = 'σ_ij n_j' # <-- this is an inadmissible boundary term else: - energy -= domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2) - args = System(energy, trial='t', test='u').solve(constrain={'t': numpy.isnan(cons['u'])}, arguments=args) + ns.t = domain.field('t', btype=btype, degree=degree, shape=[2]) + system = System(energy - domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2), trial='t', test='u') + args = system.solve(constrain={'t': numpy.isnan(cons['u'])}, arguments=args) F = domain.boundary['top'].integrate('t_i dS' @ ns, degree=degree*2, arguments=args) log.user('total clamping force:', F) diff --git a/examples/finitestrain.py b/examples/finitestrain.py index 058d5bd6b..5de553cab 100644 --- a/examples/finitestrain.py +++ b/examples/finitestrain.py @@ -51,7 +51,7 @@ def main(nelems: int = 20, ns.angle = angle * numpy.pi / 180 ns.λ = 2 * poisson ns.μ = 1 - 2 * poisson - ns.add_field('u', domain.basis(btype, degree=degree), shape=(domain.ndims,)) + ns.u = domain.field('u', btype=btype, degree=degree, shape=[domain.ndims]) ns.x_i = 'X_i + u_i' ns.ε_ij = '.5 (∇_j(u_i) + ∇_i(u_j))' ns.energy = 'λ ε_ii ε_jj + 2 μ ε_ij ε_ij' diff --git a/examples/laplace.py b/examples/laplace.py index 0c52ddfbc..783a2cd5e 100644 --- a/examples/laplace.py +++ b/examples/laplace.py @@ -49,13 +49,13 @@ def main(nelems: int = 10, # To be able to write index based tensor contractions, we need to bundle # all relevant functions together in a namespace. Here we add the geometry - # `x`, a test function `v`, and the solution `u`. The latter two are formed - # by contracting a basis with function arguments of the same name. + # `x`, a test function `v`, and the solution `u`. ns = Namespace() ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 'v'), domain.basis(btype, degree=degree)) + ns.u = domain.field('u', btype=btype, degree=degree) + ns.v = domain.field('v', btype=btype, degree=degree) # We are now ready to implement the Laplace equation. In weak form, the # solution is a scalar field `u` for which ∫_Ω ∇v·∇u - ∫_Γn v f = 0 ∀ v. diff --git a/examples/platewithhole.py b/examples/platewithhole.py index cda9bc8bf..ed86b0b32 100644 --- a/examples/platewithhole.py +++ b/examples/platewithhole.py @@ -121,7 +121,8 @@ def main(mode: Union[FCM, NURBS] = NURBS(), ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) ns.λ = 2 * poisson ns.μ = 1 - poisson - ns.add_field(('u', 'v'), basis, shape=(2,)) + ns.u = function.field('u', basis, shape=[2]) + ns.v = function.field('v', basis, shape=[2]) ns.X_i = 'x_i + u_i' ns.ε_ij = '(∇_j(u_i) + ∇_i(u_j)) / 2' ns.σ_ij = 'λ ε_kk δ_ij + 2 μ ε_ij' diff --git a/examples/poisson.py b/examples/poisson.py index be294675c..8c5dc4e8e 100644 --- a/examples/poisson.py +++ b/examples/poisson.py @@ -20,7 +20,7 @@ def main(nelems: int = 32): ''' topo, x = mesh.unitsquare(nelems, etype='square') - u = function.field('u', topo.basis('std', degree=1)) + u = topo.field('u', btype='std', degree=1) g = u.grad(x) J = function.J(x) diff --git a/examples/torsion.py b/examples/torsion.py index 9d0a84d9a..4826cc9f5 100644 --- a/examples/torsion.py +++ b/examples/torsion.py @@ -66,7 +66,7 @@ def main(length: float = 2*np.pi, zgrid = length * np.linspace(-.5, .5, round(length / elemsize)+1) θgrid = np.linspace(-np.pi, np.pi, round(2 * np.pi / elemsize)+1) cylinder, (z, θ) = mesh.rectilinear([zgrid, θgrid], periodic=(1,)) - φ = θ - (z / length * np.pi / 180) * function.Argument('φ', shape=()) + φ = θ - (z / length * np.pi / 180) * function.field('φ') if trim: cylinder = cylinder.trim(θ**2 + z**2 - trim**2, maxrefine=2) extrusion, r = mesh.line([1 - thickness/2, 1 + thickness/2], space='T') @@ -77,8 +77,7 @@ def main(length: float = 2*np.pi, ns.X = np.stack([z, r * np.sin(θ), r * np.cos(θ)]) # reference geometry ns.Xφ = np.stack([z * stretch, r * np.sin(φ), r * np.cos(φ)]) ns.define_for('X', gradient='∇', jacobians=('dV',)) - ns.add_field('u', topo.basis('spline', degree=degree, - removedofs=((0,-1),None,None)), shape=(3,)) # deformation, clamped + ns.u = topo.field('u', btype='spline', degree=degree, removedofs=((0,-1),None,None), shape=[3]) # deformation, clamped ns.x_i = 'Xφ_i + u_i' # deformed geometry ns.F_ij = '∇_j(x_i)' ns.J = np.linalg.det(ns.F) diff --git a/nutils/topology.py b/nutils/topology.py index abbb2b93c..2776f06ab 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -392,6 +392,23 @@ def basis(self, btype=None, __degree=None, **kwargs): basis = numpy.ravel(basis[..., numpy.newaxis] * other_basis) return basis + @log.withcontext + def field(self, name: str, /, *, shape=(), dtype=float, **kwargs): + '''Create a field. + + A field is a contraction of a basis with an argument of given type and + shape, the latter prefixed with the required dof axes. The following + two functions are fully equivalent on non-tensorial topologies, and + functionally equivalent on all topologies: + + ``` + self.field(name, shape=S, dtype=D, **B) + function.field(name, self.basis(**B), shape=S, dtype=D) + ``` + ''' + + return function.field(name, *self._tensorial_bases(**kwargs), shape=shape, dtype=dtype) + def sample(self, ischeme: str, degree: int) -> Sample: 'Create sample.' From 6bd65dad09e72a9212594de36a3c9c89d53a0fd0 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 3 Feb 2025 15:02:03 +0100 Subject: [PATCH 09/10] Introduce function.arguments_for This patch adds the function.arguments_for function, which can be used to retrieve array arguments. Unlike that .arguments attribute (which will be deprecated in the next commit) the returned dictionary maps to fully formed Argument objects, that can directly be used in operations such as derivative. Furthermore, the function is nutils-dispatched, so that if can be used even if the array is wrapped in an SI quantity. --- examples/cahnhilliard.py | 2 +- examples/drivencavity.py | 4 ++-- nutils/SI.py | 5 +++++ nutils/function.py | 24 ++++++++++++++++++++++++ tests/test_function.py | 27 +++++++++++++++++++++++++++ 5 files changed, 59 insertions(+), 3 deletions(-) diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index 5f218de5b..091f64c47 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -180,7 +180,7 @@ def main(size: Length = parse('10cm'), system = System(nrg / tol, trial='φ,η') numpy.random.seed(seed) - args = dict(φ=numpy.random.normal(0, .5, system.argshapes['φ'])) # initial condition + args = dict(φ=numpy.random.normal(0, .5, function.arguments_for(nrg)['φ'].shape)) # initial condition with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps: for istep in steps: diff --git a/examples/drivencavity.py b/examples/drivencavity.py index d1d28f697..cf70d42d2 100644 --- a/examples/drivencavity.py +++ b/examples/drivencavity.py @@ -119,7 +119,7 @@ def main(nelems: int = 32, # strong enforcement of non-penetrating boundary conditions sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2) cons = System(sqr, trial='u').solve_constraints(droptol=1e-15) - cons['p'] = numpy.zeros(res.argshapes['p'], dtype=bool) + cons['p'] = numpy.zeros(function.arguments_for(res)['p'].shape, dtype=bool) cons['p'].flat[0] = True # point constraint if strongbc: @@ -158,7 +158,7 @@ def postprocess(domain, ns, **arguments): # reconstruct velocity streamlines sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4) - consψ = numpy.zeros(sqr.argshapes['ψ'], dtype=bool) + consψ = numpy.zeros(function.arguments_for(sqr)['ψ'].shape, dtype=bool) consψ.flat[0] = True # point constraint arguments = System(sqr, trial='ψ').solve(arguments=arguments, constrain={'ψ': consψ}) diff --git a/nutils/SI.py b/nutils/SI.py index 245a67499..43ec2cc0c 100644 --- a/nutils/SI.py +++ b/nutils/SI.py @@ -411,6 +411,11 @@ def __field(op, *args, **kwargs): dims, args = zip(*Quantity.__unpack(*args)) # we abuse the fact that unpack str returns dimensionless return functools.reduce(operator.mul, dims).wrap(op(*args, **kwargs)) + @register('arguments_for') + def __attribute(op, *args, **kwargs): + __dims, args = zip(*Quantity.__unpack(*args)) + return op(*args, **kwargs) + del register ## DEFINE OPERATORS diff --git a/nutils/function.py b/nutils/function.py index 356386535..0ab24508e 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -2417,6 +2417,30 @@ def lower(self, args: LowerArgs) -> evaluable.Array: return evaluable.prependaxes(self._array, args.points_shape) +@nutils_dispatch +def arguments_for(*arrays) -> Dict[str, Argument]: + '''Get all arguments that array(s) depend on. + + Given any number of arrays, return a dictionary of all arguments involved, + mapping the name to the :class:`Argument` object. Raise a ``ValueError`` if + arrays have conflicting arguments, i.e. sharing a name but differing in + shape and/or dtype. + ''' + + arguments = {} + for array in arrays: + if isinstance(array, Array): + for name, (shape, dtype) in array._arguments.items(): + argument = arguments.get(name) + if argument is None: + arguments[name] = Argument(name, shape, dtype) + elif argument.shape != shape: + raise ValueError(f'inconsistent shapes for argument {name!r}') + elif argument.dtype != dtype: + raise ValueError(f'inconsistent dtypes for argument {name!r}') + return arguments + + # BASES diff --git a/tests/test_function.py b/tests/test_function.py index 2fb15eac9..8f211ad66 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -1450,3 +1450,30 @@ def test_lower_spaces(self): topo, geom = mesh.rectilinear([3]) with self.assertRaisesRegex(ValueError, 'cannot lower function with spaces \(.+\) - did you forget integral or sample?'): function.factor(geom) + + +class arguments_for(TestCase): + + def test_single(self): + x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int) + f = x**2 + self.assertEqual({a.name: (a.shape, a.dtype) for a in function.arguments_for(f).values()}, + {'x': ((3,2), int)}) + + def test_multiple(self): + x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int) + y = function.field('y', numpy.array([4,5]), dtype=float) + z = function.field('z', dtype=complex) + f = x * y + g = x**2 * z + self.assertEqual({a.name: (a.shape, a.dtype) for a in function.arguments_for(f, g).values()}, + {'x': ((3,2), int), 'y': ((2,), float), 'z': ((), complex)}) + + def test_conflict(self): + x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int) + y1 = function.field('y', numpy.array([4,5]), dtype=float) + y2 = function.field('y', dtype=complex) + f = x * y1 + g = x**2 * y2 + with self.assertRaisesRegex(ValueError, "inconsistent shapes for argument 'y'"): + function.arguments_for(f, g) From be8bc53549b338a8ef045610338873f7147756d8 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 4 Feb 2025 11:46:34 +0100 Subject: [PATCH 10/10] Deprecate Array.arguments and Array.argshapes This patch deprecates the .arguments and .argshapes attributes of the function.Array object in favour of the newly introduced function.arguments_for function. --- nutils/function.py | 64 ++++++++++++++++++++----------------- nutils/sample.py | 6 ++-- nutils/solver.py | 5 ++- tests/test_expression_v2.py | 10 +++--- tests/test_function.py | 4 +-- 5 files changed, 47 insertions(+), 42 deletions(-) diff --git a/nutils/function.py b/nutils/function.py index 0ab24508e..e4b8b9a0a 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -224,7 +224,7 @@ def __init__(self, shape: Shape, dtype: DType, spaces: FrozenSet[str], arguments self.shape = tuple(sh.__index__() for sh in shape) self.dtype = dtype self.spaces = frozenset(spaces) - self.arguments = dict(arguments) + self._arguments = dict(arguments) def lower(self, args: LowerArgs) -> evaluable.Array: raise NotImplementedError @@ -236,8 +236,8 @@ def as_evaluable_array(self) -> evaluable.Array: return self.lower(LowerArgs((), {}, {})) def __index__(self): - if self.arguments or self.spaces: - raise ValueError('cannot convert non-constant array to index: arguments={}'.format(','.join(self.arguments))) + if self._arguments or self.spaces: + raise ValueError('cannot convert non-constant array to index: arguments={}'.format(','.join(self._arguments))) elif self.ndim: raise ValueError('cannot convert non-scalar array to index: shape={}'.format(self.shape)) elif self.dtype != int: @@ -483,11 +483,17 @@ def replace(self, __arguments: Mapping[str, IntoArray]) -> 'Array': def contains(self, __name: str) -> bool: 'Test if target occurs in this function.' - return __name in self.arguments + return __name in self._arguments + + @property + def arguments(self) -> Mapping[str, Tuple[Shape, DType]]: + warnings.deprecation('array.arguments is deprecated and will be removed in Nutils 10, please use function.arguments_for(array) instead') + return self._arguments.copy() @property def argshapes(self) -> Mapping[str, Tuple[int, ...]]: - return {name: shape for name, (shape, dtype) in self.arguments.items()} + warnings.deprecation("array.argshapes[...] is deprecated and will be removed in Nutils 10, please use function.arguments_for(array)[...].shape instead") + return {name: shape for name, (shape, dtype) in self._arguments.items()} def conjugate(self): '''Return the complex conjugate, elementwise. @@ -684,7 +690,7 @@ def __init__(self, args: Iterable[Any], shape: Tuple[int], dtype: DType, npointw self._args = args self._npointwise = npointwise spaces = frozenset(space for arg in args if isinstance(arg, Array) for space in arg.spaces) - arguments = _join_arguments(arg.arguments for arg in args if isinstance(arg, Array)) + arguments = _join_arguments(arg._arguments for arg in args if isinstance(arg, Array)) super().__init__(shape=(*points_shape, *shape), dtype=dtype, spaces=spaces, arguments=arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -692,7 +698,7 @@ def lower(self, args: LowerArgs) -> evaluable.Array: add_points_shape = tuple(map(evaluable.asarray, self.shape[:self._npointwise])) points_shape = args.points_shape + add_points_shape coordinates = {space: evaluable.Transpose.to_end(evaluable.appendaxes(coords, add_points_shape), coords.ndim-1) for space, coords in args.coordinates.items()} - return _CustomEvaluable(type(self).__name__, self.evalf, self.partial_derivative, evalargs, self.shape[self._npointwise:], self.dtype, self.spaces, types.frozendict(self.arguments), LowerArgs(points_shape, types.frozendict(args.transform_chains), types.frozendict(coordinates))) + return _CustomEvaluable(type(self).__name__, self.evalf, self.partial_derivative, evalargs, self.shape[self._npointwise:], self.dtype, self.spaces, types.frozendict(self._arguments), LowerArgs(points_shape, types.frozendict(args.transform_chains), types.frozendict(coordinates))) @classmethod def evalf(cls, *args: Any) -> numpy.ndarray: @@ -833,7 +839,7 @@ class _WithoutPoints: def __init__(self, __arg: Array) -> None: self._arg = __arg self.spaces = __arg.spaces - self.arguments = __arg.arguments + self._arguments = __arg._arguments def lower(self, args: LowerArgs) -> evaluable.Array: return self._arg.lower(LowerArgs((), args.transform_chains, {})) @@ -851,7 +857,7 @@ def __init__(self, lower: Callable[..., evaluable.Array], *args: Lowerable, shap self._args = args assert all(hasattr(arg, 'lower') for arg in self._args) spaces = frozenset(space for arg in args for space in arg.spaces) - arguments = _join_arguments(arg.arguments for arg in self._args) + arguments = _join_arguments(arg._arguments for arg in self._args) super().__init__(shape, dtype, spaces, arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -916,8 +922,8 @@ def __init__(self, arg: Array, replacements: Dict[str, Array]) -> None: raise ValueError(f'replacement functions cannot be bound to a space, but replacement for Argument {old.name!r} is bound to {", ".join(new.spaces)}.') self._replacements[old.name] = new # Build arguments map with replacements. - unreplaced = {name: shape_dtype for name, shape_dtype in arg.arguments.items() if name not in replacements} - arguments = _join_arguments([unreplaced] + [replacement.arguments for replacement in self._replacements.values()]) + unreplaced = {name: shape_dtype for name, shape_dtype in arg._arguments.items() if name not in replacements} + arguments = _join_arguments([unreplaced] + [replacement._arguments for replacement in self._replacements.values()]) super().__init__(arg.shape, arg.dtype, arg.spaces, arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -950,7 +956,7 @@ def to_end(cls, array: Array, *axes: int) -> Array: def __init__(self, arg: Array, axes: Tuple[int, ...]) -> None: self._arg = arg self._axes = tuple(n.__index__() for n in axes) - super().__init__(tuple(arg.shape[axis] for axis in axes), arg.dtype, arg.spaces, arg.arguments) + super().__init__(tuple(arg.shape[axis] for axis in axes), arg.dtype, arg.spaces, arg._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: arg = self._arg.lower(args) @@ -964,7 +970,7 @@ class _Opposite(Array): def __init__(self, arg: Array, space: str) -> None: self._arg = arg self._space = space - super().__init__(arg.shape, arg.dtype, arg.spaces, arg.arguments) + super().__init__(arg.shape, arg.dtype, arg.spaces, arg._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: oppargs = LowerArgs(args.points_shape, dict(args.transform_chains), args.coordinates) @@ -1031,7 +1037,7 @@ def __init__(self, arg: Array, var: Argument) -> None: self._arg = arg self._var = var self._eval_var = evaluable.Argument(var.name, tuple(evaluable.constant(n) for n in var.shape), var.dtype) - arguments = _join_arguments((arg.arguments, var.arguments)) + arguments = _join_arguments((arg._arguments, var._arguments)) super().__init__(arg.shape+var.shape, complex if var.dtype == complex else arg.dtype, arg.spaces | var.spaces, arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -1056,7 +1062,7 @@ def __init__(self, func: Array, geom: Array) -> None: common_shape = broadcast_shapes(func.shape, geom.shape[:-1]) self._func = numpy.broadcast_to(func, common_shape) self._geom = numpy.broadcast_to(geom, (*common_shape, geom.shape[-1])) - arguments = _join_arguments((func.arguments, geom.arguments)) + arguments = _join_arguments((func._arguments, geom._arguments)) super().__init__(self._geom.shape, complex if func.dtype == complex else float, func.spaces | geom.spaces, arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -1082,7 +1088,7 @@ def __init__(self, func: Array, geom: Array) -> None: common_shape = broadcast_shapes(func.shape, geom.shape[:-1]) self._func = numpy.broadcast_to(func, common_shape) self._geom = numpy.broadcast_to(geom, (*common_shape, geom.shape[-1])) - arguments = _join_arguments((func.arguments, geom.arguments)) + arguments = _join_arguments((func._arguments, geom._arguments)) super().__init__(self._geom.shape, complex if func.dtype == complex else float, func.spaces | geom.spaces, arguments) def lower(self, args: LowerArgs) -> evaluable.Array: @@ -1112,7 +1118,7 @@ def __init__(self, geom: Array, tip_dim: Optional[int] = None) -> None: 'not greater than the dimension of the geometry.') self._tip_dim = tip_dim self._geom = geom - super().__init__((), float, geom.spaces, geom.arguments) + super().__init__((), float, geom.spaces, geom._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: geom = self._geom.lower(args) @@ -1133,7 +1139,7 @@ class _Normal(Array): def __init__(self, geom: Array) -> None: self._geom = geom assert geom.dtype == float - super().__init__(geom.shape, float, geom.spaces, geom.arguments) + super().__init__(geom.shape, float, geom.spaces, geom._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: geom = self._geom.lower(args) @@ -1168,7 +1174,7 @@ class _ExteriorNormal(Array): def __init__(self, rgrad: Array) -> None: assert rgrad.dtype == float and rgrad.shape[-2] == rgrad.shape[-1] + 1 self._rgrad = rgrad - super().__init__(rgrad.shape[:-1], float, rgrad.spaces, rgrad.arguments) + super().__init__(rgrad.shape[:-1], float, rgrad.spaces, rgrad._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: rgrad = self._rgrad.lower(args) @@ -1195,7 +1201,7 @@ def __init__(self, __arrays: Sequence[IntoArray], axis: int) -> None: shape=(*shape0[:self.axis], builtins.sum(array.shape[self.axis] for array in self.arrays), *shape0[self.axis+1:]), dtype=self.arrays[0].dtype, spaces=functools.reduce(operator.or_, (array.spaces for array in self.arrays)), - arguments=_join_arguments(array.arguments for array in self.arrays)) + arguments=_join_arguments(array._arguments for array in self.arrays)) def lower(self, args: LowerArgs) -> evaluable.Array: return util.sum(evaluable._inflate(array.lower(args), evaluable.Range(evaluable.constant(array.shape[self.axis])) + offset, evaluable.constant(self.shape[self.axis]), self.axis-self.ndim) @@ -1700,14 +1706,14 @@ def _argument_to_array(d: Any, array: Array) -> Iterable[Tuple[Argument, Array]] arg, new = item.split(':', 1) if isinstance(item, str) else item if isinstance(arg, str): - if arg not in array.arguments: + if arg not in array._arguments: continue - arg = Argument(arg, *array.arguments[arg]) + arg = Argument(arg, *array._arguments[arg]) elif not isinstance(arg, Argument): raise ValueError('Key must be string or argument') elif arg.name not in arguments: continue - elif array.arguments[arg.name] != (arg.shape, arg.dtype): + elif array._arguments[arg.name] != (arg.shape, arg.dtype): raise ValueError(f'Argument {arg.name!r} has wrong shape or dtype') if isinstance(new, str): @@ -1852,14 +1858,14 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array: arg = Array.cast(__arg) if isinstance(__var, str): - if __var not in arg.arguments: + if __var not in arg._arguments: raise ValueError('no such argument: {}'.format(__var)) - shape, dtype = arg.arguments[__var] + shape, dtype = arg._arguments[__var] __var = Argument(__var, shape, dtype=dtype) elif not isinstance(__var, Argument): raise ValueError('Expected an instance of `Argument` as second argument of `derivative` but got a `{}.{}`.'.format(type(__var).__module__, type(__var).__qualname__)) - if __var.name in arg.arguments: - shape, dtype = arg.arguments[__var.name] + if __var.name in arg._arguments: + shape, dtype = arg._arguments[__var.name] if __var.shape != shape: raise ValueError('Argument {!r} has shape {} in the function, but the derivative to {!r} with shape {} was requested.'.format(__var.name, shape, __var.name, __var.shape)) if __var.dtype != dtype: @@ -2411,7 +2417,7 @@ class factor(Array): def __init__(self, array: Array) -> None: self._array = evaluable.factor(array) - super().__init__(shape=array.shape, dtype=array.dtype, spaces=set(), arguments=array.arguments) + super().__init__(shape=array.shape, dtype=array.dtype, spaces=set(), arguments=array._arguments) def lower(self, args: LowerArgs) -> evaluable.Array: return evaluable.prependaxes(self._array, args.points_shape) @@ -2502,7 +2508,7 @@ def __init__(self, ndofs: int, nelems: int, index: Array, coords: Array) -> None self.nelems = nelems self.index = Array.cast(index, dtype=int, ndim=0) self.coords = coords - arguments = _join_arguments((index.arguments, coords.arguments)) + arguments = _join_arguments((index._arguments, coords._arguments)) super().__init__((ndofs,), float, spaces=index.spaces | coords.spaces, arguments=arguments) _index = evaluable.Argument('_index', shape=(), dtype=int) diff --git a/nutils/sample.py b/nutils/sample.py index ca992e8af..bd76ed8fb 100644 --- a/nutils/sample.py +++ b/nutils/sample.py @@ -913,7 +913,7 @@ class _Integral(function.Array): def __init__(self, integrand: function.Array, sample: Sample) -> None: self._integrand = integrand self._sample = sample - super().__init__(shape=integrand.shape, dtype=float if integrand.dtype in (bool, int) else integrand.dtype, spaces=integrand.spaces - frozenset(sample.spaces), arguments=integrand.arguments) + super().__init__(shape=integrand.shape, dtype=float if integrand.dtype in (bool, int) else integrand.dtype, spaces=integrand.spaces - frozenset(sample.spaces), arguments=integrand._arguments) def lower(self, args: function.LowerArgs) -> evaluable.Array: ielem = evaluable.loop_index('_sample_' + '_'.join(self._sample.spaces), self._sample.nelems) @@ -928,7 +928,7 @@ class _ConcatenatePoints(function.Array): def __init__(self, func: function.Array, sample: _TransformChainsSample) -> None: self._func = func self._sample = sample - super().__init__(shape=(sample.npoints, *func.shape), dtype=func.dtype, spaces=func.spaces - frozenset(sample.spaces), arguments=func.arguments) + super().__init__(shape=(sample.npoints, *func.shape), dtype=func.dtype, spaces=func.spaces - frozenset(sample.spaces), arguments=func._arguments) def lower(self, args: function.LowerArgs) -> evaluable.Array: axis = len(args.points_shape) @@ -948,7 +948,7 @@ def __init__(self, func: function.Array, indices: evaluable.Array) -> None: self._func = func self._indices = indices assert indices.ndim == 1 and func.shape[0] == indices.shape[0].__index__() - super().__init__(shape=func.shape, dtype=func.dtype, spaces=func.spaces, arguments=func.arguments) + super().__init__(shape=func.shape, dtype=func.dtype, spaces=func.spaces, arguments=func._arguments) def lower(self, args: function.LowerArgs) -> evaluable.Array: func = self._func.lower(args) diff --git a/nutils/solver.py b/nutils/solver.py index beafcb65e..f4fb54e3e 100644 --- a/nutils/solver.py +++ b/nutils/solver.py @@ -1002,9 +1002,8 @@ def _split_trial_test(target): def _target_helper(target, *args): trial, test = _split_trial_test(target) if test is not None: - arguments = function._join_arguments(arg.arguments for arg in args) - testargs = [function.Argument(t, *arguments[t]) for t in test] - args = [map(arg.derivative, testargs) for arg in args] + arguments = function.arguments_for(*args) + args = [[arg.derivative(arguments[t]) for t in test] for arg in args] elif len(args) > 1: shapes = [{f.shape for f in ziparg if f is not None} for ziparg in zip(*args)] if any(len(arg) != len(shapes) for arg in args) or any(len(shape) != 1 for shape in shapes): diff --git a/tests/test_expression_v2.py b/tests/test_expression_v2.py index 437170439..91f77345e 100644 --- a/tests/test_expression_v2.py +++ b/tests/test_expression_v2.py @@ -560,27 +560,27 @@ def test_define_for_3d(self): def test_add_single_field(self): ns = expression_v2.Namespace() ns.add_field('u', numpy.array([1,2,3])) - self.assertEqual(ns.u.argshapes, dict(u=(3,))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,))) self.assertEqual(ns.u.shape, ()) def test_add_multiple_fields(self): ns = expression_v2.Namespace() ns.add_field(('u', 'v'), numpy.array([1,2,3])) - self.assertEqual(ns.u.argshapes, dict(u=(3,))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,))) self.assertEqual(ns.u.shape, ()) - self.assertEqual(ns.v.argshapes, dict(v=(3,))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.v).values()}, dict(v=(3,))) self.assertEqual(ns.v.shape, ()) def test_add_single_field_multiple_bases(self): ns = expression_v2.Namespace() ns.add_field('u', numpy.array([1,2,3]), numpy.array([4,5,6,7])) - self.assertEqual(ns.u.argshapes, dict(u=(3,4))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,4))) self.assertEqual(ns.u.shape, ()) def test_add_single_field_with_shape(self): ns = expression_v2.Namespace() ns.add_field('u', numpy.array([1,2,3]), shape=(2,)) - self.assertEqual(ns.u.argshapes, dict(u=(3,2))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,2))) self.assertEqual(ns.u.shape, (2,)) def test_copy(self): diff --git a/tests/test_function.py b/tests/test_function.py index 8f211ad66..ba1bfb02c 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -128,7 +128,7 @@ def test_argshapes(self): a = function.Argument('a', (2, 3), dtype=int) b = function.Argument('b', (3,), dtype=int) f = (a * b[None]).sum(-1) - self.assertEqual(dict(f.argshapes), dict(a=(2, 3), b=(3,))) + self.assertEqual({a.name: a.shape for a in function.arguments_for(f).values()}, dict(a=(2, 3), b=(3,))) def test_argshapes_shape_mismatch(self): with self.assertRaises(Exception): @@ -425,7 +425,7 @@ def test(self): f = function._Unlower(e, frozenset(), arguments, function.LowerArgs((2, 3), {}, {})) self.assertEqual(f.shape, (4, 5)) self.assertEqual(f.dtype, int) - self.assertEqual(f.arguments, arguments) + self.assertEqual(f._arguments, arguments) self.assertEqual(f.lower(function.LowerArgs((2, 3), {}, {})), e) with self.assertRaises(ValueError): f.lower(function.LowerArgs((3, 4), {}, {}))