Merge pull request #610 from evalf/locate

Locate

nutils/element.py

@@ -155,6 +155,12 @@ def getischeme(self, ischeme):
     return points.coords, getattr(points, 'weights', None)
 
   def getpoints(self, ischeme, degree):
+    if ischeme == '_centroid':
+      gauss = self.getpoints('gauss', 1)
+      if gauss.npoints == 1:
+        return gauss
+      volume = gauss.weights.sum()
+      return points.CoordsUniformPoints(gauss.coords.T[_] @ gauss.weights / volume, volume)
     raise Exception('unsupported ischeme for {}: {!r}'.format(self.__class__.__name__, ischeme))
 
   def with_children(self, child_refs):
@@ -167,12 +173,13 @@ def with_children(self, child_refs):
 
   @property
   def volume(self):
-    return self.getpoints('gauss', 1).weights.sum()
+    volume, = self.getpoints('_centroid', None).weights
+    return volume
 
   @property
   def centroid(self):
-    gauss = self.getpoints('gauss', 1)
-    return types.frozenarray(gauss.coords.T.dot(gauss.weights) / gauss.weights.sum(), copy=False)
+    centroid, = self.getpoints('_centroid', None).coords
+    return centroid
 
   def trim(self, levels, maxrefine, ndivisions):
     'trim element along levelset'
@@ -774,7 +781,7 @@ class Cone(Reference):
   'cone'
 
   __slots__ = 'edgeref', 'etrans', 'tip', 'extnorm', 'height'
-  __cache__ = 'vertices', 'edge_transforms', 'edge_refs', 'volume'
+  __cache__ = 'vertices', 'edge_transforms', 'edge_refs'
 
   @types.apply_annotations
   def __init__(self, edgeref, etrans, tip:types.frozenarray):
@@ -827,17 +834,16 @@ def getpoints(self, ischeme, degree):
       tx, tw = points.gauss((degree + self.ndims - 1)//2)
       wx = tx**(self.ndims-1) * tw * self.extnorm * self.height
       return points.CoordsWeightsPoints((tx[:,_,_] * (self.etrans.apply(epoints.coords)-self.tip)[_,:,:] + self.tip).reshape(-1, self.ndims), (epoints.weights[_,:] * wx[:,_]).ravel())
-    if ischeme == 'uniform':
-      coords = numpy.concatenate([(self.etrans.apply(self.edgeref.getpoints('uniform', i+1).coords) - self.tip) * ((i+.5)/degree) + self.tip for i in range(degree)])
-      return points.CoordsUniformPoints(coords, self.volume)
+    if ischeme in ('_centroid', 'uniform'):
+      layers = [(i+1, (i+.5)/degree) for i in range(degree)] if ischeme == 'uniform' \
+          else [(None, self.ndims/(self.ndims+1))] # centroid
+      coords = numpy.concatenate([self.etrans.apply(self.edgeref.getpoints(ischeme, p).coords) * w + self.tip * (1-w) for p, w in layers])
+      volume = self.edgeref.volume * self.extnorm * self.height / self.ndims
+      return points.CoordsUniformPoints(coords, volume)
     if ischeme == 'vtk' and self.nverts == 5 and self.ndims==3: # pyramid
       return points.CoordsPoints(self.vertices[[1,2,4,3,0]])
     return points.ConePoints(self.edgeref.getpoints(ischeme, degree), self.etrans, self.tip)
 
-  @property
-  def volume(self):
-    return self.edgeref.volume * self.extnorm * self.height / self.ndims
-
   @property
   def simplices(self):
     if self.nverts == self.ndims+1 or self.edgeref.ndims == 2 and self.edgeref.nverts == 4: # simplices and square-based pyramids are ok
@@ -964,10 +970,15 @@ def subvertex(self, ichild, i):
   def getpoints(self, ischeme, degree):
     if ischeme == 'vertex':
       return self.baseref.getpoints(ischeme, degree)
+    if ischeme == '_centroid':
+      return super().getpoints(ischeme, degree)
     if ischeme == 'bezier':
-      childpoints = [points.TransformPoints(ref.getpoints('bezier', degree//2+1), trans) for trans, ref in self.children if ref]
-      return points.ConcatPoints(childpoints, points.find_duplicates(childpoints))
-    return points.ConcatPoints(points.TransformPoints(ref.getpoints(ischeme, degree), trans) for trans, ref in self.children if ref)
+      degree = degree//2+1 # modify child degree to keep (approximate) uniformity
+      dedup = True
+    else:
+      dedup = False
+    childpoints = [points.TransformPoints(ref.getpoints(ischeme, degree), trans) for trans, ref in self.children if ref]
+    return points.ConcatPoints(childpoints, points.find_duplicates(childpoints) if dedup else ())
 
   @property
   def simplices(self):
@@ -1120,8 +1131,14 @@ def simplices(self):
   def getpoints(self, ischeme, degree):
     if ischeme == 'vertex':
       return self.baseref.getpoints(ischeme, degree)
+    if ischeme == '_centroid':
+      return super().getpoints(ischeme, degree)
     subpoints = [subvol.getpoints(ischeme, degree) for subvol in self.subrefs]
     dups = points.find_duplicates(subpoints) if ischeme == 'bezier' else ()
+    # NOTE We could consider postprocessing gauss1 to a single point scheme,
+    # rather than a concatenation, but for that we would have to verify that
+    # the centroid is contained in the element. We leave this optimization for
+    # later, to be combined with a reduction of gauss schemes of any degree.
     return points.ConcatPoints(subpoints, dups)
 
   def inside(self, point, eps=0):

nutils/mesh.py

@@ -691,7 +691,7 @@ def unitsquare(nelems, etype):
         connectivity = [c-numpy.greater(c,n*2) for c in connectivity]
       topo = topology.ConnectedTopology(References.from_iter(references, 2), transformseq.PlainTransforms(transforms, 2),transformseq.PlainTransforms(transforms, 2), tuple(types.frozenarray(c, copy=False) for c in connectivity))
 
-    x, y = topo.boundary.elem_mean(function.rootcoords(2), degree=1).T
+    x, y = topo.boundary.sample('_centroid', None).eval(function.rootcoords(2)).T
     bgroups = dict(left=x==0, right=x==nelems, bottom=y==0, top=y==nelems)
     topo = topo.withboundary(**{name: topo.boundary[numpy.where(mask)[0]] for name, mask in bgroups.items()})
 

nutils/topology.py

@@ -508,7 +508,8 @@ def locate(self, geom, coords, *, tol, eps=0, maxiter=100, arguments=None, weigh
       coords = coords[...,_]
     if not geom.shape == coords.shape[1:] == (self.ndims,):
       raise Exception('invalid geometry or point shape for {}D topology'.format(self.ndims))
-    centroids = self.elem_mean(geom, geometry=geom, degree=2)
+    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)
@@ -521,23 +522,23 @@ def locate(self, geom, coords, *, tol, eps=0, maxiter=100, arguments=None, weigh
                 else sorted((dist < maxdist).nonzero()[0], key=dist.__getitem__):
           converged = False
           ref = self.references[ielem]
-          p = ref.getpoints('gauss', 1)
-          xi = p.coords
-          w = p.weights
-          xi = (numpy.dot(w,xi) / w.sum())[_] if len(xi) > 1 else xi.copy()
+          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 {})
+            (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
-            xi += numpy.linalg.solve(J_xi, coord - coord_xi)
-          if converged and ref.inside(xi[0], eps=eps):
+            try:
+              xi += numpy.linalg.solve(J_xi, coord - coord_xi)
+            except numpy.linalg.LinAlgError:
+              break
+          if converged and ref.inside(xi, eps=max(eps, *tol/abs(numpy.linalg.eigvals(J_xi)))):
             ielems[ipoint] = ielem
-            xis[ipoint], = xi
+            xis[ipoint] = xi
             break
         else:
           raise LocateError('failed to locate point: {}'.format(coord))

nutils/types.py

@@ -1271,6 +1271,7 @@ def __index__(self):
   __rsub__ = lambda self, other: self.__base.__rsub__(other)
   __mul__ = lambda self, other: self.__base.__mul__(other)
   __rmul__ = lambda self, other: self.__base.__rmul__(other)
+  __matmul__ = lambda self, other: self.__base.__matmul__(other)
   __truediv__ = lambda self, other: self.__base.__truediv__(other)
   __rtruediv__ = lambda self, other: self.__base.__rtruediv__(other)
   __floordiv__ = lambda self, other: self.__base.__floordiv__(other)