From b7d6a7e7a2590618b10a21a2d8eb27f1716d7a50 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Tue, 30 Nov 2021 14:20:46 +0100
Subject: [PATCH 01/11] add function.shape, ndim, size

This patch adds the methods function.shape, function.ndim and
function.size, as well as support for numpy.shape, numpy.ndim and
numpy.size via the array function protocol.
---
 nutils/function.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/nutils/function.py b/nutils/function.py
index d37f4fea6..fa72e2bab 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -1074,6 +1074,18 @@ def asarray(__arg: IntoArray) -> Array:
 
   return Array.cast(__arg)
 
+@implements(numpy.shape)
+def shape(__arg: IntoArray) -> Tuple[int, ...]:
+    return asarray(__arg).shape
+
+@implements(numpy.ndim)
+def ndim(__arg: IntoArray) -> int:
+    return asarray(__arg).ndim
+
+@implements(numpy.size)
+def size(__arg: IntoArray) -> int:
+    return asarray(__arg).size
+
 def zeros(shape: Shape, dtype: DType = float) -> Array:
   '''Create a new :class:`Array` of given shape and dtype, filled with zeros.
 

From a970786d017873659da1d007c9c7828b26e53913 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Sat, 30 Oct 2021 14:33:57 +0200
Subject: [PATCH 02/11] lift function.Array restriction from Namespace

This patch removes the Namespace restriction in expression_v2 that
namespace elements must be instances of function.Array, instead relying
where possible on Numpy methods for array manipulation.
---
 nutils/expression_v2.py | 103 +++++++++++++++++++++-------------------
 1 file changed, 53 insertions(+), 50 deletions(-)

diff --git a/nutils/expression_v2.py b/nutils/expression_v2.py
index b901b820e..7fdddb8c8 100644
--- a/nutils/expression_v2.py
+++ b/nutils/expression_v2.py
@@ -595,33 +595,33 @@ class Namespace:
 
   def __init__(self) -> None:
     self.opposite = function.opposite
-    self.sin = function.sin
-    self.cos = function.cos
-    self.tan = function.tan
-    self.sinh = function.sinh
-    self.cosh = function.cosh
-    self.tanh = function.tanh
-    self.arcsin = function.arcsin
-    self.arccos = function.arccos
-    self.arctan = function.arctan
-    self.arctanh = function.arctanh
-    self.exp = function.exp
-    self.abs = function.abs
-    self.ln = function.ln
-    self.log = function.ln
-    self.log2 = function.log2
-    self.log10 = function.log10
-    self.sqrt = function.sqrt
-    self.sign = function.sign
-    self.conj = function.conj
-    self.real = function.real
-    self.imag = function.imag
+    self.sin = numpy.sin
+    self.cos = numpy.cos
+    self.tan = numpy.tan
+    self.sinh = numpy.sinh
+    self.cosh = numpy.cosh
+    self.tanh = numpy.tanh
+    self.arcsin = numpy.arcsin
+    self.arccos = numpy.arccos
+    self.arctan = numpy.arctan
+    self.arctanh = numpy.arctanh
+    self.exp = numpy.exp
+    self.abs = numpy.abs
+    self.ln = numpy.log
+    self.log = numpy.log
+    self.log2 = numpy.log2
+    self.log10 = numpy.log10
+    self.sqrt = numpy.sqrt
+    self.sign = numpy.sign
+    self.conj = numpy.conj
+    self.real = numpy.real
+    self.imag = numpy.imag
 
   def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None:
     name, underscore, indices = attr.partition('_')
     if isinstance(value, (int, float, complex, numpy.ndarray)):
       value = function.Array.cast(value)
-    if isinstance(value, function.Array):
+    if hasattr(value, '__array_ufunc__') and hasattr(value, '__array_function__'):
       if underscore:
         raise AttributeError('Cannot assign an array to an attribute with an underscore.')
       super().__setattr__(name, value)
@@ -632,7 +632,7 @@ def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None:
         raise AttributeError('All indices must be unique.')
       ops = _FunctionArrayOps(self)
       array, shape, expression_indices, summed = _Parser(ops).parse_expression(_Substring(value))
-      assert array.shape == shape
+      assert numpy.shape(array) == shape
       if expression_indices != indices:
         for index in sorted(set(indices) - set(expression_indices)):
           raise AttributeError('Index {} of the namespace attribute is missing in the expression.'.format(index))
@@ -640,7 +640,7 @@ def __setattr__(self, attr: str, value: Union[function.Array, str]) -> None:
           raise AttributeError('Index {} of the expression is missing in the namespace attribute.'.format(index))
         array = ops.align(array, expression_indices, indices)
       super().__setattr__(name, array)
-    elif isinstance(value, Callable):
+    elif callable(value):
       if underscore:
         raise AttributeError('Cannot assign a function to an attribute with an underscore.')
       super().__setattr__(name, value)
@@ -652,7 +652,7 @@ def __rmatmul__(self, expression):
     parser = _Parser(ops)
     if isinstance(expression, str):
       array, shape, indices, summed = parser.parse_expression(_Substring(expression))
-      assert array.shape == shape
+      assert numpy.shape(array) == shape
       array = ops.align(array, indices, ''.join(sorted(indices)))
       return array
     elif isinstance(expression, tuple):
@@ -702,17 +702,17 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio
     if gradient:
       setattr(self, gradient, lambda arg: function.grad(arg, geom))
     if curl:
-      if geom.shape != (3,):
-        raise ValueError('The curl can only be defined for a geometry with shape (3,) but got {}.'.format(geom.shape))
+      if numpy.shape(geom) != (3,):
+        raise ValueError('The curl can only be defined for a geometry with shape (3,) but got {}.'.format(numpy.shape(geom)))
       # Definition: `curl_ki(u_...)` := `ε_kji ∇_j(u_...)`. Should be used as
       # `curl_ki(u_i)`, which is equivalent to `ε_kji ∇_j(u_i)`.
       setattr(self, curl, lambda arg: (function.levicivita(3) * function.grad(arg, geom)[...,numpy.newaxis,:,numpy.newaxis]).sum(-2))
     if normal:
       setattr(self, normal, function.normal(geom))
     for i, jacobian in enumerate(jacobians):
-      if i > geom.size:
+      if i > numpy.size(geom):
         raise ValueError('Cannot define the jacobian {!r}: dimension is negative.'.format(jacobian))
-      setattr(self, jacobian, function.jacobian(geom, geom.size - i))
+      setattr(self, jacobian, function.jacobian(geom, numpy.size(geom) - i))
 
   def copy_(self, **replacements: Mapping[str, function.Array]) -> 'Namespace':
     '''Return a copy of this namespace.
@@ -755,12 +755,12 @@ def get_variable(self, name: str, ndim: int) -> Optional[Union[Tuple[function.Ar
       array = getattr(self.namespace, name)
     except AttributeError:
       return None
-    if not isinstance(array, function.Array):
+    if callable(array):
       return None
-    if array.ndim == ndim:
-      return array, array.shape
+    elif numpy.ndim(array) == ndim:
+      return array, numpy.shape(array)
     else:
-      return _InvalidDimension(array.ndim)
+      return _InvalidDimension(numpy.ndim(array))
 
   def call(self, name: str, ngenerates: int, arg: function.Array) -> Optional[Union[Tuple[function.Array, _Shape],_InvalidDimension]]:
     try:
@@ -768,23 +768,22 @@ def call(self, name: str, ngenerates: int, arg: function.Array) -> Optional[Unio
     except AttributeError:
       return None
     array = func(arg)
-    assert isinstance(array, function.Array)
-    assert array.shape[:arg.ndim] == arg.shape
-    if array.ndim == arg.ndim + ngenerates:
-      return array, array.shape[arg.ndim:]
+    assert numpy.shape(array)[:numpy.ndim(arg)] == numpy.shape(arg)
+    if numpy.ndim(array) == numpy.ndim(arg) + ngenerates:
+      return array, numpy.shape(array)[numpy.ndim(arg):]
     else:
-      return _InvalidDimension(array.ndim - arg.ndim)
+      return _InvalidDimension(numpy.ndim(array) - numpy.ndim(arg))
 
   def get_element(self, array: function.Array, axis: int, index: int) -> function.Array:
-    assert 0 <= axis < array.ndim and 0 <= index < array.shape[axis]
-    return function.get(array, axis, index)
+    assert 0 <= axis < numpy.ndim(array) and 0 <= index < numpy.shape(array)[axis]
+    return numpy.take(array, index, axis)
 
   def transpose(self, array: function.Array, axes: Tuple[int, ...]) -> function.Array:
-    assert array.ndim == len(axes)
-    return function.transpose(array, axes)
+    assert numpy.ndim(array) == len(axes)
+    return numpy.transpose(array, axes)
 
   def trace(self, array: function.Array, axis1: int, axis2: int) -> function.Array:
-    return function.trace(array, axis1, axis2)
+    return numpy.trace(array, axis1, axis2)
 
   def scope(self, array: function.Array) -> function.Array:
     return array
@@ -796,20 +795,24 @@ def jump(self, array: function.Array) -> function.Array:
     return function.jump(array)
 
   def add(self, *args: Tuple[bool, function.Array]) -> function.Array:
-    assert all(arg.shape == args[0][1].shape for neg, arg in args[1:])
+    assert all(numpy.shape(arg) == numpy.shape(args[0][1]) for neg, arg in args[1:])
     negated = (-arg if neg else arg for neg, arg in args)
-    return functools.reduce(function.add, negated)
+    return functools.reduce(numpy.add, negated)
+
+  def append_axes(self, array, shape):
+    shuffle = numpy.concatenate([len(shape) + numpy.arange(numpy.ndim(array)), numpy.arange(len(shape))])
+    return numpy.transpose(numpy.broadcast_to(array, shape + numpy.shape(array)), shuffle)
 
   def multiply(self, *args: function.Array) -> function.Array:
     result = args[0]
     for arg in args[1:]:
-      result = function.multiply(function._append_axes(result, arg.shape), function._prepend_axes(arg, result.shape))
+      result = numpy.multiply(self.append_axes(result, numpy.shape(arg)), arg)
     return result
 
   def divide(self, numerator: function.Array, denominator: function.Array) -> function.Array:
-    assert denominator.ndim == 0
-    return function.divide(numerator, function._append_axes(denominator, numerator.shape))
+    assert numpy.ndim(denominator) == 0
+    return numpy.true_divide(numerator, denominator)
 
   def power(self, base: function.Array, exponent: function.Array) -> function.Array:
-    assert exponent.ndim == 0
-    return function.power(base, function._append_axes(exponent, base.shape))
+    assert numpy.ndim(exponent) == 0
+    return numpy.power(base, exponent)

From 1a35c1ee5106a74face3f9bcc09ae68f2d9ccebf Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Sat, 18 Dec 2021 23:07:08 +0100
Subject: [PATCH 03/11] remove some array casts from __array_ufunc__

This patch removes some array casts from the ufunc handler so that certain
special values (most notably the number 1.) can be recognized for
low-hanging-fruit-type simplifications.
---
 nutils/function.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nutils/function.py b/nutils/function.py
index fa72e2bab..431a01e27 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -141,7 +141,7 @@ def __array_ufunc__(self, ufunc, method, *inputs, out=None, **kwargs):
     if method != '__call__' or ufunc not in HANDLED_FUNCTIONS:
       return NotImplemented
     try:
-      arrays = [Array.cast(v) for v in inputs]
+      arrays = [v if isinstance(v, (Array, bool, int, float, complex, numpy.ndarray)) else Array.cast(v) for v in inputs]
     except ValueError:
       return NotImplemented
     return HANDLED_FUNCTIONS[ufunc](*arrays, **kwargs)

From bda636d1521c925ba99edcf6071de2f5b0308771 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Sat, 18 Dec 2021 16:27:09 +0100
Subject: [PATCH 04/11] add simplifications for function.multiply, divide

This patch adds simplifications for function.multiply and divide in case one of
the operands is a unit scalar.
---
 nutils/function.py     | 14 ++++++++++++--
 tests/test_function.py | 18 ++++++++++++++++++
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/nutils/function.py b/nutils/function.py
index 431a01e27..261e4b876 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -1229,7 +1229,10 @@ def multiply(__left: IntoArray, __right: IntoArray) -> Array:
   :class:`Array`
   '''
 
-  return _Wrapper.broadcasted_arrays(evaluable.multiply, __left, __right)
+  left, right = typecast_arrays(__left, __right)
+  return right if _is_unit_scalar(__left) \
+    else left if _is_unit_scalar(__right) \
+    else _Wrapper.broadcasted_arrays(evaluable.multiply, left, right)
 
 @implements(numpy.true_divide)
 def divide(__dividend: IntoArray, __divisor: IntoArray) -> Array:
@@ -1244,7 +1247,10 @@ def divide(__dividend: IntoArray, __divisor: IntoArray) -> Array:
   :class:`Array`
   '''
 
-  return multiply(__dividend, reciprocal(__divisor))
+  dividend, divisor = typecast_arrays(__dividend, __divisor, min_dtype=float)
+  return dividend if _is_unit_scalar(__divisor) \
+    else reciprocal(divisor) if _is_unit_scalar(__dividend) \
+    else _Wrapper.broadcasted_arrays(evaluable.divide, dividend, divisor)
 
 @implements(numpy.floor_divide)
 def floor_divide(__dividend: IntoArray, __divisor: IntoArray) -> Array:
@@ -3648,3 +3654,7 @@ def f_dofs_coeffs(self, index: evaluable.Array) -> Tuple[evaluable.Array,evaluab
 def Namespace(*args, **kwargs):
   from .expression_v1 import Namespace
   return Namespace(*args, **kwargs)
+
+def _is_unit_scalar(v):
+  T = type(v)
+  return T in _dtypes and v == T(1)
diff --git a/tests/test_function.py b/tests/test_function.py
index 94ef4e309..31c05b577 100644
--- a/tests/test_function.py
+++ b/tests/test_function.py
@@ -1310,3 +1310,21 @@ def test_multiple(self):
     g = function.dotarg('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))
+
+class simplifications(TestCase):
+
+  def test_multiply(self):
+    f = function.Argument('test', shape=(2,3), dtype=float)
+    self.assertIs(f * 1, f)
+    self.assertIs(f * 1., f)
+    f = function.Argument('test', shape=(2,3), dtype=int)
+    self.assertIs(f * 1, f)
+    self.assertIsNot(f * 1., f)
+
+  def test_divide(self):
+    f = function.Argument('test', shape=(2,3), dtype=float)
+    self.assertIs(f / 1, f)
+    self.assertIs(f / 1., f)
+    f = function.Argument('test', shape=(2,3), dtype=int)
+    self.assertIsNot(f / 1, f)
+    self.assertIsNot(f / 1., f)

From 4a2c15b344191094bbecb0cc9b94ef882b54206f Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 22 Oct 2021 11:18:47 +0200
Subject: [PATCH 05/11] add Array.cast_withscale

This patch adds support for dimensional values via function.Arrays' new
cast_withscale method. While primarily aimed at supporting nutils-SI
type quanties, the functionality supports any wrapper type that has an
attribute `reference_quantity` that can be used to make a value dimensionless
through division.
---
 nutils/function.py | 48 ++++++++++++++++++++++++++++++++--------------
 nutils/sample.py   | 13 ++++++++-----
 2 files changed, 42 insertions(+), 19 deletions(-)

diff --git a/nutils/function.py b/nutils/function.py
index 261e4b876..6dfee58a6 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -177,6 +177,17 @@ def cast(cls, __value: IntoArray, dtype: Optional[DType] = None, ndim: Optional[
       raise ValueError('expected an array with dimension `{}` but got `{}`'.format(ndim, value.ndim))
     return value
 
+  @classmethod
+  def cast_withscale(cls, __value: IntoArray, dtype: Optional[DType] = None, ndim: Optional[int] = None):
+    try:
+      scale = type(__value).reference_quantity
+    except AttributeError:
+      value = cls.cast(__value, dtype=dtype, ndim=ndim)
+      scale = value.dtype(1)
+    else:
+      value = cls.cast(__value / scale, dtype=dtype, ndim=ndim)
+    return value, scale
+
   def __init__(self, shape: Shape, dtype: DType, spaces: FrozenSet[str], arguments: Mapping[str, Tuple[Shape, DType]]) -> None:
     self.shape = tuple(sh.__index__() for sh in shape)
     self.dtype = dtype
@@ -2732,7 +2743,7 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array:
   :class:`Array`
   '''
 
-  arg = Array.cast(__arg)
+  arg, argscale = Array.cast_withscale(__arg)
   if isinstance(__var, str):
     if __var not in arg.arguments:
       raise ValueError('no such argument: {}'.format(__var))
@@ -2746,7 +2757,7 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array:
       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:
       raise ValueError('Argument {!r} has dtype {} in the function, but the derivative to {!r} with dtype {} was requested.'.format(__var.name, dtype.__name__ if dtype in _dtypes else dtype, __var.name, __var.dtype.__name__ if __var.dtype in _dtypes else __var.dtype))
-  return _Derivative(arg, __var)
+  return _Derivative(arg, __var) * argscale
 
 def grad(__arg: IntoArray, __geom: IntoArray, ndims: int = 0) -> Array:
   '''Return the gradient of the argument to the given geometry.
@@ -2762,21 +2773,22 @@ def grad(__arg: IntoArray, __geom: IntoArray, ndims: int = 0) -> Array:
   :class:`Array`
   '''
 
-  arg = Array.cast(__arg)
-  geom = Array.cast(__geom)
+  arg, argscale = Array.cast_withscale(__arg)
+  geom, geomscale = Array.cast_withscale(__geom)
   if geom.dtype != float:
     raise ValueError('The geometry must be real-valued.')
   if geom.ndim == 0:
-    return grad(arg, _append_axes(geom, (1,)))[...,0]
+    ret = grad(arg, _append_axes(geom, (1,)))[...,0]
   elif geom.ndim > 1:
     sh = geom.shape[-2:]
-    return unravel(grad(arg, ravel(geom, geom.ndim-2)), arg.ndim+geom.ndim-2, sh)
+    ret = unravel(grad(arg, ravel(geom, geom.ndim-2)), arg.ndim+geom.ndim-2, sh)
   elif ndims == 0 or ndims == geom.shape[0]:
-    return _Gradient(arg, geom)
+    ret = _Gradient(arg, geom)
   elif ndims == -1 or ndims == geom.shape[0] - 1:
-    return _SurfaceGradient(arg, geom)
+    ret = _SurfaceGradient(arg, geom)
   else:
     raise NotImplementedError
+  return ret * (argscale / geomscale)
 
 def curl(__arg: IntoArray, __geom: IntoArray) -> Array:
   '''Return the curl of the argument w.r.t. the given geometry.
@@ -2819,7 +2831,7 @@ def normal(__geom: IntoArray, refgeom: Optional[Array] = None) -> Array:
   :class:`Array`
   '''
 
-  geom = Array.cast(__geom)
+  geom, geomscale_ = Array.cast_withscale(__geom)
   if geom.dtype != float:
     raise ValueError('The geometry must be real-valued.')
   if geom.ndim == 0:
@@ -2888,15 +2900,22 @@ def jacobian(__geom: IntoArray, __ndims: Optional[int] = None) -> Array:
   :class:`Array`
   '''
 
-  geom = Array.cast(__geom)
+  geom, geomscale = Array.cast_withscale(__geom)
   if geom.dtype != float:
     raise ValueError('The geometry must be real-valued.')
+  if __ndims is not None:
+    scale = geomscale**__ndims
+  elif _is_unit_scalar(geomscale):
+    scale = 1.
+  else:
+    raise ValueError('scaled arrays require an explicit dimension')
   if geom.ndim == 0:
-    return jacobian(insertaxis(geom, 0, 1), __ndims)
+    J = jacobian(insertaxis(geom, 0, 1), __ndims)
   elif geom.ndim > 1:
-    return jacobian(ravel(geom, geom.ndim-2), __ndims)
+    J = jacobian(ravel(geom, geom.ndim-2), __ndims)
   else:
-    return _Jacobian(geom, __ndims)
+    J = _Jacobian(geom, __ndims)
+  return J * scale
 
 def J(__geom: IntoArray, __ndims: Optional[int] = None) -> Array:
   '''Return the absolute value of the determinant of the Jacobian matrix of the given geometry.
@@ -3044,7 +3063,8 @@ def eval(funcs: evaluable.AsEvaluableArray, **arguments: Mapping[str, numpy.ndar
   results : :class:`tuple` of arrays
   '''
 
-  return map(sparse.toarray, evaluable.eval_sparse(funcs, **arguments))
+  funcs, funcscales = zip(*map(Array.cast_withscale, funcs))
+  return tuple(sparse.toarray(data) * scale for data, scale in zip(evaluable.eval_sparse(funcs, **arguments), funcscales))
 
 def isarray(__arg: Any) -> bool:
   'Test if the argument is an instance of :class:`Array`.'
diff --git a/nutils/sample.py b/nutils/sample.py
index 469762cd7..0e155a6b5 100644
--- a/nutils/sample.py
+++ b/nutils/sample.py
@@ -189,9 +189,10 @@ def integrate(self, funcs: Iterable[function.IntoArray], arguments: Mapping[str,
         Optional arguments for function evaluation.
     '''
 
+    funcs, funcscales = zip(*map(function.Array.cast_withscale, funcs))
     datas = self.integrate_sparse(funcs, argdict(arguments))
     with log.iter.fraction('assembling', datas) as items:
-      return tuple(_convert(data, inplace=True) for data in items)
+      return tuple(_convert(data, inplace=True) * scale for data, scale in zip(items, funcscales))
 
   @util.single_or_multiple
   def integrate_sparse(self, funcs: Iterable[function.IntoArray], arguments: Optional[Mapping[str, numpy.ndarray]] = None) -> Tuple[numpy.ndarray, ...]:
@@ -216,7 +217,8 @@ def integral(self, __func: function.IntoArray) -> function.Array:
         Integrand.
     '''
 
-    return _Integral(function.Array.cast(__func), self)
+    func, funcscale = function.Array.cast_withscale(__func)
+    return _Integral(func, self) * funcscale
 
   @util.positional_only
   @util.single_or_multiple
@@ -231,9 +233,10 @@ def eval(self, funcs: Iterable[function.IntoArray], arguments: Mapping[str, nump
         Optional arguments for function evaluation.
     '''
 
+    funcs, funcscales = zip(*map(function.Array.cast_withscale, funcs))
     datas = self.eval_sparse(funcs, arguments)
     with log.iter.fraction('assembling', datas) as items:
-      return tuple(map(sparse.toarray, items))
+      return tuple(sparse.toarray(data) * funcscale for data, funcscale in zip(datas, funcscales))
 
   @util.positional_only
   @util.single_or_multiple
@@ -251,10 +254,10 @@ def eval_sparse(self, funcs: Iterable[function.IntoArray], arguments: Optional[M
     return evaluable.eval_sparse(map(self, funcs), **(arguments or {}))
 
   def __call__(self, __func: function.IntoArray) -> function.Array:
-    func = _ConcatenatePoints(function.Array.cast(__func), self)
+    func, funcscale = function.Array.cast_withscale(__func)
     ielem = evaluable.loop_index('_sample_' + '_'.join(self.spaces), self.nelems)
     indices = evaluable.loop_concatenate(evaluable._flat(self.get_evaluable_indices(ielem)), ielem)
-    return _ReorderPoints(func, indices)
+    return _ReorderPoints(_ConcatenatePoints(func, self), indices) * funcscale
 
   def basis(self) -> function.Array:
     '''Basis-like function that for every point in the sample evaluates to the

From 8b2ad09b1f1b5353f88e085dc2e3cf7c96bad435 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 2 Dec 2021 16:35:14 +0100
Subject: [PATCH 06/11] restructure export.triplot

This patch cleans up export.triplot without changing any functionality.
---
 nutils/export.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/nutils/export.py b/nutils/export.py
index 474dbe65f..f7bce9f29 100644
--- a/nutils/export.py
+++ b/nutils/export.py
@@ -54,15 +54,20 @@ def mplfigure(name, kwargs=...):
       finally:
         fig.set_canvas(None) # break circular reference
 
+def plotlines_(ax, xy, lines, **kwargs):
+  from matplotlib import collections
+  lc = collections.LineCollection(numpy.asarray(xy).T[lines], **kwargs)
+  ax.add_collection(lc)
+  return lc
+
 def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k'):
   if (tri is None) != (values is None):
     raise Exception('tri and values can only be specified jointly')
   with mplfigure(name) as fig:
-    ax = fig.add_subplot(111)
     if points.shape[1] == 1:
+      ax = fig.add_subplot(111)
       if tri is not None:
-        import matplotlib.collections
-        ax.add_collection(matplotlib.collections.LineCollection(numpy.array([points[:,0], values]).T[tri]))
+        plotlines_(ax, [points[:,0], values], tri)
       if hull is not None:
         for x in points[hull[:,0],0]:
           ax.axvline(x, color=linecolor, linewidth=linewidth)
@@ -72,15 +77,14 @@ def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=N
       else:
         ax.set_ylim(clim)
     elif points.shape[1] == 2:
-      ax.set_aspect('equal')
+      ax = fig.add_subplot(111, aspect='equal')
       if tri is not None:
-        im = ax.tripcolor(points[:,0], points[:,1], tri, values, shading='gouraud', cmap=cmap, rasterized=True)
+        im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True)
         if clim is not None:
           im.set_clim(clim)
         fig.colorbar(im)
       if hull is not None:
-        import matplotlib.collections
-        ax.add_collection(matplotlib.collections.LineCollection(points[hull], colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5))
+        plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5)
       ax.autoscale(enable=True, axis='both', tight=True)
     else:
       raise Exception('invalid spatial dimension: {}'.format(points.shape[1]))

From 63f5613afcbe9dcf585ba533e271c3433a95e3ff Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Sat, 30 Oct 2021 14:24:54 +0200
Subject: [PATCH 07/11] add plabel, vlabel arguments to export.triplot

This patch adds the plabel and vlabel arguments to export.triplot for
the coordinates and values, respectively. For 1D plots the labels are
used for the x and y axes; for 2D plots the plabel is used for both x
and y and the vlabel is used for the colorbar.
---
 nutils/export.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/nutils/export.py b/nutils/export.py
index f7bce9f29..a7894f0b9 100644
--- a/nutils/export.py
+++ b/nutils/export.py
@@ -60,12 +60,12 @@ def plotlines_(ax, xy, lines, **kwargs):
   ax.add_collection(lc)
   return lc
 
-def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k'):
+def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None):
   if (tri is None) != (values is None):
     raise Exception('tri and values can only be specified jointly')
   with mplfigure(name) as fig:
     if points.shape[1] == 1:
-      ax = fig.add_subplot(111)
+      ax = fig.add_subplot(111, xlabel=plabel, ylabel=vlabel)
       if tri is not None:
         plotlines_(ax, [points[:,0], values], tri)
       if hull is not None:
@@ -77,12 +77,12 @@ def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=N
       else:
         ax.set_ylim(clim)
     elif points.shape[1] == 2:
-      ax = fig.add_subplot(111, aspect='equal')
+      ax = fig.add_subplot(111, xlabel=plabel, ylabel=plabel, aspect='equal')
       if tri is not None:
         im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True)
         if clim is not None:
           im.set_clim(clim)
-        fig.colorbar(im)
+        fig.colorbar(im, label=vlabel)
       if hull is not None:
         plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5)
       ax.autoscale(enable=True, axis='both', tight=True)

From dc1d784aa4c6ae0c218b87d803a30426cb41ce99 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 2 Dec 2021 15:59:47 +0100
Subject: [PATCH 08/11] add numeric.simplex_grid

See the docstring of numeric.simplex_grid for more info.
---
 nutils/numeric.py     | 48 +++++++++++++++++++++++++++++++++++++++++++
 tests/test_numeric.py | 25 ++++++++++++++++++++++
 2 files changed, 73 insertions(+)

diff --git a/nutils/numeric.py b/nutils/numeric.py
index da55e0698..a9b72f0d0 100644
--- a/nutils/numeric.py
+++ b/nutils/numeric.py
@@ -145,6 +145,54 @@ def meshgrid(*args, dtype=None):
   assert n == 0
   return grid
 
+def _simplex_grid_helper(bounds):
+  if bounds.ndim != 1 or len(bounds) == 0:
+    raise ValueError
+  nd = len(bounds)
+  spacing = [numpy.sqrt((1+i/2)/(1+i)) for i in range(nd)] # simplex height orthogonal to lower dimension
+  grid = meshgrid(*[numpy.arange(bound, step=step) for step, bound in zip(spacing, bounds)])
+  out_of_bounds = []
+  for idim in range(nd-1):
+    stripes = grid[(idim,)+(slice(None),)*(idim+1)+(slice(1,None,2),)]
+    stripes += spacing[idim] * (idim+1) / (idim+2)
+    if stripes.size and stripes.flat[-1] >= bounds[idim]:
+      out_of_bounds.append(idim)
+  if out_of_bounds:
+    select = numpy.ones(grid.shape[1:], dtype=bool)
+    for idim in out_of_bounds:
+      select[(slice(None),)*(idim)+(-1,)+(slice(1,None,2),)] = False
+    points = grid[:,select].T
+  else:
+    points = grid.reshape(nd, -1).T
+  d = numpy.subtract(bounds, points.max(axis=0))
+  assert (d > 0).all()
+  points += d / 2
+  return points
+
+def simplex_grid(shape, spacing):
+  '''Multi-dimensional generator for equilateral simplex grids.
+
+  Simplex_grid generates a point cloud within an n-dimensional orthotope, which
+  ranges from zero to a specified shape. The point coordinates are spaced in
+  such a way that the nearest neighbours are at distance `spacing`, thus
+  forming vertices of regular simplices. The returned array is two-dimensional,
+  with the first axis being the spatial dimension (matching `shape`) and the
+  second a stacking of the generated points.
+
+  Parameters
+  ----------
+  shape : :class:`tuple`
+    list or tuple of dimensions of the orthotope to be filled.
+  spacing : :class:`float`
+    minimum spacing in the generated point cloud.
+
+  Returns
+  -------
+  :class:`numpy.ndarray`
+  '''
+
+  return _simplex_grid_helper(numpy.divide(shape, spacing)) * spacing
+
 def takediag(A, axis=-2, rmaxis=-1):
   axis = normdim(A.ndim, axis)
   rmaxis = normdim(A.ndim, rmaxis)
diff --git a/tests/test_numeric.py b/tests/test_numeric.py
index d18999c06..b24f725a7 100644
--- a/tests/test_numeric.py
+++ b/tests/test_numeric.py
@@ -226,6 +226,31 @@ def test_dtype(self):
     self.assertEqual(m.dtype, float)
     self.assertAllEqual(m, numpy.ones((1,)))
 
+class simplex_grid(TestCase):
+
+  def simplex_grid(self, shape, spacing):
+    coords = numeric.simplex_grid(shape, spacing)
+    self.assertEqual(coords.ndim, 2)
+    self.assertEqual(coords.shape[1], len(shape))
+    self.assertTrue((coords > 0).all())
+    self.assertTrue((coords < shape).all())
+    mindist = min(numpy.linalg.norm(c1 - c2) for i, c1 in enumerate(coords) for c2 in coords[:i])
+    self.assertAlmostEqual(mindist, spacing)
+    return coords
+
+  def test_1d(self):
+    coords = self.simplex_grid([2], .8)
+    self.assertEqual(len(coords), 3)
+    self.assertAllAlmostEqual(coords[:,0], [.2,1,1.8])
+
+  def test_2d(self):
+    coords = self.simplex_grid([2,3], .8)
+    self.assertEqual(len(coords), 13)
+
+  def test_3d(self):
+    coords = self.simplex_grid([2,3,4], .8)
+    self.assertEqual(len(coords), 82)
+
 class overlapping(TestCase):
 
   def test_pairwise(self):

From a076dfa28e14a11a5d70739b5a8f35c558300778 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 2 Dec 2021 16:24:04 +0100
Subject: [PATCH 09/11] add skip_missing argument to Topology.locate

This patch adds the skip_missing argument to Topology.locate, allowing
it to silently drop points that could not be located rather than raising
a LocateError.
---
 nutils/topology.py | 47 ++++++++++++++++++++++++++++++++--------------
 1 file changed, 33 insertions(+), 14 deletions(-)

diff --git a/nutils/topology.py b/nutils/topology.py
index 66b3fac5a..5463500a8 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -697,7 +697,7 @@ def select(self, indicator: function.Array, ischeme: str = 'bezier2', **kwargs:
     selected = types.frozenarray(numpy.unique(ielem[isactive]))
     return self[selected]
 
-  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample:
+  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample:
     '''Create a sample based on physical coordinates.
 
     In a finite element application, functions are commonly evaluated in points
@@ -745,6 +745,10 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
         element centroid above which the element is rejected immediately. If
         all points are expected to be located then this can safely be left
         unspecified.
+    skip_missing : :class:`bool` (default: False)
+        When set to true, skip points that are not found (for instance because
+        they fall outside the domain) in the returned sample. When set to false
+        (the default) missing points raise a ``LocateError``.
 
     Returns
     -------
@@ -950,7 +954,7 @@ def withgroups(self, vgroups: Mapping[str, Union[str, Topology]] = {}, bgroups:
     def indicator(self, subtopo: Union[str, Topology]) -> Topology:
       raise SkipTest('`{}` does not implement `Topology.indicator`'.format(type(self).__qualname__))
 
-    def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample:
+    def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample:
       raise SkipTest('`{}` does not implement `Topology.locate`'.format(type(self).__qualname__))
 
 else:
@@ -1269,8 +1273,8 @@ def indicator(self, subtopo: Union[str, Topology]) -> Topology:
     else:
       return super().indicator(subtopo)
 
-  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None) -> Sample:
-    return self.parent.locate(geom, coords, tol=tol, eps=eps, maxiter=maxiter, arguments=arguments, weights=weights, maxdist=maxdist, ischeme=ischeme, scale=scale)
+  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample:
+    return self.parent.locate(geom, coords, tol=tol, eps=eps, maxiter=maxiter, arguments=arguments, weights=weights, maxdist=maxdist, ischeme=ischeme, scale=scale, skip_missing=skip_missing)
 
   def boundary_spaces_unchecked(self, spaces: FrozenSet[str]) -> Topology:
     return _WithGroupAliases(self.parent.boundary_spaces_unchecked(spaces), self.bgroups, types.frozendict({}), types.frozendict({}))
@@ -1473,7 +1477,7 @@ def indicator(self, subtopo):
     return function.get(values, 0, self.f_index)
 
   @log.withcontext
-  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None):
+  def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False):
     if ischeme is not None:
       warnings.deprecation('the ischeme argument is deprecated and will be removed in future')
     if scale is not None:
@@ -1495,6 +1499,10 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
     egeom = geom.lower((), {self.space: (self.transforms.get_evaluable(_ielem), self.opposites.get_evaluable(_ielem))}, {self.space: _point})
     xJ = evaluable.Tuple((egeom, evaluable.derivative(egeom, _point))).simplified
     arguments = dict(arguments or ())
+    if skip_missing:
+      if weights is not None:
+        raise ValueError('weights and skip_missing are mutually exclusive')
+      missing = parallel.shzeros(len(coords), dtype=bool)
     with parallel.ctxrange('locating', len(coords)) as ipoints:
       for ipoint in ipoints:
         xt = coords[ipoint] # target
@@ -1528,7 +1536,13 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
               points[ipoint] = p
               break
         else:
-          raise LocateError('failed to locate point: {}'.format(xt))
+          if skip_missing:
+            missing[ipoint] = True
+          else:
+            raise LocateError('failed to locate point: {}'.format(xt))
+    if skip_missing:
+      ielems = ielems[~missing]
+      points = points[~missing]
     return self._sample(ielems, points, weights)
 
   def _sample(self, ielems, coords, weights=None):
@@ -2215,7 +2229,7 @@ def refined(self):
     axes = [axis.refined for axis in self.axes]
     return StructuredTopology(self.space, self.root, axes, self.nrefine+1, bnames=self._bnames)
 
-  def locate(self, geom, coords, *, tol, eps=0, weights=None, **kwargs):
+  def locate(self, geom, coords, *, tol, eps=0, weights=None, skip_missing=False, **kwargs):
     coords = numpy.asarray(coords, dtype=float)
     if geom.ndim == 0:
       geom = geom[_]
@@ -2225,9 +2239,9 @@ def locate(self, geom, coords, *, tol, eps=0, weights=None, **kwargs):
     geom0, scale, index = self._asaffine(geom)
     e = self.sample('uniform', 2).eval(function.norm2(geom0 + index * scale - geom)).max() # inf-norm on non-gauss sample
     if e > tol:
-      return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, **kwargs)
+      return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, skip_missing=skip_missing, **kwargs)
     log.info('locate detected linear geometry: x = {} + {} xi ~{:+.1e}'.format(geom0, scale, e))
-    return self._locate(geom0, scale, coords, eps=eps, weights=weights)
+    return self._locate(geom0, scale, coords, eps=eps, weights=weights, skip_missing=skip_missing)
 
   def _asaffine(self, geom):
     p0 = p1 = self
@@ -2243,14 +2257,19 @@ def _asaffine(self, geom):
     scale = (geom1 - geom0) / numpy.array(self.shape)
     return geom0, scale, index
 
-  def _locate(self, geom0, scale, coords, *, eps=0, weights=None):
+  def _locate(self, geom0, scale, coords, *, eps=0, weights=None, skip_missing=False):
     mincoords, maxcoords = numpy.sort([geom0, geom0 + scale * self.shape], axis=0)
-    outofbounds = numpy.less(coords, mincoords - eps) | numpy.greater(coords, maxcoords + eps)
-    if outofbounds.any():
-      raise LocateError('failed to locate {}/{} points'.format(outofbounds.any(axis=1).sum(), len(coords)))
+    missing = numpy.any(numpy.less(coords, mincoords - eps) | numpy.greater(coords, maxcoords + eps), axis=1)
+    if not skip_missing and missing.any():
+      raise LocateError('failed to locate {}/{} points'.format(missing.sum(), len(coords)))
     xi = (coords - geom0) / scale
     ielem = numpy.minimum(numpy.maximum(xi.astype(int), 0), numpy.array(self.shape)-1)
-    return self._sample(numpy.ravel_multi_index(ielem.T, self.shape), xi - ielem, weights)
+    ielems = numpy.ravel_multi_index(ielem.T, self.shape)
+    points = xi - ielem
+    if skip_missing:
+      ielems = ielems[~missing]
+      points = points[~missing]
+    return self._sample(ielems, points, weights)
 
   def __str__(self):
     'string representation'

From d0bb13e0e14395e35287e50c5ad55cec565d2c07 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Sat, 30 Oct 2021 14:23:54 +0200
Subject: [PATCH 10/11] make cahnhilliard example dimensional

This patch changes the cahnhilliard example to use dimensional values.
---
 .github/workflows/test.yaml |   2 +-
 examples/cahnhilliard.py    | 281 ++++++++++++++++++++++++------------
 2 files changed, 191 insertions(+), 92 deletions(-)

diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index a8175fa97..003423841 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -134,7 +134,7 @@ jobs:
         id: install
         run: |
           python -um pip install --upgrade wheel
-          python -um pip install --upgrade treelog stringly matplotlib scipy pillow numpy
+          python -um pip install --upgrade treelog stringly matplotlib scipy pillow numpy git+https://github.com/evalf/nutils-SI.git
           # Install Nutils from `dist` dir created in job `build-python-package`.
           python -um pip install --no-index --find-links ./dist nutils
       - name: Test
diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py
index e653cf5fa..58ad3ac4b 100644
--- a/examples/cahnhilliard.py
+++ b/examples/cahnhilliard.py
@@ -1,101 +1,203 @@
 #! /usr/bin/env python3
 #
-# In this script we solve the Cahn-Hiilliard equation, which models the
-# unmixing of two phases under the effect of surface tension.
+# In this script we solve the Cahn-Hilliard equation, which models the unmixing
+# of two phases (φ=+1 and φ=-1) under the influence of surface tension. It is a
+# mixed equation of two scalar equations for phase φ and chemical potential η::
+#
+#     dφ/dt = -div(J(η))
+#     η = ψ'(φ) σ / ε - σ ε Δ(φ)
+#
+# along with constitutive relations for the flux vector J = -M ∇η and the
+# double well potential ψ = .25 (φ² - 1)², and subject to boundary conditions
+# ∂ₙφ = -σd / σ ε and ∂ₙη = 0. Parameters are the interface thickness ε, fluid
+# surface tension σ, differential wall surface tension σd, and mobility M.
+#
+# Cahn-Hilliard is a diffuse interface model, which means that phases do not
+# separate sharply, but instead form a transition zone between the phases. The
+# transition zone has a thickness proportional to ε, as is readily confirmed in
+# one dimension, where a steady state solution on an infinite domain is formed
+# by η(x) = 0, φ(x) = tanh(x / √2 ε).
+#
+# The main driver of the unmixing process is the double well potential ψ that
+# is proportional to the mixing energy, taking its minima at the pure phases
+# φ=+1 and φ=-1. The interface between the phases adds an energy contribution
+# proportional to its length. At the wall we have a phase-dependent fluid-solid
+# energy. Over time, the system minimizes the total energy::
+#
+#     E(φ) = ∫_Ω ψ(φ) σ / ε + ∫_Ω .5 σ ε ‖∇φ‖² + ∫_Γ (σm + φ σd)
+#            \                \                  \
+#             mixing energy    interface energy   wall energy
+#
+# Proof: the time derivative of E followed by substitution of the strong form
+# and boundary conditions yields dE/dt = ∫_Ω η dφ/dt = -∫_Ω M ‖∇η‖² ≤ 0. □
+#
+# Switching to discrete time we set dφ/dt = (φ - φ0) / dt and add a stabilizing
+# perturbation term δψ(φ, φ0) to the doube well potential for reasons outlined
+# below. This yields the following time discrete system::
+#
+#     φ = φ0 - dt div(J(η))
+#     η = (ψ'(φ) + δψ'(φ, ψ0)) σ / ε - σ ε Δ(φ)
+#
+# with the equivalent weak formulation::
+#
+#     ∂/∂η ∫_Ω [ η (φ - φ0) + .5 dt J(η)·∇η ] = 0
+#     ∂/∂φ ∫_Ω [ E(φ) + δψ(φ, φ0) σ / ε - η φ ] = 0
+#
+# For stability we wish for the perturbation δψ to be such that the time
+# discrete system preserves the energy dissipation property E(φ) ≤ E(φ0) for
+# any timestep dt. To derive suitable perturbation terms to this effect, we
+# define without loss of generality δψ'(φ, φ0) = .5 (φ - φ0) f(φ, φ0) and
+# derive the following condition for unconditional stability::
+#
+#     E(φ) - E(φ0) = ∫_Ω .5 (1 - φ² - .5 (φ + φ0)² - f(φ, φ0)) (φ - φ0)² σ / ε
+#                  - ∫_Ω (.5 σ ε ‖∇φ - ∇φ0‖² + dt M ‖∇η‖²) ≤ 0
+#
+# The inequality holds true if the perturbation f is bounded from below such
+# that f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)². To keep the energy minima at the pure
+# phases we additionally impose that f(±1, φ0) = 0, and select 1 - φ² as a
+# suitable upper bound which we will call "nonlinear".
+#
+# We next observe that η is linear in φ if f(φ, φ0) = g(φ0) - φ² - (φ + φ0)²
+# for any function g, which dominates if g(φ0) ≥ 1 + .5 (φ + φ0)². While this
+# cannot be made to hold true for all φ, we make it hold for -√2 ≤ φ, φ0 ≤ √2
+# by defining g(φ0) = 2 + 2 ‖φ0‖ + φ0², which we will call "linear". This
+# scheme further satisfies a weak minima preservation f(±1, ±‖φ0‖) = 0.
+#
+# We have thus arrived at the three stabilization schemes implemented here:
+#
+# - nonlinear: f(φ, φ0) = 1 - φ²
+# - linear: f(φ, φ0) = 2 + 2 ‖φ0‖ - 2 φ (φ + φ0)
+# - none: f(φ, φ0) = 0 (not unconditionally stable)
+#
+# The stab enum in this script defines the schemes in terms of δψ to allow
+# Nutils to construct the residuals through automatic differentiation.
+#
+# NOTE: This script uses dimensional quantities and requires the nutils.SI
+# module, which is installed separate from the the core nutils.
 
-from nutils import mesh, function, solver, export, cli, testing
+from nutils import mesh, function, solver, numeric, export, cli, testing
 from nutils.expression_v2 import Namespace
-import numpy, treelog, itertools, enum, typing
+from nutils.SI import Length, Time, Density, Tension, Energy, Pressure, Velocity, parse
+import numpy, itertools, enum
+import treelog as log
 
 class stab(enum.Enum):
-  none = '0' # for educational purposes only
-  linear = '.5 dc^2 (6 - 6 c^2 + 8 c dc - 3 dc^2) / epsilon^2'
-  optimal = '.5 dc^2 (1 - dc^2 / 12) / epsilon^2'
-
-def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[float],
-         contactangle:float, timestep:float, mtol:float, seed:int, circle:bool, stab:stab):
+  nonlinear = '.25 dφ^2 (1 - φ^2 + φ dφ (2 / 3) - dφ^2 / 6)'
+  linear = '.25 dφ^2 (2 + 2 abs(φ0) - (φ + φ0)^2)'
+  none = '0' # not unconditionally stable
+
+def main(size:Length, epsilon:Length, mobility:Time/Density, stens:Tension,
+         wtensn:Tension, wtensp:Tension, nelems:int, etype:str, degree:int,
+         timestep:Time, tol:Energy/Length, endtime:Time, seed:int, circle:bool,
+         stab:stab, showflux:bool):
   '''
   Cahn-Hilliard equation on a unit square/circle.
 
   .. arguments::
 
-     nelems [20]
-       Number of elements along domain edge.
+     size [10cm]
+       Domain size.
+     epsilon [2mm]
+       Interface thickness; defaults to an automatic value based on the
+       configured mesh density if left unspecified.
+     mobility [1mL*s/kg]
+       Mobility.
+     stens [50mN/m]
+       Surface tension.
+     wtensn [30mN/m]
+       Wall surface tension for phase -1.
+     wtensp [20mN/m]
+       Wall surface tension for phase +1.
+     nelems [0]
+       Number of elements along domain edge. When set to zero a value is set
+       automatically based on the configured domain size and epsilon.
      etype [square]
        Type of elements (square/triangle/mixed).
-     btype [std]
-       Type of basis function (std/spline), with availability depending on the
-       configured element type.
      degree [2]
        Polynomial degree.
-     epsilon []
-       Interface thickness; defaults to an automatic value based on the
-       configured mesh density if left unspecified.
-     contactangle [90]
-       Wall contact angle in degrees.
-     timestep [.01]
+     timestep [.5s]
        Time step.
-     mtol [.01]
-       Threshold value for chemical potential peak to peak difference, used as
-       a stop criterion.
+     tol [1nJ/m]
+       Newton tolerance.
+     endtime [1min]
+       End of the simulation.
      seed [0]
        Random seed for the initial condition.
      circle [no]
        Select circular domain as opposed to a unit square.
      stab [linear]
-       Stabilization method (linear/optimal/none).
+       Stabilization method (linear/nonlinear/none).
+     showflux [no]
+       Overlay flux vectors on phase plot
   '''
 
-  mineps = 1./nelems
-  if epsilon is None:
-    treelog.info('setting epsilon={}'.format(mineps))
-    epsilon = mineps
-  elif epsilon < mineps:
-    treelog.warning('epsilon under crititical threshold: {} < {}'.format(epsilon, mineps))
+  nmin = round(size / epsilon)
+  if nelems <= 0:
+    nelems = nmin
+    log.info('setting nelems to {}'.format(nelems))
+  elif nelems < nmin:
+    log.warning('mesh is too coarse, consider increasing nelems to {:.0f}'.format(nmin))
+
+  log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi))
 
   domain, geom = mesh.unitsquare(nelems, etype)
-  bezier = domain.sample('bezier', 5) # sample for plotting
+  if circle:
+    angle = (geom-.5) * (numpy.pi/2)
+    geom = .5 + function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2)
+
+  bezier = domain.sample('bezier', 5) # sample for surface plots
+  grid = domain.locate(geom, numeric.simplex_grid([1,1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5) # sample for quivers
+
+  φbasis = ηbasis = domain.basis('std', degree=degree)
+  ηbasis *= stens / epsilon # basis scaling to give η the required unit
 
   ns = Namespace()
-  if not circle:
-    ns.x = geom
-  else:
-    angle = (geom-.5) * (numpy.pi/2)
-    ns.x = function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2)
+  ns.x = size * geom
   ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-  ns.epsilon = epsilon
-  ns.ewall = .5 * numpy.cos(contactangle * numpy.pi / 180)
-  ns.cbasis = ns.mbasis = domain.basis('std', degree=degree)
-  ns.c = function.dotarg('c', ns.cbasis)
-  ns.dc = ns.c - function.dotarg('c0', ns.cbasis)
-  ns.m = function.dotarg('m', ns.mbasis)
-  ns.F = '.5 (c^2 - 1)^2 / epsilon^2'
-  ns.dF = stab.value
+  ns.ε = epsilon
+  ns.σ = stens
+  ns.φ = function.dotarg('φ', φbasis)
+  ns.σmean = (wtensp + wtensn) / 2
+  ns.σdiff = (wtensp - wtensn) / 2
+  ns.σwall = 'σmean + φ σdiff'
+  ns.φ0 = function.dotarg('φ0', φbasis)
+  ns.dφ = 'φ - φ0'
+  ns.η = function.dotarg('η', ηbasis)
+  ns.ψ = '.25 (φ^2 - 1)^2'
+  ns.δψ = stab.value
+  ns.M = mobility
+  ns.J_i = '-M ∇_i(η)'
   ns.dt = timestep
 
-  nrg_mix = domain.integral('F dV' @ ns, degree=7)
-  nrg_iface = domain.integral('.5 ∇_k(c) ∇_k(c) dV' @ ns, degree=7)
-  nrg_wall = domain.boundary.integral('(abs(ewall) + c ewall) dS' @ ns, degree=7)
-  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 ∇_k(m) ∇_k(m)) dV' @ ns, degree=7)
+  nrg_mix = domain.integral('(ψ σ / ε) dV' @ ns, degree=7)
+  nrg_iface = domain.integral('.5 σ ε ∇_k(φ) ∇_k(φ) dV' @ ns, degree=7)
+  nrg_wall = domain.boundary.integral('σwall dS' @ ns, degree=7)
+  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(δψ σ / ε - η dφ + .5 dt J_k ∇_k(η)) dV' @ ns, degree=7)
 
   numpy.random.seed(seed)
-  state = dict(c=numpy.random.normal(0,.5,ns.cbasis.shape), m=numpy.random.normal(0,.5,ns.mbasis.shape)) # initial condition
+  state = dict(φ=numpy.random.normal(0, .5, φbasis.shape)) # initial condition
 
-  with treelog.iter.plain('timestep', itertools.count()) as steps:
+  with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
    for istep in steps:
 
-    E = function.eval([nrg_mix, nrg_iface, nrg_wall], **state)
-    treelog.user('energy: {0:.3f} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(sum(E), 100*numpy.array(E)/sum(E)))
-
-    x, c, m = bezier.eval(['x_i', 'c', 'm'] @ ns, **state)
-    export.triplot('phase.png', x, c, tri=bezier.tri, clim=(-1,1))
-    export.triplot('chempot.png', x, m, tri=bezier.tri)
-
-    if numpy.ptp(m) < mtol:
-      break
-
-    state['c0'] = state['c']
-    state = solver.optimize(['c', 'm'], nrg, arguments=state, tol=1e-10)
+    E = numpy.stack(function.eval([nrg_mix, nrg_iface, nrg_wall], **state))
+    log.user('energy: {0:,.0μJ/m} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(numpy.sum(E), 100*E/numpy.sum(E)))
+
+    state['φ0'] = state['φ']
+    state = solver.optimize(['φ', 'η'], nrg / tol, arguments=state, tol=1)
+
+    with export.mplfigure('phase.png') as fig:
+      ax = fig.add_subplot(aspect='equal', xlabel='[mm]', ylabel='[mm]')
+      x, φ = bezier.eval(['x_i', 'φ'] @ ns, **state)
+      im = ax.tripcolor(*(x/'mm').T, bezier.tri, φ, shading='gouraud', cmap='bwr')
+      im.set_clim(-1, 1)
+      fig.colorbar(im)
+      if showflux:
+        x, J = grid.eval(['x_i', 'J_i'] @ ns, **state)
+        log.info('largest flux: {:.1mm/h}'.format(numpy.max(numpy.hypot(J[:,0], J[:,1]))))
+        ax.quiver(*(x/'mm').T, *(J/'m/s').T, color='r')
+        ax.quiver(*(x/'mm').T, *-(J/'m/s').T, color='b')
+      ax.autoscale(enable=True, axis='both', tight=True)
 
   return state
 
@@ -114,37 +216,34 @@ def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[f
 class test(testing.TestCase):
 
   def test_initial(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=float('inf'), seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('1h'), seed=0, circle=False, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ0'], '''
       eNoBYgCd/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
       L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJMF4tqQ==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoBYgCd/w7NQModNFjLtckB0VA0rDCiM+zKBMzPygjMcMr4yJcy0MsUyXY0OcoxMFo19zE5Np/JMDTG
-      yk7KGstQzFgwvMnDNXk0Msm+Njc3SjZizebJEjbPy1w25zLsNfQzSjURLRk3Qt4uBQ==''')
 
   def test_square(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/5o2nDZANsQ1dTQO1QXQmTaaNjw2vDVVNOrPqM5mNmY2xDWwNPDRNMsZyys2KjYdNQoyh8v9
-      yffJ1zXTNSA0Ks2JyorJicluNWQ1IDJXy+/JNck3yWQ1WjX0MVHL78k3yTnJbYgt7Q==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNpVyU0KgCAABeGrK9SiQCiQCsrwB8GjzLF8uJNvNYzBsuP5yGKmWlhxXKMSmxzchPGceF4iVU55+Ck0
-      masD28JDNQ==''')
-
-  def test_contactangle(self):
-    state = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=45, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/zs2bjYbNqs19zNqzRXMazaZNkw25zVxNLXOncxDNnU28zUqNeIwTsvKyhc2TDaANdozSMwp
-      yuXJojXiNYo05c/Oyp3Jb8n7NE81iDK2ywfKO8kayXI03jQjLyrLzskZyfrIkb8vTg==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoNzCEOggAYgNFLvbPguILeQRpWgyQPoAYCuBGhWHVzGshOLegMSvffS1/5EqmF2sUnTKJylbO3l6mZ
-      pb2TZ5jLrDWObmGlUDroDWFjq43H3dcvaqdz9TCGP1tYLVU=''')
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='square', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('2h'), seed=0, circle=False, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], '''
+      eNoBYgCd/zE1EjX1NaA2+TXiMxkz0TS9NL01ajaRNZoxYNElNRM1LDUlNZQw0cqgysI1nTWcNN4xLsuk
+      ybDJvDWaNTQ07s7nysnJ6ckPNQY1CzNozKjK58kOysQ0zTQKM73M3coVyjfKR9cuPg==''')
+    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], '''
+      eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
+      xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
 
   def test_mixedcircle(self):
-    state = main(nelems=3, etype='mixed', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=True, stab=stab.linear)
-    with self.subTest('concentration'): self.assertAlmostEqual64(state['c'], '''
-      eNoBYgCd/y415zQQNUw1dzVTNB01fTS4My81dDQpNWU03DJENDo0rjEVMCfMj8vjysUzozMmz3E0VCyg
-      1I00MDMWM3bQU8t4ym7Kd8uKzoTPgc7PzmfOYs3ZzJvM6jIozunL7sqezgXMYsUuJg==''')
-    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['m'], '''
-      eNoljEkKgEAQxP7/B/EkiCCO27gOuJxy8FEGpSEUoaoTkYWTgyAjHRuTZpGzHGgpKchMMz2JRnN/l6j1
-      uctge2W08W8fdlPF5ccXGqBI7Q==''')
+    state = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
+      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
+      etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+      endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
+    with self.subTest('concentration'): self.assertAlmostEqual64(state['φ'], '''
+      eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK
+      z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''')
+    with self.subTest('chemical-potential'): self.assertAlmostEqual64(state['η'], '''
+      eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv
+      x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''')

From 45089267983376caea0d1e34dd92e9bb3f539010 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 30 Dec 2021 23:31:36 +0100
Subject: [PATCH 11/11] add nutils.SI to container image

This match adds the external SI module to the nutils container.
---
 devtools/container/build.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/devtools/container/build.py b/devtools/container/build.py
index 04a9c62f7..49a58902f 100644
--- a/devtools/container/build.py
+++ b/devtools/container/build.py
@@ -57,9 +57,10 @@
   if not examples.exists():
     raise SystemExit(f'cannot find example dir: {examples}')
 
-  container = stack.enter_context(Container.new_from(base, network='none', mounts=[Mount(src=dist, dst='/mnt')]))
+  container = stack.enter_context(Container.new_from(base, mounts=[Mount(src=dist, dst='/mnt')]))
 
   container.run('pip', 'install', '--no-cache-dir', '--no-index', '--find-links', 'file:///mnt/', 'nutils', env=dict(PYTHONHASHSEED='0'))
+  container.run('pip', 'install', '--no-cache-dir', 'https://github.com/evalf/nutils-SI/archive/main.tar.gz', env=dict(PYTHONHASHSEED='0'))
   container.add_label('org.opencontainers.image.url', 'https://github.com/evalf/nutils')
   container.add_label('org.opencontainers.image.source', 'https://github.com/evalf/nutils')
   container.add_label('org.opencontainers.image.authors', 'Evalf')