updated the interp.grid object to be a better base class

interp/grid/__init__.py

@@ -7,31 +7,123 @@ import numpy as np
 import scipy.spatial
 
 from interp.baker import run_baker
-from interp.tools import exact_func, smberror, log
+from interp.tools import exact_func, log
 from simplex      import face, contains
 from smcqdelaunay import *
 
 
 class grid(object):
-  def __init__(self, points, q):
+  def __init__(self, verts = None, q = None):
     """
       this thing eats two pre-constructed arrays of floats:
-      points = array of arrays (i will convert to numpy.array)
+
-               [[x0,y0], [x1,y1], ...]
+      verts = array of arrays (if passed in, will convert to numpy.array)
-      q      = array (1D) of important values
+              [
+                [x0,y0],
+                [x1,y1], ...
+              ]
+
+      q     = array (1D) of physical values
     """
 
-    self.points           = np.array(points)
+    if verts:
+      self.verts = np.array(verts)
+    if q:
       self.q = np.array(q)
-    self.tree             = scipy.spatial.KDTree(self.points)
+
+    self.tree           = scipy.spatial.KDTree(self.verts)
     self.faces          = {}
     self.faces_for_vert = defaultdict(list)
 
+  def get_containing_simplex(self, X):
+    if not self.faces:
+      self.construct_connectivity()
+
+    # get closest point
+    (dist, indicies) = self.tree.query(X, 2)
+    closest_point = indicies[0]
+
+    log.debug('X: %s' % X)
+    log.debug('point index: %d' % closest_point)
+    log.debug('actual point %s' % self.verts[closest_point])
+    log.debug('distance = %0.4f' % dist[0])
+
+    simplex = None
+    checked_faces = []
+    faces_to_check = self.faces_for_vert[closest_point]
+
+    attempts = 0
+    while not simplex and faces_to_check:
+      attempts += 1
+
+      # if attempts > 20:
+      #   raise Exception("probably recursing to many times")
+
+      cur_face = faces_to_check.pop(0)
+      checked_faces.append(cur_face)
+
+      if cur_face.contains(X, self):
+        simplex = cur_face
+        continue
+
+      new_facest = []
+      for neighbor in cur_face.neighbors:
+        if (neighbor not in checked_faces) and (neighbor not in faces_to_check):
+          faces_to_check.append(neighbor)
+
+    if not simplex:
+      raise AssertionError('no containing simplex found')
+
+    R = self.create_mesh(simplex.verts)
+
+    log.debug('total attempts before finding simplex: %d' % attempts)
+    return R
+
+  def create_mesh(self, indicies):
+    """
+      this function takes a list of indicies, and then creates
+      and returns a grid object (collection of verts and q).
+
+      note: the input is indicies, the grid contains verts
+    """
+    p = [self.verts[i] for i in indicies]
+    q = [self.q[i]      for i in indicies]
+    return grid(p, q)
+
+  def get_simplex_and_nearest_points(self, X, extra_points = 3, simplex_size = 3):
+    """
+      this returns two grid objects: R and S.
+
+      R is a grid object that is supposedly a containing simplex
+        around point X (it tends not to be)
+
+      S is S_j from baker's paper : some verts from all point that are not the simplex
+    """
+    raise Exception("abstract base class: override this method")
+
+  def run_baker(self, X, extra_points = 3, order = 2):
+    raise Exception("abstract base class: override this method")
+
+  def __str__(self):
+    r = ''
+    assert( len(self.verts) == len(self.q) )
+    for c, i in enumerate(zip(self.verts, self.q)):
+      r += "%d %r: %0.4f" % (c,i[0], i[1])
+      face_str = ", ".join([str(f.name) for f in self.faces_for_vert[c]])
+      r += "  faces: [%s]" % face_str
+      r += "\n"
+    if self.faces:
+      for v in self.faces.itervalues():
+        r += "%s\n" % v
+    return r
+
+
+
 class delaunay_grid(grid):
-  facet_re = re.compile(r'''
+  face_re = re.compile(r'''
-                            -\s+(?P<facet>f\d+).*?
+                            -\s+(?P<face>f\d+).*?
                             vertices:\s(?P<verts>.*?)\n.*?
-                            neighboring\s facets:\s+(?P<neigh>[\sf\d]*)
+                            neighboring\s faces:\s+(?P<neigh>[\sf\d]*)
                          ''', re.S|re.X)
 
   point_re = re.compile(r'''
@@ -43,21 +135,10 @@ class delaunay_grid(grid):
                             (p\d+)
                          ''', re.S|re.X)
 
-  def __init__(self, points, q):
+  def __init__(self, verts, q):
-    grid.__init__(self, points,q)
+    grid.__init__(self, verts,q)
 
 
-  def create_mesh(self, indicies):
-    """
-      this function takes a list of indicies, and then creates
-      and returns a grid object (collection of points and q).
-
-      note: the input is indicies, the grid contains points
-    """
-    p = [self.points[i] for i in indicies]
-    q = [self.q[i]      for i in indicies]
-    return grid(p, q)
-
   def get_containing_simplex(self, X):
     if not self.faces:
       self.construct_connectivity()
@@ -66,38 +147,38 @@ class delaunay_grid(grid):
     (dist, indicies) = self.tree.query(X, 2)
     closest_point = indicies[0]
 
-    smblog.debug('X: %s' % X)
+    log.debug('X: %s' % X)
-    smblog.debug('point index: %d' % closest_point)
+    log.debug('point index: %d' % closest_point)
-    smblog.debug('actual point %s' % self.points[closest_point])
+    log.debug('actual point %s' % self.verts[closest_point])
-    smblog.debug('distance = %0.4f' % dist[0])
+    log.debug('distance = %0.4f' % dist[0])
 
     simplex = None
-    checked_facets = []
+    checked_faces = []
-    facets_to_check = self.faces_for_vert[closest_point]
+    faces_to_check = self.faces_for_vert[closest_point]
 
     attempts = 0
-    while not simplex and facets_to_check:
+    while not simplex and faces_to_check:
       attempts += 1
       # if attempts > 20:
-      #   raise smberror("probably recursing to many times")
+      #   raise Exception("probably recursing to many times")
-      cur_facet = facets_to_check.pop(0)
+      cur_face = faces_to_check.pop(0)
-      checked_facets.append(cur_facet)
+      checked_faces.append(cur_face)
 
-      if cur_facet.contains(X, self):
+      if cur_face.contains(X, self):
-        simplex = cur_facet
+        simplex = cur_face
         continue
 
       new_facest = []
-      for neighbor in cur_facet.neighbors:
+      for neighbor in cur_face.neighbors:
-        if (neighbor not in checked_facets) and (neighbor not in facets_to_check):
+        if (neighbor not in checked_faces) and (neighbor not in faces_to_check):
-          facets_to_check.append(neighbor)
+          faces_to_check.append(neighbor)
 
     if not simplex:
       raise AssertionError('no containing simplex found')
 
     R = self.create_mesh(simplex.verts)
 
-    smblog.debug('total attempts before finding simplex: %d' % attempts)
+    log.debug('total attempts before finding simplex: %d' % attempts)
     return R
 
 
@@ -108,33 +189,33 @@ class delaunay_grid(grid):
       R is a grid object that is supposedly a containing simplex
         around point X (it tends not to be)
 
-      S is S_j from baker's paper : some points from all point that are not the simplex
+      S is S_j from baker's paper : some verts from all point that are not the simplex
     """
-    smblog.debug(inspect.stack()[1][3])
+    log.debug(inspect.stack()[1][3])
-    smblog.debug("extra points: %d" % extra_points)
+    log.debug("extra verts: %d" % extra_points)
-    smblog.debug("simplex size: %d" % simplex_size)
+    log.debug("simplex size: %d" % simplex_size)
 
     r_mesh = self.get_containing_simplex(X)
-    # smblog.debug("R:\n%s" % r_mesh)
+    # log.debug("R:\n%s" % r_mesh)
 
-    # and some UNIQUE extra points
+    # and some UNIQUE extra verts
     (dist, indicies) = self.tree.query(X, simplex_size + extra_points)
 
     unique_indicies = []
     for index in indicies:
-      if self.points[index] not in r_mesh.points:
+      if self.verts[index] not in r_mesh.verts:
         unique_indicies.append(index)
 
-    smblog.debug("indicies: %s" % ",".join([str(i) for i in indicies]))
+    log.debug("indicies: %s" % ",".join([str(i) for i in indicies]))
-    smblog.debug("indicies: %s" % ",".join([str(i) for i in unique_indicies]))
+    log.debug("indicies: %s" % ",".join([str(i) for i in unique_indicies]))
     s_mesh = self.create_mesh(unique_indicies)# indicies[simplex_size:])
 
     # TODO: eventually remove this test:
-    for point in s_mesh.points:
+    for point in s_mesh.verts:
-      if point in r_mesh.points:
+      if point in r_mesh.verts:
-        smblog.error("ERROR")
+        log.error("ERROR")
-        smblog.error("\n%s\nin\n%s" % (point, r_mesh))
+        log.error("\n%s\nin\n%s" % (point, r_mesh))
-        raise smberror("repeating point S and R")
+        raise Exception("repeating point S and R")
 
     return (r_mesh, s_mesh)
 
@@ -143,16 +224,16 @@ class delaunay_grid(grid):
       this returns two grid objects: R and S.
 
       this function differes from the get_simplex_and_nearest_points
-      function in that it builds up the extra points based on
+      function in that it builds up the extra verts based on
       connectivity information, not just nearest-neighbor.
       in theory, this will work much better for situations like
-      points near a short edge in a boundary layer cell where the
+      verts near a short edge in a boundary layer cell where the
-      nearest points would all be colinear
+      nearest verts would all be colinear
 
       also, it guarantees that we find a containing simplex
 
       R is a grid object that is the (a) containing simplex around point X
-      S is a connectivity-based nearest-neighbor lookup, limited to 3 extra points
+      S is a connectivity-based nearest-neighbor lookup, limited to 3 extra verts
     """
     if not self.faces:
       self.construct_connectivity()
@@ -161,9 +242,9 @@ class delaunay_grid(grid):
     (dist, indicies) = self.tree.query(X, 2)
 
     simplex = None
-    for facet in self.faces_for_vert[indicies[0]]:
+    for face in self.faces_for_vert[indicies[0]]:
-      if facet.contains(X, self):
+      if face.contains(X, self):
-        simplex = facet
+        simplex = face
         break
 
     if not simplex:
@@ -184,11 +265,11 @@ class delaunay_grid(grid):
 
     try:
       (R, S) = self.get_simplex_and_nearest_points(X)
-      if not contains(X, R.points):
+      if not contains(X, R.verts):
-        raise smberror("run_baker with get_simplex_and_nearest_points returned non-containing simplex")
+        raise Exception("run_baker with get_simplex_and_nearest_points returned non-containing simplex")
       answer = run_baker(X, R, S, order)
-    except smberror, e:
+    except Exception, e:
-      smblog.error("caught error: %s, trying with connectivity-based mesh" % e)
+      log.error("caught error: %s, trying with connectivity-based mesh" % e)
       (R, S) = self.get_points_conn(X)
       answer = run_baker(X, R, S, order)
 
@@ -202,43 +283,29 @@ class delaunay_grid(grid):
 
       this is part of the __init__ for a rect_grid, but can be called from any grid object
     """
-    smblog.debug('start')
+    log.debug('start')
     qdelaunay_string = get_qdelaunay_dump_str(self)
-    facet_to_facets = []
+    face_to_faces = []
-    for matcher in grid.facet_re.finditer(qdelaunay_string):
+    for matcher in grid.face_re.finditer(qdelaunay_string):
       d = matcher.groupdict()
 
-      facet_name          = d['facet']
+      face_name          = d['face']
       verticies           = d['verts']
-      neighboring_facets  = d['neigh']
+      neighboring_faces  = d['neigh']
 
-      cur_face = face(facet_name)
+      cur_face = face(face_name)
-      self.faces[facet_name] = cur_face
+      self.faces[face_name] = cur_face
 
       for v in grid.vert_re.findall(verticies):
         vertex_index = int(v[1:])
         cur_face.add_vert(vertex_index)
         self.faces_for_vert[vertex_index].append(cur_face)
 
-      nghbrs = [(facet_name, i) for i in  neighboring_facets.split()]
+      nghbrs = [(face_name, i) for i in  neighboring_faces.split()]
-      facet_to_facets.extend(nghbrs)
+      face_to_faces.extend(nghbrs)
 
-    for rel in facet_to_facets:
+    for rel in face_to_faces:
       if rel[1] in self.faces:
         self.faces[rel[0]].add_neighbor(self.faces[rel[1]])
 
-    smblog.debug('end')
+    log.debug('end')
-
-
-  def __str__(self):
-    r = ''
-    assert( len(self.points) == len(self.q) )
-    for c, i in enumerate(zip(self.points, self.q)):
-      r += "%d %r: %0.4f" % (c,i[0], i[1])
-      facet_str = ", ".join([f.name for f in self.faces_for_vert[c]])
-      r += "  faces: [%s]" % facet_str
-      r += "\n"
-    if self.faces:
-      for v in self.faces.itervalues():
-        r += "%s\n" % v
-    return r