From a3d72ed2002945f608fc4c187f24e1091f458d72 Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Sun, 21 Feb 2010 20:42:32 -0700 Subject: [PATCH] wrote the method that returns grid objects based on connectivity I'm going to hurry and clean up this method, but for now it is verbose, and uses sets. i'm going to just check it in for future reference --- bin/grid_regular.py | 3 +- lib/baker.py | 6 +-- lib/grid.py | 100 +++++++++++++++++++++++++++----------------- 3 files changed, 66 insertions(+), 43 deletions(-) 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.