Backports WIP

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/__init__.py

@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9.0'
+__version__ = version = '9.0.1'
 version_name = 'jook-sing'

nutils/evaluable.py

@@ -501,6 +501,10 @@ class Array(Evaluable):
     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 @@ def dependencies(self):
         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 evalf(arr, symmetric, w_dtype, vt_dtype):
 
     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 @@ def evalf(dtype):
         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 @@ def _unaligned(self):
         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 @@ def _compile(self, builder):
         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(
@@ -7061,7 +7060,7 @@ def new_empty_array_for_evaluable(self, array: Array) -> _pyast.Variable:
         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:

nutils/function.py

@@ -2121,6 +2121,8 @@ def grad(arg: IntoArray, geom: IntoArray, /, ndims: int = 0, *, spaces: Optional
     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,19 @@ def assemble_block_csr(blocks):
     values = []
     rowptr = [0]
     colidx = []
+    dtype = None
     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'
+            if dtype is None:
+                dtype = block_values.dtype
+            elif dtype != block_values.dtype:
+                raise ValueError('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,17 +148,19 @@ def assemble_block_csr(blocks):
                     colidx.append(block_colidx[i:j])
                     ptr += j - i
                 rowptr.append(ptr)
-    return assemble_csr(numpy.concatenate(values), numpy.array(rowptr), numpy.concatenate(colidx), ncols)
+    if not values:
+        return empty((len(rowptr)-1, ncols), dtype)
+    return assemble_csr(numpy.concatenate(values, dtype=dtype), numpy.array(rowptr), numpy.concatenate(colidx), ncols)
 
 
 def fromsparse(data, inplace=False):
     (rowidx, colidx), values, (nrows, ncols) = sparse.extract(sparse.prune(sparse.dedup(data, inplace=inplace), inplace=True))
     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

@@ -1249,7 +1249,18 @@ def optimize(target, functional: evaluable.asarray, *, tol: float = 0., argument
     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 \

tests/test_matrix.py

@@ -73,10 +73,26 @@ def test_eye(self):
         self.assertEqual(ncols, 3)
 
     def test_assemble_block_csr(self):
-        A00 = [10., 20., 30.], [0, 2, 3], [0, 1, 0], 2 # 2x2
-        A01 = [], [0, 0, 0], [], 2 # 2x2
-        A10 = [40., 60.], [0, 1, 2], [0, 0], 1 # 2x1
-        A11 = [50.], [0, 1, 1], [2], 3 # 2x3
+        A00 = (
+            numpy.array([10., 20., 30.], dtype=float),
+            numpy.array([0, 2, 3], dtype=int),
+            numpy.array([0, 1, 0], dtype=int),
+            2) # 2x2
+        A01 = (
+            numpy.array([], dtype=float),
+            numpy.array([0, 0, 0], dtype=int),
+            numpy.array([], dtype=int),
+            2) # 2x2
+        A10 = (
+            numpy.array([40., 60.], dtype=float),
+            numpy.array([0, 1, 2], dtype=int),
+            numpy.array([0, 0], dtype=int),
+            1) # 2x1
+        A11 = (
+            numpy.array([50.], dtype=float),
+            numpy.array([0, 1, 1], dtype=int),
+            numpy.array([2], dtype=int),
+            3) # 2x3
         values, rowptr, colidx, ncols = matrix.assemble_block_csr([[A00, A01], [A10, A11]])
         self.assertEqual(values.tolist(), [10., 20., 30., 40., 50., 60.])
         self.assertEqual(rowptr.tolist(), [0, 2, 3, 5, 6])

tests/test_solver.py

@@ -271,9 +271,13 @@ def test_nanres(self):
         numpy.testing.assert_almost_equal(dofs, .75)
 
     def test_unknowntarget(self):
-        err = self.domain.integral('(sqrt(1 - u) - .5)^2' @ self.ns, degree=2)
-        with self.assertRaises(KeyError):
+        err = self.domain.integral('(u - .5)^2' @ self.ns, degree=2)
+        with self.assertRaises(ValueError):
             dofs = solver.optimize(['dofs', 'other'], err, tol=1e-10)
+        dofs = solver.optimize(['dofs', 'other'], err, tol=1e-10, droptol=1e-10)
+        self.assertEqual(list(dofs), ['dofs'])
+        dofs = solver.optimize(['other'], err, tol=1e-10, droptol=1e-10)
+        self.assertEqual(dofs, {})
         with self.assertRaises(KeyError):
             dofs = solver.optimize('other', err, tol=1e-10, droptol=1e-10)