From bcf1622f3580ce995b9e723eca4d0ae6b093e86f Mon Sep 17 00:00:00 2001
From: Stephen Mardson McQuay <stephen@mardson.com>
Date: Fri, 29 Oct 2010 12:43:18 -0600
Subject: [PATCH] added a generic get simplex and points function to the grid

---
 interp/grid/__init__.py | 42 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 35 insertions(+), 7 deletions(-)

diff --git a/interp/grid/__init__.py b/interp/grid/__init__.py
index 692f9ed..c8d134c 100644
--- a/interp/grid/__init__.py
+++ b/interp/grid/__init__.py
@@ -28,16 +28,16 @@ class grid(object):
 
     if verts:
       self.verts = np.array(verts)
+      self.tree           = scipy.spatial.KDTree(self.verts)
     if q:
       self.q = np.array(q)
 
-    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()
+      raise Exception("face connectivity is not set up")
 
     # get closest point
     (dist, indicies) = self.tree.query(X, 2)
@@ -86,7 +86,7 @@ class grid(object):
 
       note: the input is indicies, the grid contains verts
     """
-    p = [self.verts[i] for i in indicies]
+    p = [self.verts[i]  for i in indicies]
     q = [self.q[i]      for i in indicies]
     return grid(p, q)
 
@@ -95,20 +95,48 @@ class grid(object):
       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)
+        around point X (TODO: 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")
+    log.debug("extra verts: %d" % extra_points)
+    log.debug("simplex size: %d" % simplex_size)
+
+    r_mesh = self.get_containing_simplex(X)
+
+    # and some UNIQUE extra verts
+    (dist, indicies) = self.tree.query(X, simplex_size + extra_points)
+    log.info("extra indicies: %s" % indicies)
+
+    unique_indicies = []
+    for index in indicies:
+      close_point_in_R = False
+      for rvert in r_mesh.verts:
+        if all(rvert == self.verts[index]):
+          close_point_in_R = True
+          break
+
+      if not close_point_in_R:
+        unique_indicies.append(index)
+      else:
+        log.debug('throwing out %s: %s' % (index, self.verts[index]))
+
+    log.debug("indicies: %s" % indicies)
+    log.debug("unique indicies: %s" % unique_indicies)
+    s_mesh = self.create_mesh(unique_indicies)
+
+    return (r_mesh, s_mesh)
 
   def run_baker(self, X, extra_points = 3, order = 2):
-    raise Exception("abstract base class: override this method")
+    (R, S) = self.get_simplex_and_nearest_points(X, extra_points)
+    answer = run_baker(X, R, S, order)
+    return answer
 
   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])
+      r += "%d vert(%s): q(%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"