Merge pull request #328 from gertjanvanzwieten/fixlinear-backport

Fixlinear backport

nutils/function.py

@@ -332,8 +332,7 @@ def __init__(self, ndims:int, trans):
     super().__init__(args=[trans], todims=trans.todims)
 
   def evalf(self, trans):
-    head, tail = transform.promote(trans, self.ndims)
-    return head + tail
+    return transform.promote(trans, self.ndims)
 
 class TailOfTransform(TransformChain):
 

nutils/transform.py

@@ -64,36 +64,40 @@ def canonical(chain):
       i += 1
   return tuple(items)
 
-def promote(chain, ndims):
-  # split chain into chain1 and chain2 such that chain == chain1 << chain2 and
-  # chain1.fromdims == chain2.todims == ndims, where chain1 is canonical and
-  # chain2 climbs to ndims as fast as possible.
+def uppermost(chain):
+  # bring to highest ndims possible
   n = n_ascending(chain)
-  assert ndims >= chain[n-1].fromdims
+  if n < 2:
+    return tuple(chain)
   items = list(chain)
   i = n
-  while items[i-1].fromdims < ndims:
+  while items[i-1].todims < items[0].todims:
     swapped = items[i-2].swapup(items[i-1])
     if swapped:
       items[i-2:i] = swapped
       i += i < n
     else:
       i -= 1
-  assert items[i-1].fromdims == ndims
-  return canonical(items[:i]), tuple(items[i:])
+  return tuple(items)
+
+def promote(chain, ndims):
+  # swap transformations such that ndims is reached as soon as possible, and
+  # then maintained as long as possible (i.e. proceeds as canonical).
+  for i, item in enumerate(chain): # NOTE possible efficiency gain using bisection
+    if item.fromdims == ndims:
+      return canonical(chain[:i+1]) + uppermost(chain[i+1:])
+  return chain # NOTE at this point promotion essentially failed, maybe it's better to raise an exception
 
 def lookup(chain, transforms):
   if not transforms:
     return
   for trans in transforms:
     ndims = trans[-1].fromdims
     break
-  head, tail = promote(chain, ndims)
-  while head:
-    if head in transforms:
-      return head, tail
-    tail = head[-1:] + tail
-    head = head[:-1]
+  chain = promote(chain, ndims)
+  for i in range(len(chain), 0, -1):
+    if chain[i-1].fromdims == ndims and chain[:i] in transforms:
+      return chain[:i], chain[i:]
 
 def lookup_item(chain, transforms):
   head_tail = lookup(chain, transforms)
@@ -103,26 +107,21 @@ def lookup_item(chain, transforms):
   item = transforms[head] if isinstance(transforms, collections.Mapping) else transforms.index(head)
   return item, tail
 
-def linearfrom(chain, ndims):
-  if chain and ndims < chain[-1].fromdims:
-    for i in reversed(range(len(chain))):
-      if chain[i].todims == ndims:
-        chain = chain[:i]
-        break
-    else:
-      raise Exception('failed to find {}D coordinate system'.format(ndims))
+def linearfrom(chain, fromdims):
+  todims = chain[0].todims if chain else fromdims
+  while chain and fromdims < chain[-1].fromdims:
+    chain = chain[:-1]
   if not chain:
-    return numpy.eye(ndims)
+    assert todims == fromdims
+    return numpy.eye(fromdims)
   linear = numpy.eye(chain[-1].fromdims)
-  for trans in reversed(chain):
-    linear = numpy.dot(trans.linear, linear)
-    if trans.todims == trans.fromdims + 1:
-      linear = numpy.concatenate([linear, trans.ext[:,_]], axis=1)
-  n, m = linear.shape
-  if m >= ndims:
-    return linear[:,:ndims]
-  return numpy.concatenate([linear, numpy.zeros((n,ndims-m))], axis=1)
-
+  for transitem in reversed(uppermost(chain)):
+    linear = numpy.dot(transitem.linear, linear)
+    if transitem.todims == transitem.fromdims + 1:
+      linear = numpy.concatenate([linear, transitem.ext[:,_]], axis=1)
+  assert linear.shape[0] == todims
+  return linear[:,:fromdims] if linear.shape[1] >= fromdims \
+    else numpy.concatenate([linear, numpy.zeros((todims, fromdims-linear.shape[1]))], axis=1)
 
 ## TRANSFORM ITEMS
 

tests/test_topology.py

@@ -302,6 +302,30 @@ def test_basistype(self):
   def test_dofcount(self):
     self.assertEqual(len(self.basis), 3*5)
 
+_refined_refs = dict(
+  line=element.LineReference(),
+  quadrilateral=element.LineReference()**2,
+  hexahedron=element.LineReference()**3,
+  triangle=element.TriangleReference(),
+  tetrahedron=element.TetrahedronReference())
+
+@parametrize
+class refined(TestCase):
+
+  def test_boundary_gradient(self):
+    ref = _refined_refs[self.etype]
+    elem = element.Element(ref, (transform.RootTrans('test', (0,)*ref.ndims),))
+    domain = topology.UnstructuredTopology(ref.ndims, (elem,))
+    geom = function.rootcoords(ref.ndims)
+    basis = domain.basis('std', degree=1)
+    u = domain.projection(geom.sum(), onto=basis, geometry=geom, degree=2)
+    g = domain.refine(self.ref1).boundary.refine(self.ref2).elem_eval(u.grad(geom), ischeme='uniform1')
+    numpy.testing.assert_allclose(g, 1)
+
+for etype in _refined_refs:
+  for ref1 in 0, 1:
+    for ref2 in 0, 1:
+      refined(etype=etype, ref1=ref1, ref2=ref2)
 
 @parametrize
 class general(TestCase):