Backported fixes from up to v10a3

examples/finitestrain.py

@@ -54,7 +54,7 @@ def main(nelems: int = 20,
     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'
+    ns.energy = '.5 λ ε_ii ε_jj + μ ε_ij ε_ij'
 
     sqr = domain.boundary['left'].integral('u_k u_k dS' @ ns, degree=degree*2)
     sqr += domain.boundary['right'].integral('((u_0 - X_1 sin(2 angle) - cos(angle) + 1)^2 + (u_1 - X_1 (cos(2 angle) - 1) + sin(angle))^2) dS' @ ns, degree=degree*2)

examples/turek.geo

@@ -41,7 +41,7 @@ _() = BooleanFragments{ Surface{3,4,5}; Point{A,B}; Delete; }{};
 
 bnd_cylinder() = Abs(Boundary{ Surface{3}; });
 bnd_structure() = Abs(Boundary{ Surface{4}; });
-tmp = bnd_structure;
+tmp = bnd_structure();
 bnd_structure -= bnd_cylinder();
 bnd_cylinder -= tmp();
 bnd_fluid() = Abs(Boundary{ Surface{5}; });

nutils/_util.py

@@ -208,15 +208,15 @@ def loadlib(name):
     if sys.platform == 'linux':
         libsubdir = 'lib'
         libname = f'lib{name}.so'
-        versioned = f'lib{name}.so.(\d+)'
+        versioned = f'lib{name}.so.(\\d+)'
     elif sys.platform == 'darwin':
         libsubdir = 'lib'
         libname = f'lib{name}.dylib'
-        versioned = f'lib{name}.(\d+).dylib'
+        versioned = f'lib{name}.(\\d+).dylib'
     elif sys.platform == 'win32':
         libsubdir = r'Library\bin'
         libname = f'{name}.dll'
-        versioned = f'{name}.(\d+).dll'
+        versioned = f'{name}.(\\d+).dll'
     else:
         return
 

nutils/evaluable.py

@@ -501,6 +501,10 @@
     def ndim(self):
         return len(self.shape)
 
+    @property
+    def ast_dtype(self):
+        return _pyast.LiteralStr({bool: 'bool', int: 'int64', float: 'float64', complex: 'complex128'}[self.dtype])
+
     def __getitem__(self, item):
         if not isinstance(item, tuple):
             item = item,
@@ -2917,7 +2921,7 @@
         return self.arg,
 
     def _compile_expression(self, arg):
-        return _pyast.Variable('numpy').get_attr('array').call(arg, dtype=_pyast.Variable(self.dtype.__name__))
+        return _pyast.Variable('numpy').get_attr('array').call(arg, dtype=self.ast_dtype)
 
     def _simplified(self):
         if iszero(self.arg):
@@ -3201,12 +3205,7 @@
 
     def _compile_expression(self, array):
         evalf = _pyast.Variable('evaluable').get_attr('Eig').get_attr('evalf')
-        return evalf.call(
-            array,
-            _pyast.LiteralBool(self.symmetric),
-            _pyast.Variable(self._w_dtype.__name__),
-            _pyast.Variable(self._vt_dtype.__name__),
-        )
+        return evalf.call(array, _pyast.LiteralBool(self.symmetric), self[0].ast_dtype, self[1].ast_dtype)
 
 
 class ArrayFromTuple(Array):
@@ -3263,7 +3262,7 @@
         return numpy.array(numpy.nan, dtype)
 
     def _compile_expression(self):
-        return _pyast.Variable('evaluable').get_attr('Singular').get_attr('evalf').call(_pyast.Variable(self.dtype.__name__))
+        return _pyast.Variable('evaluable').get_attr('Singular').get_attr('evalf').call(self.ast_dtype)
 
 
 class Zeros(Array):
@@ -3281,7 +3280,7 @@
         return Zeros((), self.dtype), ()
 
     def _compile_expression(self, *shape):
-        return _pyast.Variable('numpy').get_attr('zeros').call(_pyast.Tuple(shape), dtype=_pyast.Variable(self.dtype.__name__))
+        return _pyast.Variable('numpy').get_attr('zeros').call(_pyast.Tuple(shape), dtype=self.ast_dtype)
 
     def _node(self, cache, subgraph, times, unique_loop_ids):
         if self.ndim:
@@ -3933,7 +3932,7 @@
         shape = builder.compile(self.shape)
         out = builder.get_variable_for_evaluable(self)
         block = builder.get_block_for_evaluable(self)
-        block.assign_to(out, _pyast.Variable('numpy').get_attr('asarray').call(builder.get_argument(self.name), dtype=_pyast.Variable(self.dtype.__name__)))
+        block.assign_to(out, _pyast.Variable('numpy').get_attr('asarray').call(builder.get_argument(self.name), dtype=self.ast_dtype))
         block.if_(_pyast.BinOp(shape, '!=', out.get_attr('shape'))).raise_(
             _pyast.Variable('ValueError').call(
                 _pyast.LiteralStr('argument {!r} has the wrong shape: expected {}, got {}').get_attr('format').call(
@@ -6855,12 +6854,6 @@
             _pyast.If(first_run, main, main_rerun),
         ])
 
-    if compile_parallel:
-        main = _pyast.Block([
-            _pyast.Assign(_pyast.Variable('lock'), _pyast.Variable('multiprocessing').get_attr('Lock').call()),
-            main,
-        ])
-
     if stats == 'log':
         main = _pyast.Block([
             _pyast.Assign(_pyast.Variable('stats'), _pyast.Variable('collections').get_attr('defaultdict').call(_pyast.Variable('Stats'))),
@@ -6990,7 +6983,7 @@
         self._evaluable_block_map = evaluable_block_map
         self._evaluable_deps = evaluable_deps
         self._origin = origin
-        self._shared_arrays = set()
+        self._shared_arrays = {}
         self._new_eid = map('e{}'.format, new_index).__next__
         self.new_var = lambda: _pyast.Variable('v{}'.format(next(new_index)))
 
@@ -7056,12 +7049,15 @@
             # `array` might initiate an inplace add inside a loop, in which
             # case `out` must be a shared memory array or we'll only see
             # updates to `out` from the parent process.
-            self._shared_arrays.add(out)
+            assert out not in self._shared_arrays
+            lock = self.get_lock_for_evaluable(array)
+            self._shared_arrays[out] = lock
+            self._blocks[0,].append(_pyast.Assign(lock, _pyast.Variable('multiprocessing').get_attr('Lock').call()))
             py_alloc = _pyast.Variable('parallel').get_attr('shempty')
         else:
             py_alloc = _pyast.Variable('numpy').get_attr('empty')
         alloc_block = self.get_block_for_evaluable(array, block_id=alloc_block_id, comment='alloc')
-        alloc_block.assign_to(out, py_alloc.call(shape, dtype=_pyast.Variable(array.dtype.__name__)))
+        alloc_block.assign_to(out, py_alloc.call(shape, dtype=array.ast_dtype))
         return out, out_block_id
 
     def add_constant(self, value) -> _pyast.Variable:
@@ -7121,6 +7117,10 @@
         # somewhere else.
         return _pyast.Variable('v{}'.format(self._get_evaluable_index(evaluable)))
 
+    def get_lock_for_evaluable(self, evaluable: Evaluable) -> _pyast.Variable:
+        # Returns the lock `lock{id}` where `id` is the unique index of `evaluable`.
+        return _pyast.Variable('lock{}'.format(self._get_evaluable_index(evaluable)))
+
     def _get_evaluable_index(self, evaluable: Evaluable) -> int:
         # Assigns a unique index to `evaluable` if not already assigned and
         # returns the index.
@@ -7160,30 +7160,37 @@
         self._block = block
         self.new_var = parent.new_var
 
-    def _needs_lock(self, /, *args, **kwargs):
-        # Returns true if any of the arguments references variables that
-        # require a lock.
-        variables = frozenset().union(*(arg.variables for arg in args), *(arg.variables for arg in kwargs.values()))
-        return not self._parent._shared_arrays.isdisjoint(variables)
+    def _needs_lock(self, condition):
+        # Returns true if condition references variables that require a lock.
+        return any(var in self._parent._shared_arrays for var in condition.variables)
+
+    def _iter_locks(self, /, *args, **kwargs):
+        variables = dict.fromkeys(var for args_ in (args, kwargs.values()) for arg in args_ for var in arg.variables) # stable union
+        return filter(None, map(self._parent._shared_arrays.get, variables))
 
     def _block_for(self, /, *args, **kwargs):
         # If any of the arguments references variables that require a lock,
         # return a new block that is enclosed in a `with lock` block. Otherwise
         # simply return our block.
-        if self._needs_lock(*args, **kwargs):
-            block = _pyast.Block()
-            self._block.append(_pyast.With(_pyast.Variable('lock'), block))
-            return block
-        else:
-            return self._block
+        block = self._block
+        for lock in self._iter_locks(*args, **kwargs):
+            with_block = _pyast.Block()
+            block.append(_pyast.With(lock, with_block))
+            block = with_block
+        return block
 
     def exec(self, expression: _pyast.Expression):
         # Appends the statement `{expression}`. Returns nothing.
         self._block_for(expression).append(_pyast.Exec(expression))
 
     def assign_to(self, lhs: _pyast.Expression, rhs: _pyast.Expression) -> _pyast.Variable:
         # Appends the statement `{lhs} = {rhs}` and returns `lhs`.
-        self._block_for(lhs, rhs).append(_pyast.Assign(lhs, rhs))
+        if isinstance(lhs, _pyast.Variable):
+            # Assignments to variables (not sliced) never need locks.
+            block = self._block_for(rhs)
+        else:
+            block = self._block_for(lhs, rhs)
+        block.append(_pyast.Assign(lhs, rhs))
         return lhs
 
     def eval(self, expression: _pyast.Expression) -> _pyast.Variable:

nutils/export.py

@@ -218,8 +218,8 @@
 
     Args
     ----
-    name : :class:`str`
-      Destination file name (without vtk extension).
+    name : :class:`str` or file object
+      Destination file name (without vtk extension) or file object.
     tri : :class:`int` array
       Triangulation.
     x : :class:`float` array
@@ -270,8 +270,14 @@
     t_cells[:, 0] = nverts
     t_cells[:, 1:] = cells
 
-    name_vtk = name + '.vtk'
-    with log.userfile(name_vtk, 'wb') as vtk:
+    if isinstance(name, str):
+        opener = log.userfile(name + '.vtk', 'wb')
+    elif hasattr(name, 'write'):
+        opener = contextlib.nullcontext(name)
+    else:
+        raise ValueError(f'cannot write vtk data to {name!r}')
+
+    with opener as vtk:
         vtk.write(b'# vtk DataFile Version 3.0\nvtk output\nBINARY\nDATASET UNSTRUCTURED_GRID\n')
         vtk.write('POINTS {} {}\n'.format(npoints, vtkdtype[points.dtype]).encode('ascii'))
         points.tofile(vtk)

nutils/function.py

@@ -2121,6 +2121,8 @@
     geom = Array.cast(geom)
     if geom.dtype != float:
         raise ValueError('The geometry must be real-valued.')
+    if arg.dtype in (int, bool):
+        return zeros(arg.shape + geom.shape, float)
     if spaces is None:
         spaces = geom.spaces
     else:

nutils/matrix/__init__.py

@@ -119,14 +119,16 @@ def assemble_block_csr(blocks):
     values = []
     rowptr = [0]
     colidx = []
+    dtype = blocks[0][0][0].dtype
     for row in blocks:
         nrows = len(row[0][1]) - 1
         block_data = []
         col_offset = 0
         for block_values, block_rowptr, block_colidx, block_ncols in row:
             assert len(block_rowptr) - 1 == nrows, 'sparse blocks have inconsistent row sizes'
+            assert block_values.dtype == dtype, 'blocks must all have the same dtype'
             if len(block_values):
-                block_data.append((numpy.asarray(block_values), numpy.asarray(block_rowptr), numpy.add(block_colidx, col_offset)))
+                block_data.append((block_values, block_rowptr, block_colidx + col_offset))
             col_offset += block_ncols
         assert col_offset == ncols, 'sparse blocks have inconsistent column sizes'
         ptr = rowptr[-1]
@@ -143,6 +145,9 @@ def assemble_block_csr(blocks):
                     colidx.append(block_colidx[i:j])
                     ptr += j - i
                 rowptr.append(ptr)
+    if not values:
+        # Required shortcut because numpy.concatenate cannot deal with empty lists
+        return empty((len(rowptr)-1, ncols), dtype)
     return assemble_csr(numpy.concatenate(values), numpy.array(rowptr), numpy.concatenate(colidx), ncols)
 
 
@@ -151,9 +156,9 @@ def fromsparse(data, inplace=False):
     return assemble_coo(values, rowidx, nrows, colidx, ncols)
 
 
-def empty(shape):
+def empty(shape, dtype=float):
     nrows, ncols = shape
-    return assemble_csr(numpy.zeros(0, dtype=float), numpy.zeros(nrows+1, dtype=int), numpy.zeros(0, dtype=int), ncols)
+    return assemble_csr(numpy.zeros(0, dtype=dtype), numpy.zeros(nrows+1, dtype=int), numpy.zeros(0, dtype=int), ncols)
 
 
 def diag(d):

nutils/solver.py

@@ -264,6 +264,10 @@
     def __nutils_hash__(self):
         return types.nutils_hash(('System', self.trials, self.__value if self.is_symmetric else self.__block_residual))
 
+    def __reduce__(self):
+        initargs = self.__value if self.is_symmetric else self.__block_residual, self.trials
+        return System, initargs
+
     def deconstruct(self, arguments, constrain):
         arguments = arguments.copy()
         xparts = []
@@ -648,7 +652,7 @@
 
     def __call__(self, system, *, arguments: ArrayDict = {}, constrain: ArrayDict = {}) -> System.MethodIter:
         arguments, x = system.deconstruct(arguments, constrain)
-        linargs = _copy_with_defaults(self.linargs, rtol=1-3, symmetric=system.is_symmetric)
+        linargs = _copy_with_defaults(self.linargs, rtol=1e-3, symmetric=system.is_symmetric)
         while True:
             jac, res = system.assemble_jacobian_residual(arguments, x)
             yield system.construct(arguments, x), numpy.linalg.norm(res)
@@ -681,7 +685,7 @@
 
     def __call__(self, system, *, arguments: ArrayDict = {}, constrain: ArrayDict = {}) -> System.MethodIter:
         arguments, x = system.deconstruct(arguments, constrain)
-        linargs = _copy_with_defaults(self.linargs, rtol=1-3, symmetric=system.is_symmetric)
+        linargs = _copy_with_defaults(self.linargs, rtol=1e-3, symmetric=system.is_symmetric)
 
         res = system.assemble_residual(arguments, x)
         resnorm = numpy.linalg.norm(res)
@@ -741,7 +745,7 @@
 
     def __call__(self, system, *, arguments: ArrayDict = {}, constrain: ArrayDict = {}) -> System.MethodIter:
         arguments, x = system.deconstruct(arguments, constrain)
-        linargs = _copy_with_defaults(self.linargs, rtol=1-3, symmetric=system.is_symmetric)
+        linargs = _copy_with_defaults(self.linargs, rtol=1e-3, symmetric=system.is_symmetric)
         jac, res = system.assemble_jacobian_residual(arguments, x)
         relax = self.relax0
         while True: # newton iterations
@@ -789,7 +793,7 @@
             raise ValueError('problem is not symmetric')
         arguments, x = system.deconstruct(arguments, constrain)
         jac, res, val = system.assemble_jacobian_residual_value(arguments, x)
-        linargs = _copy_with_defaults(self.linargs, rtol=1-3, symmetric=system.is_symmetric)
+        linargs = _copy_with_defaults(self.linargs, rtol=1e-3, symmetric=system.is_symmetric)
         relax = 0.
         while True:
             yield system.construct(arguments, x), numpy.linalg.norm(res)
@@ -916,7 +920,7 @@
     def __call__(self, system, *, arguments: ArrayDict = {}, constrain: ArrayDict = {}) -> System.MethodIter:
         arguments, x = system.deconstruct(arguments, constrain)
         free = numpy.concatenate([numpy.isnan(arguments[t]).ravel() for t in system.trials])
-        linargs = _copy_with_defaults(self.linargs, rtol=1-3)
+        linargs = _copy_with_defaults(self.linargs, rtol=1e-3)
         djac = matrix.assemble_block_csr(evaluable.eval_once([[evaluable.as_csr(evaluable._flat(evaluable.derivative(evaluable._flat(res), arg).simplified, 2)) for arg in system.trial_args] for res in self.inertia], arguments=arguments)).submatrix(free, free)
 
         first = _First()
@@ -1249,7 +1253,18 @@
     linargs = _strip(kwargs, 'lin')
     if kwargs:
         raise TypeError('unexpected keyword arguments: {}'.format(', '.join(kwargs)))
-    system = System(functional, *_split_trial_test(target))
+    trial, test = _split_trial_test(target)
+    args = function.arguments_for(functional)
+    rm = {i for i, t in enumerate(trial) if t not in args}
+    if rm:
+        if not droptol:
+            raise ValueError(f'target {", ".join(trial[i] for i in rm)} does not occur in integrand; consider setting droptol>0')
+        trial = [t for i, t in enumerate(trial) if i not in rm]
+        if test is not None:
+            test = [t for i, t in enumerate(test) if i not in rm]
+    if not trial:
+        return {}
+    system = System(functional, trial, test)
     if droptol is not None:
         return system.solve_constraints(arguments=arguments, constrain=constrain or {}, linargs=linargs, droptol=droptol)
     method = Direct(**linargs) if system.is_linear \

nutils/topology.py

@@ -486,7 +486,7 @@
                 bfun = onto * fun
             elif len(onto.shape) == 2:
                 Afun = numpy.einsum('ik,jk', onto, onto)
-                bfun = function.sum(onto * fun, -1)
+                bfun = numpy.sum(onto * fun, -1)
                 if fun2.ndim:
                     fun2 = fun2.sum(-1)
             else:
@@ -512,8 +512,8 @@
                 ufun = onto * fun
                 afun = onto
             elif len(onto.shape) == 2:
-                ufun = function.sum(onto * fun, axis=-1)
-                afun = function.norm2(onto)
+                ufun = numpy.sum(onto * fun, axis=-1)
+                afun = numpy.linalg.norm(onto, axis=-1)
             else:
                 raise Exception
             J = function.J(geometry)