diff --git a/bin/grid_regular.py b/bin/grid_regular.py index 22772ef..23efd5c 100755 --- a/bin/grid_regular.py +++ b/bin/grid_regular.py @@ -17,5 +17,6 @@ if __name__ == '__main__': g = simple_rect_grid(resolution, resolution) print g - print g.get_points_conn([0.4, 0]) + X = [0.4, 0.001] + print g.get_points_conn(X) open(qfile, 'w').write(g.for_qhull()) diff --git a/lib/baker.py b/lib/baker.py index 4adc1a2..89ec8e4 100644 --- a/lib/baker.py +++ b/lib/baker.py @@ -1,7 +1,7 @@ import numpy as np import sys -def get_phis(X, r): +def get_phis(X, R): """ The get_phis function is used to get barycentric coordonites for a point on a triangle. @@ -16,8 +16,8 @@ def get_phis(X, r): # baker: eq 7 A = np.array([ [ 1, 1, 1], - [r[0][0], r[1][0], r[2][0]], - [r[0][1], r[1][1], r[2][1]], + [R[0][0], R[1][0], R[2][0]], + [R[0][1], R[1][1], R[2][1]], ]) b = np.array([ 1, X[0], diff --git a/lib/grid.py b/lib/grid.py index fd4a192..11226e1 100755 --- a/lib/grid.py +++ b/lib/grid.py @@ -18,28 +18,21 @@ class face(object): self.neighbors = [] def add_vert(self, v): + """ + v should be an index into grid.points + """ self.verts.append(v) def add_neighbor(self, n): + """ + reference to another face object + """ self.neighbors.append(n) - def __str__(self): - neighbors = [i.name for i in self.neighbors] - return '%s: points: %s neighbors: [%s]' %\ - ( - self.name, - self.verts, - ", ".join(neighbors) - ) + def contains(self, X, grid): + R = [grid.points[i] for i in self.verts] - - -class simplex(object): - def __init__(self, points): - self.R = points - - def contains(self, X): - phis = get_phis(X, self.R) + phis = get_phis(X, R) r = True if [i for i in phis if i < 0.0]: @@ -47,8 +40,13 @@ class simplex(object): return r def __str__(self): - print type(self.R) - return "simplex: %d: [%s]" % (len(self.R), ", ".join([str(i) for i in self.R]),) + neighbors = [i.name for i in self.neighbors] + return '%s: verts: %s neighbors: [%s]' %\ + ( + self.name, + self.verts, + ", ".join(neighbors) + ) @@ -85,6 +83,13 @@ class grid(object): self.faces = {} self.facets_for_point = defaultdict(list) + + def create_mesh(self, indicies): + p = [self.points[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. @@ -96,35 +101,52 @@ class grid(object): # get the containing simplex - R = [self.points[i] for i in indicies[:simplex_size] ] - Rq = [self.q[i] for i in indicies[:simplex_size] ] - r_mesh = grid(R, Rq) - + r_mesh = self.create_mesh(indicies[:simplex_size]) # and some extra points - S = [self.points[i] for i in indicies[simplex_size:] ] - Sq = [self.q[i] for i in indicies[simplex_size:] ] - s_mesh = grid(S, Sq) + s_mesh = self.create_mesh(indicies[simplex_size:]) return (r_mesh, s_mesh) + def get_points_conn(self, X): + """ + this returns two grid objects: R and S. + + 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 + """ + if not self.faces: + self.construct_connectivity() + + (dist, indicies) = self.tree.query(X, 2) + smplx = None + for i in self.facets_for_point[indicies[0]]: + if i.contains(X, self): + smplx = i + break + + if not smplx: + raise AssertionError('no containing simplex found') + print "containing verts:", smplx.verts + smplx_set = set(smplx.verts) + R = self.create_mesh(smplx.verts) + + + s = [] + for c,i in enumerate(smplx.neighbors): + print "neighbor %d: %s: %s" % (c, i.name, i.verts) + st = set(i.verts) + s.append((st - smplx_set).pop()) +# s.append( [guy for guy in i.verts if not guy in smplx.verts]) + print s + S = self.create_mesh(s) + + print S + return R, S + def run_baker(self, X): (R, S) = self.get_simplex_and_nearest_points(X) return run_baker(X, R, S) - def get_points_conn(self, X): - if not self.faces: - self.construct_connectivity() - - (dist, indicies) = self.tree.query(X, 3) - print dist, indicies, [(self.points[i][0], self.points[i][1]) for i in indicies] - print [i.name for i in self.facets_for_point[indicies[0]]] - print self.facets_for_point[indicies[0]] - - smplx = simplex([self.points[i] for i in indicies]) - - if not smplx.contains(X): - raise AssertionError('simplex should contain point') - def construct_connectivity(self): """ a call to this method prepares the internal connectivity structure.