Merge pull request #665 from evalf/transform-updates

Transform updates

nutils/element.py

@@ -778,8 +778,11 @@ def edge_refs(self):
   def getpoints(self, ischeme, degree):
     if ischeme == 'gauss':
       if self.nverts == self.ndims+1: # use optimal gauss schemes for simplex-like cones
-        trans = transform.Square((self.etrans.apply(self.edgeref.vertices) - self.tip).T, self.tip)
-        return points.TransformPoints(getsimplex(self.ndims).getpoints(ischeme, degree), trans)
+        pnts = getsimplex(self.ndims).getpoints(ischeme, degree)
+        linearT = self.etrans.apply(self.edgeref.vertices) - self.tip
+        coords = pnts.coords @ linearT + self.tip
+        weights = pnts.weights * abs(numpy.linalg.det(linearT))
+        return points.CoordsWeightsPoints(coords, weights)
       epoints = self.edgeref.getpoints('gauss', degree)
       tx, tw = points.gauss((degree + self.ndims - 1)//2)
       wx = tx**(self.ndims-1) * tw * self.extnorm * self.height

nutils/function.py

@@ -279,9 +279,9 @@ def normalized(self, __axis: int = -1) -> 'Array':
     'See :func:`normalized`.'
     return normalized(self, __axis)
 
-  def normal(self, exterior: bool = False) -> 'Array':
+  def normal(self, refgeom: Optional['Array'] = None) -> 'Array':
     'See :func:`normal`.'
-    return normal(self, exterior)
+    return normal(self, refgeom)
 
   def curvature(self, ndims: int = -1) -> 'Array':
     'See :func:`curvature`.'
@@ -999,21 +999,16 @@ def lower(self, points_shape: _PointsShape, transform_chains: _TransformChainsMa
 
 class _ExteriorNormal(Array):
 
-  def __init__(self, geom: Array) -> None:
-    assert geom.dtype == float
-    self._geom = geom
-    super().__init__(geom.shape, float, geom.spaces, geom.arguments)
+  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)
 
   def lower(self, points_shape: _PointsShape, transform_chains: _TransformChainsMap, coordinates: _CoordinatesMap) -> evaluable.Array:
-    geom = self._geom.lower(points_shape, transform_chains, coordinates)
-    ref_dim = builtins.sum(transform_chains[space][0].fromdims for space in self._geom.spaces)
-    if self._geom.shape[-1] != ref_dim + 1:
-      raise ValueError('For the exterior normal the dimension of the geometry must be one larger than that of the tip coordinate system, but got {} and {} respectively.'.format(self._geom.shape[-1], ref_dim))
-    refs = tuple((_root_derivative_target if chain.todims == chain.fromdims else _tip_derivative_target)(space, chain.fromdims) for space, (chain, opposite) in transform_chains.items() if space in self._geom.spaces)
-    rgrad = evaluable.concatenate([evaluable.derivative(geom, ref) for ref in refs], axis=-1)
-    if self._geom.shape[-1] == 2:
+    rgrad = self._rgrad.lower(points_shape, transform_chains, coordinates)
+    if self._rgrad.shape[-2] == 2:
       normal = evaluable.stack([rgrad[...,1,0], -rgrad[...,0,0]], axis=-1)
-    elif self._geom.shape[-1] == 3:
+    elif self._rgrad.shape[-2] == 3:
       i = evaluable.asarray([1, 2, 0])
       j = evaluable.asarray([2, 0, 1])
       normal = evaluable.Take(rgrad[...,0], i) * evaluable.Take(rgrad[...,1], j) - evaluable.Take(rgrad[...,1], i) * evaluable.Take(rgrad[...,0], j)
@@ -2789,13 +2784,17 @@ def curl(__arg: IntoArray, __geom: IntoArray) -> Array:
     raise ValueError('Expected a function with a trailing axis of length 3 but got {}.'.format(arg.shape[-1]))
   return (levicivita(3).T * _append_axes(grad(arg, geom), (3,))).sum((-3, -2))
 
-def normal(__geom: IntoArray, exterior: bool = False) -> Array:
+def normal(__geom: IntoArray, refgeom: Optional[Array] = None) -> Array:
   '''Return the normal of the geometry.
 
   Parameters
   ----------
   geom : :class:`Array` or something that can be :meth:`~Array.cast` into one
-  exterior : :class:`bool`
+  refgeom : :class:`Array`, optional`
+      The reference geometry. If ``None``, the reference geometry is the tip
+      coordinate system of the spaces on which ``geom`` is defined. The
+      dimension of the reference geometry must be exactly one smaller than the
+      dimension of the geometry.
 
   Returns
   -------
@@ -2806,14 +2805,18 @@ def normal(__geom: IntoArray, exterior: bool = False) -> Array:
   if geom.dtype != float:
     raise ValueError('The geometry must be real-valued.')
   if geom.ndim == 0:
-    return normal(insertaxis(geom, 0, 1), exterior)[...,0]
+    return normal(insertaxis(geom, 0, 1), refgeom)[...,0]
   elif geom.ndim > 1:
     sh = geom.shape[-2:]
-    return unravel(normal(ravel(geom, geom.ndim-2), exterior), geom.ndim-2, sh)
-  elif not exterior:
+    return unravel(normal(ravel(geom, geom.ndim-2), refgeom), geom.ndim-2, sh)
+  elif refgeom is None:
     return _Normal(geom)
+  elif refgeom.dtype != float:
+    raise ValueError('The reference geometry must be real-valued.')
+  elif refgeom.shape != (geom.shape[0]-1,):
+    raise ValueError('The reference geometry must have shape ({},) but got {}.'.format(geom.shape[0]-1, refgeom.shape))
   else:
-    return _ExteriorNormal(geom)
+    return _ExteriorNormal(grad(geom, refgeom))
 
 def dotnorm(__arg: IntoArray, __geom: IntoArray, axis: int = -1) -> Array:
   '''Return the inner product of an array with the normal of the given geometry.

nutils/mesh.py

@@ -28,87 +28,72 @@
 
 from . import topology, function, util, element, numeric, transform, transformseq, warnings, types, cache
 from .elementseq import References
+from .transform import TransformItem
+from .topology import Topology
+from typing import Optional, Sequence, Tuple, Union
 import numpy, os, itertools, re, math, treelog as log, io, contextlib
 _ = numpy.newaxis
 
 # MESH GENERATORS
 
 @log.withcontext
-def rectilinear(richshape, periodic=(), name='rect', space='X'):
+def rectilinear(richshape: Sequence[Union[int, Sequence[float]]], periodic: Sequence[int] = (), name: Optional[str] = None, space: str = 'X', root: Optional[TransformItem] = None) -> Tuple[Topology, function.Array]:
   'rectilinear mesh'
 
-  ndims = len(richshape)
-  shape = []
-  offset = []
-  scale = []
-  uniform = True
-  for v in richshape:
-    if numeric.isint(v):
-      assert v > 0
-      shape.append(v)
-      scale.append(1)
-      offset.append(0)
-    elif numpy.equal(v, numpy.linspace(v[0],v[-1],len(v))).all():
-      shape.append(len(v)-1)
-      scale.append((v[-1]-v[0]) / float(len(v)-1))
-      offset.append(v[0])
-    else:
-      shape.append(len(v)-1)
-      uniform = False
-
-  root = transform.Identifier(ndims, name)
+  verts = [numpy.arange(v + 1) if numeric.isint(v) else v for v in richshape]
+  shape = [len(v) - 1 for v in verts]
+  ndims = len(shape)
+
+  if name is not None:
+    warnings.deprecation('Argument `name` is deprecated; use `root` with a `TransformItem` instead.')
+    if root is not None:
+      raise ValueError('Arguments `name` and `root` cannot be used simultaneously.')
+    root = transform.Index(hash(name))
+  elif root is None:
+    root = transform.Index(ndims, 0)
+
   axes = [transformseq.DimAxis(i=0, j=n, mod=n if idim in periodic else 0, isperiodic=idim in periodic) for idim, n in enumerate(shape)]
   topo = topology.StructuredTopology(space, root, axes)
 
-  if uniform:
-    if all(o == offset[0] for o in offset[1:]):
-      offset = offset[0]
-    if all(s == scale[0] for s in scale[1:]):
-      scale = scale[0]
-    geom = function.rootcoords(space, ndims) * scale + offset
-  else:
-    funcsp = topo.basis('spline', degree=1, periodic=())
-    coords = numeric.meshgrid(*richshape).reshape(ndims, -1)
-    geom = (funcsp * coords).sum(-1)
+  funcsp = topo.basis('spline', degree=1, periodic=())
+  coords = numeric.meshgrid(*verts).reshape(ndims, -1)
+  geom = (funcsp * coords).sum(-1)
 
   return topo, geom
 
 _oldrectilinear = rectilinear # reference for internal unittests
 
-def line(nodes, periodic=False, bnames=None, *, name: str = 'line', space: str = 'X'):
+def line(nodes: Union[int, Sequence[float]], periodic: bool = False, bnames: Optional[Sequence[Tuple[str, str]]] = None, *, name: Optional[str] = None, space: str = 'X', root: Optional[TransformItem] = None) -> Tuple[Topology, function.Array]:
+  if name is not None:
+    warnings.deprecation('Argument `name` is deprecated; use `root` with a `transform.transformitem` instead.')
+    if root is not None:
+      raise ValueError('Arguments `name` and `root` cannot be used simultaneously.')
+    root = transform.Index(hash(name))
+  elif root is None:
+    root = transform.Index(1, 0)
   if isinstance(nodes, int):
-    uniform = True
-    assert nodes > 0
-    nelems = nodes
-    scale = 1
-    offset = 0
-  else:
-    nelems = len(nodes)-1
-    scale = (nodes[-1]-nodes[0]) / nelems
-    offset = nodes[0]
-    uniform = numpy.equal(nodes, offset + numpy.arange(nelems+1) * scale).all()
-  root = transform.Identifier(1, name)
-  domain = topology.StructuredLine(space, root, 0, nelems, periodic=periodic, bnames=bnames)
-  geom = function.rootcoords(space, 1)[0] * scale + offset if uniform else domain.basis('std', degree=1, periodic=[]).dot(nodes)
+    nodes = numpy.arange(nodes + 1)
+  domain = topology.StructuredLine(space, root, 0, len(nodes) - 1, periodic=periodic, bnames=bnames)
+  geom = domain.basis('std', degree=1, periodic=[]).dot(nodes)
   return domain, geom
 
-def newrectilinear(nodes, periodic=None, name='rect', bnames=[['left','right'],['bottom','top'],['front','back']], spaces=None):
+def newrectilinear(nodes: Sequence[Union[int, Sequence[float]]], periodic: Optional[Sequence[int]] = None, name: Optional[str] = None, bnames=[['left','right'],['bottom','top'],['front','back']], spaces: Optional[Sequence[str]] = None, root: Optional[TransformItem] = None) -> Tuple[Topology, function.Array]:
   if periodic is None:
     periodic = []
   if not spaces:
     spaces = 'XYZ' if len(nodes) <= 3 else map('R{}'.format, range(len(nodes)))
   else:
     assert len(spaces) == len(nodes)
-  domains, geoms = zip(*(line(nodesi, i in periodic, bnamesi, name=name, space=spacei) for i, (nodesi, bnamesi, spacei) in enumerate(zip(nodes, tuple(bnames)+(None,)*len(nodes), spaces))))
+  domains, geoms = zip(*(line(nodesi, i in periodic, bnamesi, name=name, space=spacei, root=root) for i, (nodesi, bnamesi, spacei) in enumerate(zip(nodes, tuple(bnames)+(None,)*len(nodes), spaces))))
   return util.product(domains), function.stack(geoms)
 
 if os.environ.get('NUTILS_TENSORIAL'):
-  def rectilinear(richshape, periodic=(), name='rect', space='X'):
+  def rectilinear(richshape: Sequence[Union[int, Sequence[float]]], periodic: Sequence[int] = (), name: Optional[str] = None, space: str = 'X', root: Optional[TransformItem] = None) -> Tuple[Topology, function.Array]:
     spaces = tuple(space+str(i) for i in range(len(richshape)))
-    return newrectilinear(richshape, periodic, name, spaces=spaces)
+    return newrectilinear(richshape, periodic, name=name, spaces=spaces, root=root)
 
 @log.withcontext
-def multipatch(patches, nelems, patchverts=None, name='multipatch'):
+def multipatch(patches, nelems, patchverts=None, name: Optional[str] = None):
   '''multipatch rectilinear mesh generator
 
   Generator for a :class:`~nutils.topology.MultipatchTopology` and geometry.
@@ -194,32 +179,17 @@ def multipatch(patches, nelems, patchverts=None, name='multipatch'):
   >>> from nutils import function
   >>> topo, cube = multipatch(
   ...   patches=[
-  ...     [0,1,2,3], # left,   normal: x
-  ...     [4,5,6,7], # right,  normal: x
-  ...     [0,1,4,5], # bottom, normal: -y
-  ...     [2,3,6,7], # top,    normal: -y
-  ...     [0,2,4,6], # front,  normal: z
-  ...     [1,3,5,7], # back,   normal: z
+  ...     [0,1,2,3], # x=-1
+  ...     [4,5,6,7], # x= 1
+  ...     [0,1,4,5], # y=-1
+  ...     [2,3,6,7], # y= 1
+  ...     [0,2,4,6], # z=-1
+  ...     [1,3,5,7], # z= 1
   ...   ],
   ...   patchverts=tuple(itertools.product(*([[-1,1]]*3))),
   ...   nelems=1)
   >>> sphere = function.normalized(cube)
 
-  The normals of the patches are determined by the order of the vertex numbers.
-  An outward normal for the cube is obtained by flipping the left, top and
-  front faces:
-
-  >>> cubenormal = cube.normal(exterior=True) * topo.basis('patch').dot([-1,1,1,-1,-1,1])
-
-  At the centroids of the faces the outward normal should equal the cube geometry:
-
-  >>> numpy.testing.assert_allclose(*topo.sample('gauss', 1).eval([cubenormal, cube]))
-
-  Similarly, the outward normal of the sphere is obtained by:
-
-  >>> spherenormal = sphere.normal(exterior=True) * topo.basis('patch').dot([-1,1,1,-1,-1,1])
-  >>> numpy.testing.assert_allclose(*topo.sample('gauss', 1).eval([spherenormal, cube]))
-
   Args
   ----
   patches:
@@ -241,6 +211,9 @@ def multipatch(patches, nelems, patchverts=None, name='multipatch'):
       if ``patchverts`` is not specified.
   '''
 
+  if name is not None:
+    warnings.deprecation('Argument `name` is deprecated and can safely be omitted.')
+
   patches = numpy.array(patches)
   if patches.dtype != int:
     raise ValueError('`patches` should be an array of ints.')
@@ -300,7 +273,7 @@ def multipatch(patches, nelems, patchverts=None, name='multipatch'):
         raise ValueError('duplicate number of elements specified for patch {} in dimension {}'.format(i, dim))
       shape.append(nelems_sides[0])
     # create patch topology
-    topos.append(rectilinear(shape, name='{}{}'.format(name, i))[0])
+    topos.append(rectilinear(shape, root=transform.Index(ndims, i))[0])
     # compute patch geometry
     patchcoords = [numpy.linspace(0, 1, n+1) for n in shape]
     patchcoords = numeric.meshgrid(*patchcoords).reshape(ndims, -1)
@@ -549,7 +522,7 @@ def simplex(nodes, cnodes, coords, tags, btags, ptags, name='simplex', *, space=
   assert numpy.greater(nodes[:,1:], nodes[:,:-1]).all(), 'nodes must be sorted'
   assert ncnodes == _comb(ndims + degree, degree), 'number of coordinate nodes does not correspond to uniformly refined simplex'
 
-  transforms = transformseq.IdentifierTransforms(ndims=ndims, name=name, length=nelems)
+  transforms = transformseq.IndexTransforms(ndims=ndims, length=nelems)
   topo = topology.SimplexTopology(space, nodes, transforms, transforms)
   coeffs = element.getsimplex(ndims).get_poly_coeffs('lagrange', degree=degree)
   basis = function.PlainBasis([coeffs] * nelems, cnodes, nverts, topo.f_index, topo.f_coords)
@@ -577,7 +550,7 @@ def simplex(nodes, cnodes, coords, tags, btags, ptags, name='simplex', *, space=
 
   pgroups = {}
   if ptags:
-    ptrans = [transform.Matrix(linear=numpy.zeros(shape=(ndims,0)), offset=offset) for offset in numpy.eye(ndims+1)[:,1:]]
+    ptrans = [transform.Point(offset) for offset in numpy.eye(ndims+1)[:,1:]]
     pmap = {inode: numpy.array(numpy.equal(nodes, inode).nonzero()).T for inode in set.union(*map(set, ptags.values()))}
     for pname, inodes in ptags.items():
       ptransforms = transformseq.PlainTransforms([topo.transforms[ielem] + (ptrans[ivertex],) for inode in inodes for ielem, ivertex in pmap[inode]], ndims, 0)
@@ -669,38 +642,40 @@ def unitsquare(nelems, etype):
   space = 'X'
 
   if etype == 'square':
-    topo, geom = rectilinear([nelems, nelems], name='unitsquare', space=space)
+    topo, geom = rectilinear([nelems, nelems], space=space)
     return topo, geom / nelems
 
   elif etype in ('triangle', 'mixed'):
-    root = transform.Identifier(2, 'unitsquare')
     simplices = numpy.concatenate([
       numpy.take([i*(nelems+1)+j, i*(nelems+1)+j+1, (i+1)*(nelems+1)+j, (i+1)*(nelems+1)+j+1], [[0,1,2],[1,2,3]] if i%2==j%2 else [[0,1,3],[0,2,3]], axis=0)
         for i in range(nelems) for j in range(nelems)])
 
     v = numpy.arange(nelems+1, dtype=float)
     coords = numeric.meshgrid(v, v).reshape(2,-1).T
-    transforms = transformseq.PlainTransforms([(root, transform.Square((c[1:]-c[0]).T, c[0])) for c in coords[simplices]], 2, 2)
+    transforms = transformseq.IndexTransforms(2, len(simplices))
     topo = topology.SimplexTopology(space, simplices, transforms, transforms)
 
     if etype == 'mixed':
       references = list(topo.references)
-      transforms = list(topo.transforms)
       square = element.getsimplex(1)**2
       connectivity = list(topo.connectivity)
       isquares = [i * nelems + j for i in range(nelems) for j in range(nelems) if i%2==j%3]
+      dofs = list(simplices)
       for n in sorted(isquares, reverse=True):
         i, j = divmod(n, nelems)
         references[n*2:(n+1)*2] = square,
-        transforms[n*2:(n+1)*2] = (root, transform.Shift([float(i),float(j)])),
         connectivity[n*2:(n+1)*2] = numpy.concatenate(connectivity[n*2:(n+1)*2])[[3,2,4,1] if i%2==j%2 else [3,2,0,5]],
         connectivity = [c-numpy.greater(c,n*2) for c in connectivity]
-      topo = topology.ConnectedTopology(space, References.from_iter(references, 2), transformseq.PlainTransforms(transforms, 2, 2), transformseq.PlainTransforms(transforms, 2, 2), connectivity)
-
-    x, y = topo.boundary.sample('_centroid', None).eval(function.rootcoords(space, 2)).T
-    bgroups = dict(left=x==0, right=x==nelems, bottom=y==0, top=y==nelems)
+        dofs[n*2:(n+1)*2] = numpy.unique([*dofs[n*2], *dofs[n*2+1]]),
+      coords = coords[numpy.argsort(numpy.unique(numpy.concatenate(dofs), return_index=True)[1])]
+      transforms = transformseq.IndexTransforms(2, len(connectivity))
+      topo = topology.ConnectedTopology(space, References.from_iter(references, 2), transforms, transforms, connectivity)
+
+    geom = (topo.basis('std', degree=1) * coords.T).sum(-1)
+    x, y = topo.boundary.sample('_centroid', None).eval(geom).T
+    bgroups = dict(left=x<.1, right=x>nelems-.1, bottom=y<.1, top=y>nelems-.1)
     topo = topo.withboundary(**{name: topo.boundary[numpy.where(mask)[0]] for name, mask in bgroups.items()})
-    return topo, function.rootcoords(space, 2) / nelems
+    return topo, geom / nelems
 
   else:
     raise Exception('invalid element type {!r}'.format(etype))

nutils/topology.py

@@ -713,7 +713,7 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
     >>> from . import mesh
     >>> domain, geom = mesh.unitsquare(nelems=3, etype='mixed')
     >>> sample = domain.locate(geom, [[.9, .4]], tol=1e-12)
-    >>> sample.eval(geom).tolist()
+    >>> sample.eval(geom).round(5).tolist()
     [[0.9, 0.4]]
 
     Locate requires a geometry function, an array of coordinates, and at least
@@ -2230,11 +2230,18 @@ def locate(self, geom, coords, *, tol, eps=0, weights=None, **kwargs):
     return self._locate(geom0, scale, coords, eps=eps, weights=weights)
 
   def _asaffine(self, geom):
-    index = function.rootcoords(self.space, self.transforms.todims)[[axis.isdim for axis in self.axes]] * 2**self.nrefine - [axis.i for axis in self.axes if axis.isdim]
-    basis = function.concatenate([function.eye(self.ndims), function.diagonalize(index)], axis=0)
-    A, b = map(sparse.toarray, self.sample('gauss', 2).integrate_sparse([(basis[:,_,:] * basis[_,:,:]).sum(-1), (basis * geom).sum(-1)]))
-    x = numpy.linalg.solve(A, b)
-    return x[:self.ndims], x[self.ndims:], index
+    p0 = p1 = self
+    for (b0, b1), axis in zip(self._bnames, self.axes):
+      if axis.isdim:
+        p0 = p0[:].boundary[b0]
+        p1 = p1[:].boundary[b1]
+    geom0, = p0.sample('gauss', 0).eval(geom)
+    geom1, = p1.sample('gauss', 0).eval(geom)
+    funcsp = self.basis('std', degree=1, periodic=())
+    verts = numeric.meshgrid(*map(numpy.arange, numpy.array(self.shape)+1)).reshape(self.ndims, -1)
+    index = (funcsp * verts).sum(-1)
+    scale = (geom1 - geom0) / numpy.array(self.shape)
+    return geom0, scale, index
 
   def _locate(self, geom0, scale, coords, *, eps=0, weights=None):
     mincoords, maxcoords = numpy.sort([geom0, geom0 + scale * self.shape], axis=0)