From 2cda5201c436d9f4060150da33d3477a9adc2636 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Sat, 12 Dec 2020 21:07:10 +0100 Subject: [PATCH 1/3] support ndim!=2 coordinates in lower --- nutils/evaluable.py | 5 +++-- nutils/function.py | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index e43f519b0..f5be14571 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -1446,7 +1446,7 @@ class ApplyTransforms(Array): def __init__(self, trans:types.strict[TransformChain], points, todims:types.strictint): self.trans = trans self._todims = todims - super().__init__(args=[points, trans], shape=[points.shape[0], todims], dtype=float) + super().__init__(args=[points, trans], shape=points.shape[:-1]+(todims,), dtype=float) def evalf(self, points, chain): return transform.apply(chain, points) @@ -2300,6 +2300,7 @@ class Sampled(Array): @types.apply_annotations def __init__(self, points:asarray, expect:asarray): + assert points.ndim == 2 super().__init__(args=[points, expect], shape=(points.shape[0], expect.shape[0]), dtype=int) def evalf(self, points, expect): @@ -2840,7 +2841,7 @@ def evalf(self, evalargs): raise ValueError('argument {!r} missing'.format(self._name)) else: value = numpy.asarray(value) - assert value.shape == self.shape + assert value.shape == self.shape, 'invalid argument shape: got {}, expected {}'.format(value.shape, self.shape) value = value.astype(self.dtype, casting='safe', copy=False) return value diff --git a/nutils/function.py b/nutils/function.py index 432b494c9..41f938a33 100644 --- a/nutils/function.py +++ b/nutils/function.py @@ -50,7 +50,7 @@ def lower(self, *, transform_chains: Tuple[evaluable.TransformChain] = (), coord def _lower(self, **kwargs): result = self._ArrayMeta__lower(**kwargs) assert isinstance(result, evaluable.Array) - offset = 1 if kwargs.get('coordinates', ()) and not type(self).__name__ == '_WithoutPoints' else 0 + offset = kwargs['coordinates'][0].ndim-1 if kwargs.get('coordinates', ()) else 0 assert result.ndim == self.ndim + offset for n, m in zip(result.shape[offset:], self.shape): if isinstance(m, int): @@ -578,7 +578,8 @@ def __getnewargs__(self): return self._arg, self._axes def lower(self, **kwargs: Any) -> evaluable.Array: - axes = (0, *(i+1 for i in self._axes)) if kwargs.get('coordinates', ()) else self._axes + offset = kwargs['coordinates'][0].ndim-1 if kwargs.get('coordinates', ()) else 0 + axes = (*range(offset), *(i+offset for i in self._axes)) return evaluable.Transpose(self._arg.lower(**kwargs), axes) class _Opposite(Array): From be904b8dbd1a2c8bd4b68ddd3462bf2d6055fe84 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 6 Jan 2021 14:03:27 +0100 Subject: [PATCH 2/3] repurpose locate eps argument as newton tolerance --- nutils/topology.py | 95 +++++++++++++++++++++++++++------------------- 1 file changed, 57 insertions(+), 38 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index e64c6d8a3..c2629e46a 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -448,7 +448,7 @@ def select(self, indicator, ischeme='bezier2', **kwargs): return self[selected] @log.withcontext - def locate(self, geom, coords, *, tol, eps=0, maxiter=100, 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): '''Create a sample based on physical coordinates. In a finite element application, functions are commonly evaluated in points @@ -467,9 +467,10 @@ def locate(self, geom, coords, *, tol, eps=0, maxiter=100, arguments=None, weigh >>> sample.eval(geom).tolist() [[0.9, 0.4]] - Locate has a long list of arguments that can be used to steer the nonlinear - search process, but the default values should be fine for reasonably - standard situations. + Locate requires a geometry function, an array of coordinates, and at least + one of ``tol`` and ``eps`` to set the tolerance in physical of element + space, respectively; if both are specified the least restrictive takes + precedence. Args ---- @@ -477,21 +478,24 @@ def locate(self, geom, coords, *, tol, eps=0, maxiter=100, arguments=None, weigh Geometry function of length ``ndims``. coords : 2-dimensional :class:`float` array Array of coordinates with ``ndims`` columns. - tol : :class:`float` - Maximum allowed distance between original and located coordinate. + tol : :class:`float` (default: 0) + Maximum allowed distance in physical coordinates between target and + located point. eps : :class:`float` (default: 0) - Epsilon radius around element within which a point is considered to be - inside. - maxiter : :class:`int` (default: 100) - Maximum allowed number of Newton iterations. + Maximum allowed distance in element coordinates between target and + located point. + maxiter : :class:`int` (default: 0) + Maximum allowed number of Newton iterations, or 0 for unlimited. arguments : :class:`dict` (default: None) Arguments for function evaluation. weights : :class:`float` array (default: None) - Optional weights, in case ``coords`` are quadrature points. + Optional weights, in case ``coords`` are quadrature points, making the + resulting sample suitable for integration. maxdist : :class:`float` (default: None) - Speed up failure by setting a distance between point and element - centroid above which the element is rejected immediately. If all points - are expected to be located then this can safely be left unspecified. + Speed up failure by setting a physical distance between point and + element centroid above which the element is rejected immediately. If + all points are expected to be located then this can safely be left + unspecified. Returns ------- @@ -502,47 +506,62 @@ def locate(self, geom, coords, *, tol, eps=0, maxiter=100, arguments=None, weigh warnings.deprecation('the ischeme argument is deprecated and will be removed in future') if scale is not None: warnings.deprecation('the scale argument is deprecated and will be removed in future') + if max(tol, eps) <= 0: + raise ValueError('locate requires either tol or eps to be strictly positive') coords = numpy.asarray(coords, dtype=float) if geom.ndim == 0: geom = geom[_] coords = coords[...,_] if not geom.shape == coords.shape[1:] == (self.ndims,): - raise Exception('invalid geometry or point shape for {}D topology'.format(self.ndims)) + raise ValueError('invalid geometry or point shape for {}D topology'.format(self.ndims)) centroids = self.sample('_centroid', None).eval(geom) assert len(centroids) == len(self) ielems = parallel.shempty(len(coords), dtype=int) - xis = parallel.shempty((len(coords),len(geom)), dtype=float) - J = function.localgradient(geom, self.ndims) - geom_J = evaluable.Tuple((geom.prepare_eval(ndims=self.ndims), J.prepare_eval(ndims=self.ndims))).simplified + points = parallel.shempty((len(coords),len(geom)), dtype=float) + _ielem = evaluable.Argument('_locate_ielem', shape=(), dtype=int) + _point = evaluable.Argument('_locate_point', shape=(self.ndims,)) + lower_args = dict( + transform_chains = ( + evaluable.TransformChainFromSequence(self.transforms, _ielem), + evaluable.TransformChainFromSequence(self.opposites, _ielem)), + coordinates = (_point, _point)) + xJ = evaluable.Tuple((geom.lower(**lower_args), function.localgradient(geom, self.ndims).lower(**lower_args))).simplified + arguments = dict(arguments or ()) with parallel.ctxrange('locating', len(coords)) as ipoints: for ipoint in ipoints: - coord = coords[ipoint] - dist = numpy.linalg.norm(centroids - coord, axis=1) + xt = coords[ipoint] # target + dist = numpy.linalg.norm(centroids - xt, axis=1) for ielem in numpy.argsort(dist) if maxdist is None \ else sorted((dist < maxdist).nonzero()[0], key=dist.__getitem__): - converged = False ref = self.references[ielem] - xi = numpy.array(ref.centroid) - for iiter in range(maxiter): - (coord_xi,), (J_xi,) = geom_J.eval(_transforms=(self.transforms[ielem], self.opposites[ielem]), _points=points.CoordsPoints(xi[_]), **arguments or {}) - err = numpy.linalg.norm(coord - coord_xi) - if err < tol: - converged = True - break - if iiter and err > prev_err: - break - prev_err = err + arguments['_locate_ielem'] = ielem + arguments['_locate_point'] = p = numpy.array(ref.centroid) + ex = ep = numpy.inf + iiter = 0 + while ex > tol and ep > eps: # newton loop + if iiter > maxiter > 0: + break # maximum number of iterations reached + iiter += 1 + xp, Jp = xJ.eval(**arguments) + dx = xt - xp + ex0 = ex + ex = numpy.linalg.norm(dx) + if ex >= ex0: + break # newton is diverging try: - xi += numpy.linalg.solve(J_xi, coord - coord_xi) + dp = numpy.linalg.solve(Jp, dx) except numpy.linalg.LinAlgError: + break # jacobian is singular + ep = numpy.linalg.norm(dp) + p += dp # NOTE: modifies arguments['_locate_point'] in place + else: + if ref.inside(p, max(eps, ep)): + ielems[ipoint] = ielem + points[ipoint] = p break - if converged and ref.inside(xi, eps=max(eps, *tol/abs(numpy.linalg.eigvals(J_xi)))): - ielems[ipoint] = ielem - xis[ipoint] = xi - break else: - raise LocateError('failed to locate point: {}'.format(coord)) - return self._sample(ielems, xis, weights) + raise LocateError('failed to locate point: {}'.format(xt)) + return self._sample(ielems, points, weights) def _sample(self, ielems, coords, weights=None): uielems = numpy.unique(ielems) From 4e140ae3acedea804ac99c95ed2777dda057dfe1 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 6 Jan 2021 14:30:57 +0100 Subject: [PATCH 3/3] add changelog entry --- CHANGELOG.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index de479a817..8edcbdb11 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -9,6 +9,17 @@ features in inverse chronological order. New in v7.0 (in development) ---------------------------- +- Changed: locate arguments + + The :func:`nutils.topology.Topology.locate` method now allows ``tol`` to be + left unspecified if ``eps`` is specified instead, which is repurposed as stop + criterion for distances in element coordinates. Conversely, if only ``tol`` + is specified, a corresponding minimal ``eps`` value is set automatically to + match points near element edges. The ``ischeme`` and ``scale`` arguments are + deprecated and replaced by ``maxdist``, which can be left unspecified in + general. The optional ``weights`` argument results in a sample that is + suitable for integration. + - Moved: unit from types to separate module The ``unit`` type has been moved into its own :mod:`nutils.unit` module, with