diff --git a/nutils/element.py b/nutils/element.py index 9071204c6..7b4a8cad2 100644 --- a/nutils/element.py +++ b/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 diff --git a/nutils/function.py b/nutils/function.py index 4441e90ec..d37f4fea6 100644 --- a/nutils/function.py +++ b/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. diff --git a/nutils/mesh.py b/nutils/mesh.py index 6f58cf421..f49327d2b 100644 --- a/nutils/mesh.py +++ b/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)) diff --git a/nutils/topology.py b/nutils/topology.py index 553975c1a..66b3fac5a 100644 --- a/nutils/topology.py +++ b/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) diff --git a/nutils/transform.py b/nutils/transform.py index 19f1da432..e6e6e8873 100644 --- a/nutils/transform.py +++ b/nutils/transform.py @@ -212,36 +212,7 @@ def transform_poly(self, coeffs): self._transform_matrix[degree] = M return types.frozenarray(numpy.einsum('jk,ik', M.reshape([(degree+1)**self.fromdims]*2), coeffs.reshape(coeffs.shape[0],-1)).reshape(coeffs.shape), copy=False) -class Shift(Square): - '''Shift transformation :math:`x ↦ x + b` - - Parameters - ---------- - offset : :class:`numpy.ndarray` - The offset :math:`b`. - ''' - - __slots__ = () - - det = 1. - - @types.apply_annotations - def __init__(self, offset:types.arraydata): - assert offset.ndim == 1 and offset.dtype == float - super().__init__(numpy.eye(offset.shape[0]), offset) - - @types.lru_cache - def apply(self, points): - return types.frozenarray(points + self.offset, copy=False) - - @types.lru_cache - def invapply(self, points): - return types.frozenarray(points - self.offset, copy=False) - - def __str__(self): - return '{}+x'.format(util.obj2str(self.offset)) - -class Identity(Shift): +class Identity(Square): '''Identity transformation :math:`x ↦ x` Parameters @@ -252,8 +223,10 @@ class Identity(Shift): __slots__ = () + det = 1. + def __init__(self, ndims): - super().__init__(numpy.zeros(ndims)) + super().__init__(numpy.eye(ndims), numpy.zeros(ndims)) def apply(self, points): return points @@ -264,45 +237,23 @@ def invapply(self, points): def __str__(self): return 'x' -class Scale(Square): - '''Affine transformation :math:`x ↦ a x + b`, with :math:`a` a scalar +class Index(Identity): + '''Identity transform with index - Parameters - ---------- - scale : :class:`float` - The scalar :math:`a`. - offset : :class:`numpy.ndarray` - The offset :math:`b`. + This transformation serves as an element-specific or topology-specific index + to form the basis of transformation lookups. Otherwise, the transform behaves + like an identity. ''' - __slots__ = 'scale', + __slots__ = 'index' @types.apply_annotations - def __init__(self, scale:float, offset:types.arraydata): - assert offset.ndim == 1 and offset.dtype == float - self.scale = scale - super().__init__(numpy.eye(offset.shape[0]) * scale, offset) - - @types.lru_cache - def apply(self, points): - return types.frozenarray(self.scale * points + self.offset, copy=False) - - @types.lru_cache - def invapply(self, points): - return types.frozenarray((points - self.offset) / self.scale, copy=False) - - @property - def det(self): - return self.scale**self.todims - - def __str__(self): - return '{}+{}*x'.format(util.obj2str(self.offset), self.scale) + def __init__(self, ndims:int, index:int): + self.index = index + super().__init__(ndims) - def __mul__(self, other): - assert isinstance(other, Matrix) and self.fromdims == other.todims - if isinstance(other, Scale): - return Scale(self.scale * other.scale, self.apply(other.offset)) - return super().__mul__(other) + def __repr__(self): + return 'Index({}, {})'.format(self.todims, self.index) class Updim(Matrix): '''Affine transformation :math:`x ↦ A x + b`, with :math:`A` an :math:`n×(n-1)` matrix @@ -512,23 +463,11 @@ def __init__(self, trans1, trans2): def det(self): return self.trans1.det * self.trans2.det -class Identifier(Identity): - '''Generic identifier - - This transformation serves as an element-specific or topology-specific token - to form the basis of transformation lookups. Otherwise, the transform behaves - like an identity. - ''' - - __slots__ = 'token' +class Point(Matrix): @types.apply_annotations - def __init__(self, ndims:int, token): - self.token = token - super().__init__(ndims) - - def __str__(self): - return ':'.join(map(str, self._args)) + def __init__(self, offset:types.arraydata): + super().__init__(numpy.zeros((offset.shape[0],0)), offset) ## EVALUABLE TRANSFORM CHAIN diff --git a/nutils/transformseq.py b/nutils/transformseq.py index 32f0c96c8..17d88a4c8 100644 --- a/nutils/transformseq.py +++ b/nutils/transformseq.py @@ -146,21 +146,21 @@ def index_with_tail(self, trans): Example ------- - Consider the following plain sequence of two shift transforms: + Consider the following plain sequence of two index transforms: - >>> from nutils.transform import Shift, Scale - >>> transforms = PlainTransforms([(Shift([0.]),), (Shift([1.]),)], 1, 1) + >>> from nutils.transform import Index, SimplexChild + >>> transforms = PlainTransforms([(Index(1, 0),), (Index(1, 1),)], 1, 1) Calling :meth:`index_with_tail` with the first transform gives index ``0`` and no tail: - >>> transforms.index_with_tail((Shift([0.]),)) + >>> transforms.index_with_tail((Index(1, 0),)) (0, ()) Calling with an additional scale gives: - >>> transforms.index_with_tail((Shift([0.]), Scale(0.5, [0.]))) - (0, (Scale([0]+0.5*x),)) + >>> transforms.index_with_tail((Index(1, 0), SimplexChild(1, 0))) + (0, (SimplexChild([0]+[.5]*x0),)) ''' raise NotImplementedError @@ -191,23 +191,23 @@ def index(self, trans): Example ------- - Consider the following plain sequence of two shift transforms: + Consider the following plain sequence of two index transforms: - >>> from nutils.transform import Shift, Scale - >>> transforms = PlainTransforms([(Shift([0.]),), (Shift([1.]),)], 1, 1) + >>> from nutils.transform import Index, SimplexChild + >>> transforms = PlainTransforms([(Index(1, 0),), (Index(1, 1),)], 1, 1) Calling :meth:`index` with the first transform gives index ``0``: - >>> transforms.index((Shift([0.]),)) + >>> transforms.index((Index(1, 0),)) 0 Calling with an additional scale raises an exception, because the transform is not present in ``transforms``. - >>> transforms.index((Shift([0.]), Scale(0.5, [0.]))) + >>> transforms.index((Index(1, 0), SimplexChild(1, 0))) Traceback (most recent call last): ... - ValueError: (Shift([0]+x), Scale([0]+0.5*x)) not in sequence of transforms + ValueError: (Index(1, 0), SimplexChild([0]+[.5]*x0)) not in sequence of transforms ''' index, tail = self.index_with_tail(trans) @@ -419,46 +419,42 @@ def index_with_tail(self, trans): raise ValueError('{!r} not in sequence of transforms'.format(orig_trans)) return self._indices[i], trans[len(match):] -class IdentifierTransforms(Transforms): - '''A sequence of :class:`nutils.transform.Identifier` singletons. - - Every identifier is instantiated with three arguments: the dimension, the - name string, and an integer index matching its position in the sequence. +class IndexTransforms(Transforms): + '''A sequence of :class:`nutils.transform.Index` singletons. Parameters ---------- ndims : :class:`int` Dimension of the transformation. - name : :class:`str` - Identifying name string. length : :class:`int` Length of the sequence. + offset : :class:`int` + The index of the first :class:`nutils.transform.Index` in this sequence. ''' - __slots__ = '_name', '_length' + __slots__ = '_length', '_offset' @types.apply_annotations - def __init__(self, ndims:types.strictint, name:str, length:int): - self._name = name + def __init__(self, ndims:types.strictint, length:int, offset:int = 0): self._length = length + self._offset = offset super().__init__(ndims, ndims) def __getitem__(self, index): if not numeric.isint(index): return super().__getitem__(index) - index = int(index) # make sure that index is a Python integer rather than numpy.intxx - return transform.Identifier(self.fromdims, (self._name, numeric.normdim(self._length, index))), + return transform.Index(self.fromdims, self._offset + numeric.normdim(self._length, index.__index__())), def get_evaluable(self, index: evaluable.Array) -> EvaluableTransformChain: - return _EvaluableIdentifierChain(self.fromdims, self._name, evaluable.InRange(index, self._length)) + return _EvaluableIndexChain(self.fromdims, self._offset + evaluable.InRange(index, self._length)) def __len__(self): return self._length def index_with_tail(self, trans): root = trans[0] - if root.fromdims == self.fromdims and isinstance(root, transform.Identifier) and isinstance(root.token, tuple) and len(root.token) == 2 and root.token[0] == self._name and 0 <= root.token[1] < self._length: - return root.token[1], trans[1:] + if root.fromdims == self.fromdims and isinstance(root, transform.Index) and 0 <= root.index - self._offset < self._length: + return root.index - self._offset, trans[1:] raise ValueError class Axis(types.Singleton): @@ -598,23 +594,23 @@ def __getitem__(self, index): for i in range(self._nrefine): indices, r = divmod(indices, self._ctransforms.shape) ctransforms.insert(0, self._ctransforms[tuple(r)]) - trans0 = transform.Shift(types.frozenarray(indices, dtype=float, copy=False)) - return (self._root, trans0, *ctransforms, *self._etransforms) + trans0 = (transform.Index(len(self._axes), index) for index in indices) + return (self._root, *trans0, *ctransforms, *self._etransforms) def __len__(self): return util.product(map(len, self._axes)) def index_with_tail(self, trans): - if len(trans) < 2 + self._nrefine + len(self._etransforms): + if len(trans) < 1 + len(self._axes) + self._nrefine + len(self._etransforms): raise ValueError - root, shift, tail = trans[0], trans[1], transform.uppermost(trans[2:]) + root, indices, tail = trans[0], trans[1:1+len(self._axes)], transform.uppermost(trans[1+len(self._axes):]) if root != self._root: raise ValueError - if not isinstance(shift, transform.Shift) or len(shift.offset) != len(self._axes) or not numpy.equal(shift.offset.astype(int), shift.offset).all(): + if not all(isinstance(index, transform.Index) and index.todims == len(self._axes) for index in indices): raise ValueError - indices = numpy.array(shift.offset, dtype=int) + indices = numpy.array([index.index for index in indices], dtype=int) # Match child transforms. for item in tail[:self._nrefine]: @@ -970,17 +966,16 @@ def index_with_tail_in(self, __sequence) -> Tuple[evaluable.Array, EvaluableTran else: return super().index_with_tail_in(__sequence) -class _EvaluableIdentifierChain(EvaluableTransformChain): +class _EvaluableIndexChain(EvaluableTransformChain): - __slots__ = '_ndim', '_name' + __slots__ = '_ndim' - def __init__(self, ndim: int, name: str, index: evaluable.Array) -> None: + def __init__(self, ndim: int, index: evaluable.Array) -> None: self._ndim = ndim - self._name = name super().__init__((index,), ndim, ndim) def evalf(self, index: numpy.ndarray) -> TransformChain: - return transform.Identifier(self._ndim, (self._name, index.__index__())), + return transform.Index(self._ndim, index.__index__()), def apply(self, points: evaluable.Array) -> evaluable.Array: return points diff --git a/tests/test_basis.py b/tests/test_basis.py index fdc1abd8a..e3cc014dd 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -268,10 +268,9 @@ def setUp(self): simplices = numeric.overlapping(numpy.arange(nverts), n=self.ndims+1) coords = numpy.random.normal(size=(nverts, self.ndims)) space = 'test' - root = transform.Identifier(self.ndims, 'test') - transforms = transformseq.PlainTransforms([(root, transform.Square((c[1:]-c[0]).T, c[0])) for c in coords[simplices]], self.ndims, self.ndims) + transforms = transformseq.IndexTransforms(self.ndims, len(simplices)) domain = topology.SimplexTopology(space, simplices, transforms, transforms) - geom = function.rootcoords(space, self.ndims) + geom = (domain.basis('std', degree=1) * coords.T).sum(-1) else: raise NotImplementedError self.domain = domain diff --git a/tests/test_expression_v2.py b/tests/test_expression_v2.py index d4db5a5f1..700b2504d 100644 --- a/tests/test_expression_v2.py +++ b/tests/test_expression_v2.py @@ -546,7 +546,7 @@ def test_define_for_3d(self): ns.ε = function.levicivita(3) ns.f = function.Array.cast([['x', '-z', 'y'], ['0', 'x z', '0']] @ ns) smpl = topo.sample('gauss', 5) - assertEvalAlmostEqual = lambda *args: self.assertAllAlmostEqual(*(smpl(f).as_evaluable_array.eval() for f in args)) + assertEvalAlmostEqual = lambda *args: self.assertAllAlmostEqual(*(smpl(f).as_evaluable_array.simplified.eval() for f in args)) assertEvalAlmostEqual('curl_ij(y δ_j0 - x δ_j1 + z δ_j2)' @ ns, '-2 δ_i2' @ ns) assertEvalAlmostEqual('curl_ij(-x^2 δ_j1)' @ ns, '-2 x δ_i2' @ ns) assertEvalAlmostEqual('curl_ij((x δ_j0 - z δ_j1 + y δ_j2) δ_k0 + x z δ_j1 δ_k1)' @ ns, '2 δ_i0 δ_k0 - x δ_i0 δ_k1 + z δ_i2 δ_k1' @ ns) diff --git a/tests/test_function.py b/tests/test_function.py index 6b3a952bc..94ef4e309 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -1138,7 +1138,7 @@ def test_lower(self): class PlainBasis(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.PlainTransforms([(transform.Identifier(0,k),) for k in 'abcd'], 0, 0) + self.checktransforms = transformseq.IndexTransforms(0, 4) index, coords = self.mk_index_coords(0, self.checktransforms) self.checkcoeffs = [[1.],[2.,3.],[4.,5.],[6.]] self.checkdofs = [[0],[2,3],[1,3],[2]] @@ -1149,7 +1149,7 @@ def setUp(self): class DiscontBasis(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.PlainTransforms([(transform.Identifier(0,k),) for k in 'abcd'], 0, 0) + self.checktransforms = transformseq.IndexTransforms(0, 4) index, coords = self.mk_index_coords(0, self.checktransforms) self.checkcoeffs = [[1.],[2.,3.],[4.,5.],[6.]] self.basis = function.DiscontBasis(self.checkcoeffs, index, coords) @@ -1160,7 +1160,7 @@ def setUp(self): class LegendreBasis(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.IdentifierTransforms(1, 'test', 3) + self.checktransforms = transformseq.IndexTransforms(1, 3) index, coords = self.mk_index_coords(0, self.checktransforms) self.checkcoeffs = [[[1,0,0,0],[-1,2,0,0],[1,-6,6,0],[-1,12,-30,20]]]*3 self.basis = function.LegendreBasis(3, 3, index, coords) @@ -1172,7 +1172,7 @@ def setUp(self): class MaskedBasis(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.PlainTransforms([(transform.Identifier(0,k),) for k in 'abcd'], 0, 0) + self.checktransforms = transformseq.IndexTransforms(0, 4) index, coords = self.mk_index_coords(0, self.checktransforms) parent = function.PlainBasis([[1.],[2.,3.],[4.,5.],[6.]], [[0],[2,3],[1,3],[2]], 4, index, coords) self.basis = function.MaskedBasis(parent, [0,2]) @@ -1184,7 +1184,7 @@ def setUp(self): class PrunedBasis(CommonBasis, TestCase): def setUp(self): - parent_transforms = transformseq.PlainTransforms([(transform.Identifier(0,k),) for k in 'abcd'], 0, 0) + parent_transforms = transformseq.IndexTransforms(0, 4) parent_index, parent_coords = self.mk_index_coords(0, parent_transforms) indices = types.frozenarray([0,2]) self.checktransforms = parent_transforms[indices] @@ -1199,7 +1199,7 @@ def setUp(self): class StructuredBasis1D(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.StructuredTransforms(transform.Identifier(1, 'test'), [transformseq.DimAxis(0,4,0,False)], 0) + self.checktransforms = transformseq.IndexTransforms(1, 4) index, coords = self.mk_index_coords(1, self.checktransforms) self.basis = function.StructuredBasis([[[[1],[2]],[[3],[4]],[[5],[6]],[[7],[8]]]], [[0,1,2,3]], [[2,3,4,5]], [5], [4], index, coords) self.checkcoeffs = [[[1.],[2.]],[[3.],[4.]],[[5.],[6.]],[[7.],[8.]]] @@ -1210,7 +1210,7 @@ def setUp(self): class StructuredBasis1DPeriodic(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.StructuredTransforms(transform.Identifier(1, 'test'), [transformseq.DimAxis(0,4,4,True)], 0) + self.checktransforms = transformseq.IndexTransforms(1, 4) index, coords = self.mk_index_coords(1, self.checktransforms) self.basis = function.StructuredBasis([[[[1],[2]],[[3],[4]],[[5],[6]],[[7],[8]]]], [[0,1,2,3]], [[2,3,4,5]], [4], [4], index, coords) self.checkcoeffs = [[[1.],[2.]],[[3.],[4.]],[[5.],[6.]],[[7.],[8.]]] @@ -1221,7 +1221,7 @@ def setUp(self): class StructuredBasis2D(CommonBasis, TestCase): def setUp(self): - self.checktransforms = transformseq.StructuredTransforms(transform.Identifier(2, 'test'), [transformseq.DimAxis(0,2,0,False),transformseq.DimAxis(0,2,0,False)], 0) + self.checktransforms = transformseq.IndexTransforms(2, 4) index, coords = self.mk_index_coords(2, self.checktransforms) self.basis = function.StructuredBasis([[[[1],[2]],[[3],[4]]],[[[5],[6]],[[7],[8]]]], [[0,1],[0,1]], [[2,3],[2,3]], [3,3], [2,2], index, coords) self.checkcoeffs = [[[[5.]],[[6.]],[[10.]],[[12.]]],[[[7.]],[[8.]],[[14.]],[[16.]]],[[[15.]],[[18.]],[[20.]],[[24.]]],[[[21.]],[[24.]],[[28.]],[[32.]]]] @@ -1244,19 +1244,22 @@ def setUp(self): topo, (x, y) = mesh.unitsquare(nelems=2, etype=self.etype) self.u = x * y * (1-y) self.manifold = topo.boundary['right'] + refgeom = None else: if self.etype == 'line': self.manifold, y = mesh.line(2) self.u = y * (2-y) + refgeom = numpy.stack([y]) else: self.manifold, (y, z) = mesh.unitsquare(nelems=2, etype=self.etype) self.u = y * (1-y) * z * (1-z) + refgeom = numpy.stack([y,z]) x = 1 # geometry describes a circle/sphere with curvature K self.geom = (x/self.K) * function.stack( (function.cos(y), function.sin(y)) if self.manifold.ndims == 1 else (function.cos(y), function.sin(y) * function.cos(z), function.sin(y) * function.sin(z))) - self.normal = function.normal(self.geom, exterior=not self.boundary) + self.normal = function.normal(self.geom, refgeom=refgeom) @property def P(self): diff --git a/tests/test_sample.py b/tests/test_sample.py index 3d421c01d..cc7b35ee0 100644 --- a/tests/test_sample.py +++ b/tests/test_sample.py @@ -3,7 +3,7 @@ from nutils.testing import * from nutils.sample import Sample from nutils.pointsseq import PointsSequence -from nutils.transformseq import IdentifierTransforms +from nutils.transformseq import IndexTransforms class Common: @@ -166,8 +166,8 @@ def setUp(self): super().setUp() points1 = [element.getsimplex(2).getpoints('bezier', 2)]*3 points2 = [element.getsimplex(2).getpoints('bezier', 3)]*2 - transforms1 = IdentifierTransforms(2, 'a', 3) - transforms2 = IdentifierTransforms(2, 'b', 2) + transforms1 = IndexTransforms(2, 3, 0) + transforms2 = IndexTransforms(2, 2, 3) sample1 = Sample.new('a', (transforms1, transforms1), PointsSequence.from_iter(points1, 2)) sample2 = Sample.new('a', (transforms2, transforms2), PointsSequence.from_iter(points2, 2)) self.sample = sample1 + sample2 @@ -200,8 +200,8 @@ def setUp(self): super().setUp() points1 = [element.getsimplex(1).getpoints('bezier', 2)]*2 points2 = [(element.getsimplex(1)**2).getpoints('bezier', 2)] + [element.getsimplex(2).getpoints('bezier', 2)]*2 - transforms1 = IdentifierTransforms(1, 'a', 2) - transforms2 = IdentifierTransforms(2, 'b', 3) + transforms1 = IndexTransforms(1, 2, 0) + transforms2 = IndexTransforms(2, 3, 2) sample1 = Sample.new('a', (transforms1, transforms1), PointsSequence.from_iter(points1, 1)) sample2 = Sample.new('b', (transforms2, transforms2), PointsSequence.from_iter(points2, 2)) self.sample = sample1 * sample2 @@ -218,8 +218,8 @@ def setUp(self): super().setUp() points1 = [element.getsimplex(0).getpoints('bezier', 2)] points2 = [(element.getsimplex(1)**2).getpoints('bezier', 2)] + [element.getsimplex(2).getpoints('bezier', 2)]*2 - transforms1 = IdentifierTransforms(0, 'a', 1) - transforms2 = IdentifierTransforms(2, 'b', 3) + transforms1 = IndexTransforms(0, 1, 0) + transforms2 = IndexTransforms(2, 3, 1) sample1 = Sample.new('a', (transforms1, transforms1), PointsSequence.from_iter(points1, 0)) sample2 = Sample.new('b', (transforms2, transforms2), PointsSequence.from_iter(points2, 2)) self.sample = sample1 * sample2 @@ -236,8 +236,8 @@ def setUp(self): super().setUp() points1 = [(element.getsimplex(1)**2).getpoints('bezier', 2)] + [element.getsimplex(2).getpoints('bezier', 2)]*2 points2 = [element.getsimplex(0).getpoints('bezier', 2)] - transforms1 = IdentifierTransforms(2, 'a', 3) - transforms2 = IdentifierTransforms(0, 'b', 1) + transforms1 = IndexTransforms(2, 3, 0) + transforms2 = IndexTransforms(0, 1, 3) sample1 = Sample.new('a', (transforms1, transforms1), PointsSequence.from_iter(points1, 2)) sample2 = Sample.new('b', (transforms2, transforms2), PointsSequence.from_iter(points2, 0)) self.sample = sample1 * sample2 @@ -295,9 +295,9 @@ def setUp(self): points1 = [element.getsimplex(1).getpoints('bezier', 2)]*2 points2 = [element.getsimplex(1).getpoints('bezier', 2), element.getsimplex(1).getpoints('bezier', 3), element.getsimplex(1).getpoints('bezier', 4)] points3 = [element.getsimplex(2).getpoints('bezier', 2)] + [(element.getsimplex(1)**2).getpoints('bezier', 2)] - transforms1 = IdentifierTransforms(1, 'a', 2) - transforms2 = IdentifierTransforms(1, 'b', 3) - transforms3 = IdentifierTransforms(2, 'c', 2) + transforms1 = IndexTransforms(1, 2, 0) + transforms2 = IndexTransforms(1, 3, 2) + transforms3 = IndexTransforms(2, 2, 5) sample1 = Sample.new('a', (transforms1, transforms1), PointsSequence.from_iter(points1, 1)) sample2 = Sample.new('b', (transforms2, transforms2), PointsSequence.from_iter(points2, 1)) sample3 = Sample.new('c', (transforms3, transforms3), PointsSequence.from_iter(points3, 2)) @@ -327,7 +327,7 @@ def setUp(self): line = element.getsimplex(1) triangle = element.getsimplex(2) points = [ref.getpoints('bezier', 2) for ref in (line**2, triangle, line**2)] - self.transforms = IdentifierTransforms(2, 'test', 3) + self.transforms = IndexTransforms(2, 3) self.sample = Sample.new('a', (self.transforms, self.transforms), PointsSequence.from_iter(points, 2)) self.desired_spaces = 'a', self.desired_ndims = 2 @@ -354,7 +354,7 @@ def setUp(self): line = element.getsimplex(1) triangle = element.getsimplex(2) points = [ref.getpoints('bezier', 2) for ref in (line**2, triangle, line**2)] - self.transforms = IdentifierTransforms(2, 'test', 3) + self.transforms = IndexTransforms(2, 3) self.desired_indices = [5,10,4,9],[2,0,6],[7,8,3,1] self.sample = Sample.new('a', (self.transforms, self.transforms), PointsSequence.from_iter(points, 2), tuple(map(numpy.array, self.desired_indices))) self.desired_spaces = 'a', diff --git a/tests/test_topology.py b/tests/test_topology.py index c6aee1163..cf27ed817 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -401,15 +401,15 @@ def test_opposite(self): self.assertEqual(len(~self.topo), 0) def test_intersection(self): - atrans = transformseq.IdentifierTransforms(1, 'a', 1) - btrans = transformseq.IdentifierTransforms(2, 'b', 1) + atrans = transformseq.IndexTransforms(1, 1, 0) + btrans = transformseq.IndexTransforms(2, 1, 1) other = topology.SimplexTopology('a', numpy.array([[0,1]]), atrans, atrans) * topology.SimplexTopology('b', numpy.array([[0,1,2]]), btrans, btrans) self.assertEqual(self.topo & other, self.topo) self.assertEqual(other & self.topo, self.topo) def test_union(self): - atrans = transformseq.IdentifierTransforms(1, 'a', 1) - btrans = transformseq.IdentifierTransforms(2, 'b', 1) + atrans = transformseq.IndexTransforms(1, 1, 0) + btrans = transformseq.IndexTransforms(2, 1, 1) other = topology.SimplexTopology('a', numpy.array([[0,1]]), atrans, atrans) * topology.SimplexTopology('b', numpy.array([[0,1,2]]), btrans, btrans) self.assertEqual(self.topo | other, other) self.assertEqual(other | self.topo, other) @@ -616,7 +616,7 @@ def assertConnectivity(self, domain, geom): imask = numpy.zeros(len(interfaces), dtype=int) coordinates = evaluable.Points(evaluable.NPoints(), boundary.ndims) transform_chain = transform.EvaluableTransformChain.from_argument('trans', domain.transforms.todims, boundary.ndims) - lowered_geom = geom.lower(coordinates.shape[:-1], {domain.space: (transform_chain,)*2}, {domain.space: coordinates}) + lowered_geom = geom.lower(coordinates.shape[:-1], {domain.space: (transform_chain,)*2}, {domain.space: coordinates}).simplified for ielem, ioppelems in enumerate(domain.connectivity): for iedge, ioppelem in enumerate(ioppelems): etrans, eref = domain.references[ielem].edges[iedge] @@ -794,7 +794,7 @@ class refined(TestCase): def test_boundary_gradient(self): ref = _refined_refs[self.etype] space = 'test' - domain = topology.ConnectedTopology(space, References.uniform(ref, 1), transformseq.IdentifierTransforms(ref.ndims, 'root', 1), transformseq.IdentifierTransforms(ref.ndims, 'root', 1), ((-1,)*ref.nedges,)).refine(self.ref0) + domain = topology.ConnectedTopology(space, References.uniform(ref, 1), transformseq.IndexTransforms(ref.ndims, 1), transformseq.IndexTransforms(ref.ndims, 1), ((-1,)*ref.nedges,)).refine(self.ref0) geom = function.rootcoords(space, ref.ndims) basis = domain.basis('std', degree=1) u = domain.projection(geom.sum(), onto=basis, geometry=geom, degree=2) @@ -1169,7 +1169,7 @@ class TransformChainsTopology(TestCase, CommonTests, TransformChainsTests): def setUp(self): super().setUp() references = References.uniform(element.PointReference(), 1) - transforms = transformseq.IdentifierTransforms(0, 'test', 1) + transforms = transformseq.IndexTransforms(0, 1) self.topo = topology.TransformChainsTopology('test', references, transforms, transforms) self.geom = self.topo.f_coords self.desired_nelems = 1 diff --git a/tests/test_transform.py b/tests/test_transform.py index 3ef0f7f5b..ee8b345c6 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -3,11 +3,11 @@ class specialcases(TestCase): - def test_tensoredge_swapup_identifier(self): + def test_tensoredge_swapup_index(self): lineedge = transform.SimplexEdge(1, 0, False) for edge in transform.TensorEdge1(lineedge, 1), transform.TensorEdge2(1, lineedge): with self.subTest(type(edge).__name__): - idnt = transform.Identifier(1, 'test') + idnt = transform.Index(1, 0) self.assertEqual(edge.swapup(idnt), None) class TestTransform(TestCase): @@ -60,20 +60,15 @@ class Qquare(TestInvertible): def setUp(self): super().setUp(trans=transform.Square([[1.,2],[1,3]], [5.,6]), linear=[[1,2],[1,3]], offset=[5,6]) -class Shift(TestInvertible): - - def setUp(self): - super().setUp(trans=transform.Shift([1.,2]), linear=[[1,0],[0,1]], offset=[1,2]) - class Identity(TestInvertible): def setUp(self): super().setUp(trans=transform.Identity(2), linear=[[1,0],[0,1]], offset=[0,0]) -class Scale(TestInvertible): +class Index(TestInvertible): def setUp(self): - super().setUp(trans=transform.Scale(2, offset=[1.,2]), linear=[[2,0],[0,2]], offset=[1,2]) + super().setUp(trans=transform.Index(2, 3), linear=[[1,0],[0,1]], offset=[0,0]) class SimplexEdge(TestUpdim): @@ -85,6 +80,11 @@ class SimplexChild(TestInvertible): def setUp(self): super().setUp(trans=transform.SimplexChild(3, 1), linear=numpy.eye(3)/2, offset=[.5,0,0]) +class Point(TestTransform): + + def setUp(self): + super().setUp(trans=transform.Point(numpy.array([1.,2.,3.])), linear=numpy.zeros((3,0)), offset=[1.,2.,3.]) + del TestTransform, TestInvertible, TestUpdim diff --git a/tests/test_transformseq.py b/tests/test_transformseq.py index 326496d28..21a0fe5ae 100644 --- a/tests/test_transformseq.py +++ b/tests/test_transformseq.py @@ -175,26 +175,22 @@ def test_edges(self): square = line*line triangle = nutils.element.TriangleReference() -l1, x1, r1 = sorted([nutils.transform.Identifier(1, n) for n in ('l1', 'x1', 'r1')], key=id) - -s0 = nutils.transform.Shift([0.]) -s1 = nutils.transform.Shift([1.]) -s2 = nutils.transform.Shift([2.]) -s3 = nutils.transform.Shift([3.]) -s4 = nutils.transform.Shift([4.]) +l1, x1, r1 = nutils.transform.Index(1, 100), nutils.transform.Index(1, 101), nutils.transform.Index(1, 102) + +i10 = nutils.transform.Index(1, 0) +i11 = nutils.transform.Index(1, 1) +i12 = nutils.transform.Index(1, 2) +i13 = nutils.transform.Index(1, 3) +i14 = nutils.transform.Index(1, 4) +i20 = nutils.transform.Index(2, 0) +i21 = nutils.transform.Index(2, 1) +i22 = nutils.transform.Index(2, 2) +i23 = nutils.transform.Index(2, 3) c0,c1 = line.child_transforms e0,e1 = line.edge_transforms -l2, x2, r2 = sorted([nutils.transform.Identifier(2, n) for n in ('l2', 'x2', 'r2')], key=id) -s00 = nutils.transform.Shift([0.,0.]) -s01 = nutils.transform.Shift([0.,1.]) -s02 = nutils.transform.Shift([0.,2.]) -s03 = nutils.transform.Shift([0.,3.]) -s10 = nutils.transform.Shift([1.,0.]) -s11 = nutils.transform.Shift([1.,1.]) -s12 = nutils.transform.Shift([1.,2.]) -s13 = nutils.transform.Shift([1.,3.]) +l2, x2, r2 = nutils.transform.Index(2, 100), nutils.transform.Index(2, 101), nutils.transform.Index(2, 102) c00,c01,c10,c11 = square.child_transforms e00,e01,e10,e11 = square.edge_transforms @@ -204,7 +200,7 @@ def setUp(self): super().setUp() self.seq = nutils.transformseq.EmptyTransforms(todims=1, fromdims=1) self.check = () - self.checkmissing = (l1,s0),(x1,s4),(r1,s0) + self.checkmissing = (l1,i10),(x1,i14),(r1,i10) self.checkrefs = References.empty(1) self.checktodims = 1 self.checkfromdims = 1 @@ -212,9 +208,9 @@ def setUp(self): class PlainTransforms1D(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1),(x1,s2),(x1,s3)], 1, 1) - self.check = (x1,s0),(x1,s1),(x1,s2),(x1,s3) - self.checkmissing = (l1,s0),(x1,s4),(r1,s0) + self.seq = nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11),(x1,i12),(x1,i13)], 1, 1) + self.check = (x1,i10),(x1,i11),(x1,i12),(x1,i13) + self.checkmissing = (l1,i10),(x1,i14),(r1,i10) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -222,9 +218,9 @@ def setUp(self): class PlainTransforms2D(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.PlainTransforms([(x2,s00),(x2,s01),(x2,s10),(x2,s11)], 2, 2) - self.check = (x2,s00),(x2,s01),(x2,s10),(x2,s11) - self.checkmissing = (l2,s00),(x2,s02),(x2,s12),(r2,s00) + self.seq = nutils.transformseq.PlainTransforms([(x2,i20,i20),(x2,i20,i21),(x2,i21,i20),(x2,i21,i21)], 2, 2) + self.check = (x2,i20,i20),(x2,i20,i21),(x2,i21,i20),(x2,i21,i21) + self.checkmissing = (l2,i20,i20),(x2,i20,i22),(x2,i21,i22),(r2,i20,i20) self.checkrefs = References.from_iter((square,square,triangle,triangle), 2) self.checktodims = 2 self.checkfromdims = 2 @@ -232,9 +228,9 @@ def setUp(self): class MaskedTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.MaskedTransforms(nutils.transformseq.PlainTransforms([(x2,s00),(x2,s01),(x2,s10),(x2,s11)], 2, 2), [0,2]) - self.check = (x2,s00),(x2,s10) - self.checkmissing = (l2,s00),(x2,s01),(x2,s11),(x2,s02),(x2,s12),(r2,s00) + self.seq = nutils.transformseq.MaskedTransforms(nutils.transformseq.PlainTransforms([(x2,i20,i20),(x2,i20,i21),(x2,i21,i20),(x2,i21,i21)], 2, 2), [0,2]) + self.check = (x2,i20,i20),(x2,i21,i20) + self.checkmissing = (l2,i20,i20),(x2,i20,i21),(x2,i21,i21),(x2,i20,i22),(x2,i21,i22),(r2,i20,i20) self.checkrefs = References.from_iter((square,triangle), 2) self.checktodims = 2 self.checkfromdims = 2 @@ -242,9 +238,9 @@ def setUp(self): class ReorderedTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.ReorderedTransforms(nutils.transformseq.PlainTransforms([(x2,s00),(x2,s01),(x2,s10),(x2,s11)], 2, 2), [0,2,3,1]) - self.check = (x2,s00),(x2,s10),(x2,s11),(x2,s01) - self.checkmissing = (l2,s00),(x2,s02),(x2,s12),(r2,s00) + self.seq = nutils.transformseq.ReorderedTransforms(nutils.transformseq.PlainTransforms([(x2,i20,i20),(x2,i20,i21),(x2,i21,i20),(x2,i21,i21)], 2, 2), [0,2,3,1]) + self.check = (x2,i20,i20),(x2,i21,i20),(x2,i21,i21),(x2,i20,i21) + self.checkmissing = (l2,i20,i20),(x2,i20,i22),(x2,i21,i22),(r2,i20,i20) self.checkrefs = References.uniform(square, 4) self.checktodims = 2 self.checkfromdims = 2 @@ -252,9 +248,9 @@ def setUp(self): class DerivedTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.DerivedTransforms(nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1), References.uniform(line, 2), 'child_transforms', 1) - self.check = (x1,s0,c0),(x1,s0,c1),(x1,s1,c0),(x1,s1,c1) - self.checkmissing = (l1,s0),(x1,s0),(x1,s1),(r1,s0) + self.seq = nutils.transformseq.DerivedTransforms(nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1), References.uniform(line, 2), 'child_transforms', 1) + self.check = (x1,i10,c0),(x1,i10,c1),(x1,i11,c0),(x1,i11,c1) + self.checkmissing = (l1,i10),(x1,i10),(x1,i11),(r1,i10) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -262,9 +258,9 @@ def setUp(self): class UniformDerivedTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.UniformDerivedTransforms(nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1), line, 'child_transforms', 1) - self.check = (x1,s0,c0),(x1,s0,c1),(x1,s1,c0),(x1,s1,c1) - self.checkmissing = (l1,s0),(x1,s0),(x1,s1),(r1,s0) + self.seq = nutils.transformseq.UniformDerivedTransforms(nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1), line, 'child_transforms', 1) + self.check = (x1,i10,c0),(x1,i10,c1),(x1,i11,c0),(x1,i11,c1) + self.checkmissing = (l1,i10),(x1,i10),(x1,i11),(r1,i10) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -272,9 +268,9 @@ def setUp(self): class ChainedTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.ChainedTransforms([nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1), nutils.transformseq.PlainTransforms([(x1,s2),(x1,s3)], 1, 1)]) - self.check = (x1,s0),(x1,s1),(x1,s2),(x1,s3) - self.checkmissing = (l1,s0),(x1,s4),(r1,s0) + self.seq = nutils.transformseq.ChainedTransforms([nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1), nutils.transformseq.PlainTransforms([(x1,i12),(x1,i13)], 1, 1)]) + self.check = (x1,i10),(x1,i11),(x1,i12),(x1,i13) + self.checkmissing = (l1,i10),(x1,i14),(r1,i10) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -283,8 +279,8 @@ class StructuredTransforms1D(TestCase, Common, Edges): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.DimAxis(0,4,0,False)], 0) - self.check = (x1,s0),(x1,s1),(x1,s2),(x1,s3) - self.checkmissing = (l1,s0),(x1,s4),(r1,s0),(x1,c1) + self.check = (x1,i10),(x1,i11),(x1,i12),(x1,i13) + self.checkmissing = (l1,i10),(x1,),(x1,i14),(r1,i10),(x1,c1) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -293,8 +289,8 @@ class StructuredTransforms1DRefined(TestCase, Common, Edges): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.DimAxis(0,4,0,False)], 1) - self.check = (x1,s0,c0),(x1,s0,c1),(x1,s1,c0),(x1,s1,c1) - self.checkmissing = (l1,s0),(x1,s0),(x1,s1),(x1,s0,s1),(r1,s0) + self.check = (x1,i10,c0),(x1,i10,c1),(x1,i11,c0),(x1,i11,c1) + self.checkmissing = (l1,i10),(x1,),(x1,i10),(x1,i11),(x1,i10,i11),(r1,i10) self.checkrefs = References.uniform(line, 4) self.checktodims = 1 self.checkfromdims = 1 @@ -303,8 +299,8 @@ class StructuredTransforms1DLeft(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(3,4,9,0,False)], 0) - self.check = (x1,s3,e1), - self.checkmissing = (x1,s0,e0),(x1,s2,e0),(x1,s4,e0) + self.check = (x1,i13,e1), + self.checkmissing = (x1,i10,e0),(x1,i12,e0),(x1,i14,e0) self.checkrefs = References.uniform(point, 1) self.checktodims = 1 self.checkfromdims = 0 @@ -313,8 +309,8 @@ class StructuredTransforms1DRight(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(2,3,9,0,True)], 0) - self.check = (x1,s2,e0), - self.checkmissing = (x1,s0,e0),(x1,s3,e1),(x1,s4,e0) + self.check = (x1,i12,e0), + self.checkmissing = (x1,i10,e0),(x1,i13,e1),(x1,i14,e0) self.checkrefs = References.uniform(point, 1) self.checktodims = 1 self.checkfromdims = 0 @@ -323,8 +319,8 @@ class StructuredTransforms1DInterfacesLeft(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(1,4,9,0,False)], 0) - self.check = (x1,s1,e1),(x1,s2,e1),(x1,s3,e1) - self.checkmissing = (x1,s0,e1),(x1,s0,e0),(x1,s1,e0),(x1,s2,e0),(x1,s3,e0) + self.check = (x1,i11,e1),(x1,i12,e1),(x1,i13,e1) + self.checkmissing = (x1,i10,e1),(x1,i10,e0),(x1,i11,e0),(x1,i12,e0),(x1,i13,e0) self.checkrefs = References.uniform(point, 3) self.checktodims = 1 self.checkfromdims = 0 @@ -333,8 +329,8 @@ class StructuredTransforms1DInterfacesRight(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(0,3,9,0,True)], 0) - self.check = (x1,s0,e0),(x1,s1,e0),(x1,s2,e0) - self.checkmissing = (x1,s3,e0),(x1,s0,e1),(x1,s1,e1),(x1,s2,e1),(x1,s3,e1) + self.check = (x1,i10,e0),(x1,i11,e0),(x1,i12,e0) + self.checkmissing = (x1,i13,e0),(x1,i10,e1),(x1,i11,e1),(x1,i12,e1),(x1,i13,e1) self.checkrefs = References.uniform(point, 3) self.checktodims = 1 self.checkfromdims = 0 @@ -343,8 +339,8 @@ class StructuredTransforms1DPeriodicInterfacesLeft(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(1,5,4,0,False)], 0) - self.check = (x1,s1,e1),(x1,s2,e1),(x1,s3,e1),(x1,s0,e1) - self.checkmissing = (x1,s0,e0),(x1,s1,e0),(x1,s2,e0),(x1,s3,e0),(x1,s4,e0) + self.check = (x1,i11,e1),(x1,i12,e1),(x1,i13,e1),(x1,i10,e1) + self.checkmissing = (x1,i10,e0),(x1,i11,e0),(x1,i12,e0),(x1,i13,e0),(x1,i14,e0) self.checkrefs = References.uniform(point, 4) self.checktodims = 1 self.checkfromdims = 0 @@ -353,8 +349,8 @@ class StructuredTransforms1DPeriodicInterfacesRight(TestCase, Common): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x1, [nutils.transformseq.IntAxis(0,4,4,0,True)], 0) - self.check = (x1,s0,e0),(x1,s1,e0),(x1,s2,e0),(x1,s3,e0) - self.checkmissing = (x1,s0,e1),(x1,s1,e1),(x1,s2,e1),(x1,s3,e1),(x1,s4,e1) + self.check = (x1,i10,e0),(x1,i11,e0),(x1,i12,e0),(x1,i13,e0) + self.checkmissing = (x1,i10,e1),(x1,i11,e1),(x1,i12,e1),(x1,i13,e1),(x1,i14,e1) self.checkrefs = References.uniform(point, 4) self.checktodims = 1 self.checkfromdims = 0 @@ -363,8 +359,8 @@ class StructuredTransforms2D(TestCase, Common, Edges): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x2, [nutils.transformseq.DimAxis(0,2,0,False),nutils.transformseq.DimAxis(2,4,0,False)], 0) - self.check = (x2,s02),(x2,s03),(x2,s12),(x2,s13) - self.checkmissing = (x2,s00),(x2,s01),(x2,s10),(x2,s11) + self.check = (x2,i20,i22),(x2,i20,i23),(x2,i21,i22),(x2,i21,i23) + self.checkmissing = (x2,),(x2,i20),(x2,i20,i20),(x2,i20,i21),(x2,i21,i20),(x2,i21,i21) self.checkrefs = References.uniform(square, 4) self.checktodims = 2 self.checkfromdims = 2 @@ -373,18 +369,18 @@ class StructuredTransforms2DRefined(TestCase, Common, Edges): def setUp(self): super().setUp() self.seq = nutils.transformseq.StructuredTransforms(x2, [nutils.transformseq.DimAxis(0,2,0,False),nutils.transformseq.DimAxis(2,4,0,False)], 1) - self.check = (x2,s01,c00),(x2,s01,c01),(x2,s01,c10),(x2,s01,c11) - self.checkmissing = (x2,s00,c00), + self.check = (x2,i20,i21,c00),(x2,i20,i21,c01),(x2,i20,i21,c10),(x2,i20,i21,c11) + self.checkmissing = (x2,i20,i20,c00), self.checkrefs = References.uniform(square, 4) self.checktodims = 2 self.checkfromdims = 2 -class IdentifierTransforms(TestCase, Common, Edges): +class IndexTransforms(TestCase, Common, Edges): def setUp(self): super().setUp() - self.seq = nutils.transformseq.IdentifierTransforms(ndims=2, name='foo', length=4) - self.check = [(nutils.transform.Identifier(2, ('foo', i)),) for i in range(4)] - self.checkmissing = (nutils.transform.Identifier(1, ('foo', 0)),), (nutils.transform.Identifier(2, ('foo', -1)),), (nutils.transform.Identifier(2, ('foo', 4)),), (nutils.transform.Identifier(2, ('bar', 0)),) + self.seq = nutils.transformseq.IndexTransforms(ndims=2, length=4, offset=2) + self.check = [(nutils.transform.Index(2, i),) for i in range(2, 6)] + self.checkmissing = (nutils.transform.Index(1, 0),), (nutils.transform.Index(2, -1),), (nutils.transform.Index(2, 6,),) self.checkrefs = References.uniform(triangle, 4) self.checktodims = 2 self.checkfromdims = 2 @@ -407,22 +403,22 @@ def test_PlainTransforms_invalid_fromdims(self): def test_PlainTransforms_multiple_fromdims(self): with self.assertRaisesRegex(ValueError, 'expected transforms with fromdims=2, but got .*'): - nutils.transformseq.PlainTransforms([(x2,c00,e00),(x2,c01,s00)], 2, 2) + nutils.transformseq.PlainTransforms([(x2,c00,e00),(x2,c01,i20,i20)], 2, 2) def test_DerivedTransforms_length_mismatch(self): - transforms = nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1) + transforms = nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1) references = References.uniform(line, 3) with self.assertRaisesRegex(ValueError, '`parent` and `parent_references` should have the same length'): nutils.transformseq.DerivedTransforms(transforms, references, 'child_transforms', 1) def test_DerivedTransforms_ndims_mismatch(self): - transforms = nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1) + transforms = nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1) references = References.uniform(square, 2) with self.assertRaisesRegex(ValueError, '`parent` and `parent_references` have different dimensions'): nutils.transformseq.DerivedTransforms(transforms, references, 'child_transforms', 1) def test_UniformDerivedTransforms_ndims_mismatch(self): - transforms = nutils.transformseq.PlainTransforms([(x1,s0),(x1,s1)], 1, 1) + transforms = nutils.transformseq.PlainTransforms([(x1,i10),(x1,i11)], 1, 1) with self.assertRaisesRegex(ValueError, '`parent` and `parent_reference` have different dimensions'): nutils.transformseq.UniformDerivedTransforms(transforms, square, 'child_transforms', 1)