From a2d7b3f063123788d8149afa6601e6d4563b310a Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Mon, 8 Mar 2010 14:30:22 -0700 Subject: [PATCH] updated library files so that i don't have to edit the __init__.py files --HG-- rename : lib/baker/__init__.py => lib/baker/baker.py rename : lib/grid/__init__.py => lib/grid/grid.py --- .hgignore | 1 + lib/baker/__init__.py | 164 +----------------------- lib/baker/baker.py | 163 ++++++++++++++++++++++++ lib/grid/__init__.py | 288 +----------------------------------------- lib/grid/grid.py | 287 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 453 insertions(+), 450 deletions(-) create mode 100644 lib/baker/baker.py create mode 100755 lib/grid/grid.py diff --git a/.hgignore b/.hgignore index 205d1c5..43b0880 100644 --- a/.hgignore +++ b/.hgignore @@ -1,2 +1,3 @@ \.pyc$ \.blend$ +egg diff --git a/lib/baker/__init__.py b/lib/baker/__init__.py index 57f52d8..1cca1de 100644 --- a/lib/baker/__init__.py +++ b/lib/baker/__init__.py @@ -1,163 +1 @@ -import numpy as np -import sys - -from baker.tools import smberror - -def get_phis(X, R): - """ - The get_phis function is used to get barycentric coordonites for a point on a triangle. - - X -- the destination point (2D) - X = [0,0] - r -- the three points that make up the triangular simplex (2D) - r = [[-1, -1], [0, 2], [1, -1]] - - this will return [0.333, 0.333, 0.333] - """ - - # 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]], - ]) - b = np.array([ 1, - X[0], - X[1] - ]) - try: - phi = np.linalg.solve(A,b) - except: - print >> sys.stderr, "warning: get_phis: calculation of phis yielded a linearly dependant system" - raise smberror('get_phis') - phi = np.dot(np.linalg.pinv(A), b) - - return phi - -def get_phis_3D(X, r): - """ - The get_phis function is used to get barycentric coordonites for a point on a triangle. - - X -- the destination point (3D) - X = [0,0,0] - r -- the four points that make up the tetrahedron (3D) - r = [[-1, -1], [0, 2], [1, -1]] - - this will return [0.333, 0.333, 0.333] - """ - - # baker: eq 7 - A = np.array([ - [ 1, 1, 1, 1 ], - [r[0][0], r[1][0], r[2][0], r[3][0]], - [r[0][1], r[1][1], r[2][1], r[3][1]], - [r[0][2], r[1][2], r[2][2], r[3][2]], - ]) - b = np.array([ 1, - X[0], - X[1], - X[2] - ]) - try: - phi = np.linalg.solve(A,b) - except: - print >> sys.stderr, "warning: get_phis_3D: calculation of phis yielded a linearly dependant system" - phi = np.dot(np.linalg.pinv(A), b) - - return phi - - -def qlinear(X, R): - """ - this calculates the linear portion of q from X to r - - also, this is baker eq 3 - - X = destination point - R = simplex points - q = CFD quantities of interest at the simplex points - """ - - phis = get_phis(X, R.points) - qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) - return phis, qlin - -def qlinear_3D(X, R, q): - """ - this calculates the linear portion of q from X to r - - X = destination point - R = simplex points - q = CFD quantities of interest at the simplex points(R) - """ - - phis = get_phis_3D(X, R) - qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) - return phis, qlin - -def run_baker(X, R, S): - """ - This is the main function to call to get an interpolation to X from the input meshes - - X -- the destination point (2D) - X = [0,0] - - R = Simplex - - S = extra points - - """ - - # calculate values only for the triangle - phi, qlin = qlinear (X, R) - - if len(S.points) == 0: - answer = { - 'a': None, - 'b': None, - 'c': None, - 'qlin': qlin, - 'error': None, - 'final': None, - } - return answer - - B = [] # baker eq 9 - w = [] # baker eq 11 - - for (s, q) in zip(S.points, S.q): - cur_phi, cur_qlin = qlinear(s, R) - (phi1, phi2, phi3) = cur_phi - - B.append([phi1 * phi2, phi2 * phi3, phi3 * phi1]) - w.append(q - cur_qlin) - - B = np.array(B) - w = np.array(w) - - A = np.dot(B.T, B) - b = np.dot(B.T, w) - - # baker solve eq 10 - try: - (a, b, c) = np.linalg.solve(A,b) - except: - print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv" - (a, b, c) = np.dot(np.linalg.pinv(A), b) - - error_term = a * phi[0] * phi[1]\ - + b * phi[1] * phi[2]\ - + c * phi[2] * phi[0] - - q_final = qlin + error_term - - answer = { - 'a': a, - 'b': b, - 'c': c, - 'qlin': qlin, - 'error': error_term, - 'final': q_final, - } - - return answer +from baker import * diff --git a/lib/baker/baker.py b/lib/baker/baker.py new file mode 100644 index 0000000..d2289fb --- /dev/null +++ b/lib/baker/baker.py @@ -0,0 +1,163 @@ +import numpy as np +import sys + +from tools import smberror + +def get_phis(X, R): + """ + The get_phis function is used to get barycentric coordonites for a point on a triangle. + + X -- the destination point (2D) + X = [0,0] + r -- the three points that make up the triangular simplex (2D) + r = [[-1, -1], [0, 2], [1, -1]] + + this will return [0.333, 0.333, 0.333] + """ + + # 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]], + ]) + b = np.array([ 1, + X[0], + X[1] + ]) + try: + phi = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: get_phis: calculation of phis yielded a linearly dependant system" + raise smberror('get_phis') + phi = np.dot(np.linalg.pinv(A), b) + + return phi + +def get_phis_3D(X, r): + """ + The get_phis function is used to get barycentric coordonites for a point on a triangle. + + X -- the destination point (3D) + X = [0,0,0] + r -- the four points that make up the tetrahedron (3D) + r = [[-1, -1], [0, 2], [1, -1]] + + this will return [0.333, 0.333, 0.333] + """ + + # baker: eq 7 + A = np.array([ + [ 1, 1, 1, 1 ], + [r[0][0], r[1][0], r[2][0], r[3][0]], + [r[0][1], r[1][1], r[2][1], r[3][1]], + [r[0][2], r[1][2], r[2][2], r[3][2]], + ]) + b = np.array([ 1, + X[0], + X[1], + X[2] + ]) + try: + phi = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: get_phis_3D: calculation of phis yielded a linearly dependant system" + phi = np.dot(np.linalg.pinv(A), b) + + return phi + + +def qlinear(X, R): + """ + this calculates the linear portion of q from X to r + + also, this is baker eq 3 + + X = destination point + R = simplex points + q = CFD quantities of interest at the simplex points + """ + + phis = get_phis(X, R.points) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)]) + return phis, qlin + +def qlinear_3D(X, R, q): + """ + this calculates the linear portion of q from X to r + + X = destination point + R = simplex points + q = CFD quantities of interest at the simplex points(R) + """ + + phis = get_phis_3D(X, R) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) + return phis, qlin + +def run_baker(X, R, S): + """ + This is the main function to call to get an interpolation to X from the input meshes + + X -- the destination point (2D) + X = [0,0] + + R = Simplex + + S = extra points + + """ + + # calculate values only for the triangle + phi, qlin = qlinear (X, R) + + if len(S.points) == 0: + answer = { + 'a': None, + 'b': None, + 'c': None, + 'qlin': qlin, + 'error': None, + 'final': None, + } + return answer + + B = [] # baker eq 9 + w = [] # baker eq 11 + + for (s, q) in zip(S.points, S.q): + cur_phi, cur_qlin = qlinear(s, R) + (phi1, phi2, phi3) = cur_phi + + B.append([phi1 * phi2, phi2 * phi3, phi3 * phi1]) + w.append(q - cur_qlin) + + B = np.array(B) + w = np.array(w) + + A = np.dot(B.T, B) + b = np.dot(B.T, w) + + # baker solve eq 10 + try: + (a, b, c) = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv" + (a, b, c) = np.dot(np.linalg.pinv(A), b) + + error_term = a * phi[0] * phi[1]\ + + b * phi[1] * phi[2]\ + + c * phi[2] * phi[0] + + q_final = qlin + error_term + + answer = { + 'a': a, + 'b': b, + 'c': c, + 'qlin': qlin, + 'error': error_term, + 'final': q_final, + } + + return answer diff --git a/lib/grid/__init__.py b/lib/grid/__init__.py index ead5d24..985eee0 100755 --- a/lib/grid/__init__.py +++ b/lib/grid/__init__.py @@ -1,287 +1 @@ -#!/usr/bin/python - -import sys -import re -from collections import defaultdict - -import numpy as np -import scipy.spatial - -from baker import run_baker, get_phis -from baker.tools import exact_func, smberror -from grid.smcqdelaunay import * - -class face(object): - def __init__(self, name): - self.name = name - self.verts = [] - 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 contains(self, X, grid): - R = [grid.points[i] for i in self.verts] - - phis = get_phis(X, R) - - r = True - if [i for i in phis if i < 0.0]: - r = False - return r - - def __str__(self): - neighbors = [i.name for i in self.neighbors] - return '%s: verts: %s neighbors: [%s]' %\ - ( - self.name, - self.verts, - ", ".join(neighbors) - ) - - - - - -class grid(object): - facet_re = re.compile(r''' - -\s+(?Pf\d+).*? - vertices:\s(?P.*?)\n.*? - neighboring\s facets:\s+(?P[\sf\d]*) - ''', re.S|re.X) - - point_re = re.compile(r''' - -\s+(?Pp\d+).*? - neighbors:\s+(?P[\sf\d]*) - ''', re.S|re.X) - - vert_re = re.compile(r''' - (p\d+) - ''', re.S|re.X) - - - def __init__(self, points, q): - """ - this thing eats two pre-constructed arrays of stuff: - points = array of arrays (i will convert to numpy.array) - [[x0,y0], [x1,y1], ...] - q = array (1D) of important values - """ - - self.points = np.array(points) - self.q = np.array(q) - self.tree = scipy.spatial.KDTree(self.points) - 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. - - R is a grid object that is the (a) containing simplex around point X - S is S_j from baker's paper : some points from all point that are not the simplex - """ - (dist, indicies) = self.tree.query(X, 3 + extra_points) - - - # get the containing simplex - r_mesh = self.create_mesh(indicies[:simplex_size]) - # and some extra points - 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. - - this function differes from the get_simplex_and_nearest_points - function in that it builds up the extra points 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 - nearest points would all be colinear - - 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() - - # get closest point - (dist, indicies) = self.tree.query(X, 2) - - simplex = None - for i in self.facets_for_point[indicies[0]]: - if i.contains(X, self): - simplex = i - break - - if not simplex: - raise AssertionError('no containing simplex found') - - R = self.create_mesh(simplex.verts) - - - s = [] - for c,i in enumerate(simplex.neighbors): - s.extend([guy for guy in i.verts if not guy in simplex.verts]) - S = self.create_mesh(s) - - return R, S - - def run_baker(self, X): - answer = None - - try: - (R, S) = self.get_simplex_and_nearest_points(X) - answer = run_baker(X, R, S) - except smberror as e: - print "caught error: %s, trying with connectivity-based mesh" % e - (R, S) = self.get_points_conn(X) - answer = run_baker(X, R, S) - - return answer - - - - def construct_connectivity(self): - """ - a call to this method prepares the internal connectivity structure. - - this is part of the __init__ for a simple_rect_grid, but can be called from any grid object - """ - qdelaunay_string = get_qdelaunay_dump_str(self) - facet_to_facets = [] - for matcher in grid.facet_re.finditer(qdelaunay_string): - d = matcher.groupdict() - - facet_name = d['facet'] - verticies = d['verts'] - neighboring_facets = d['neigh'] - - cur_face = face(facet_name) - self.faces[facet_name] = cur_face - - for v in grid.vert_re.findall(verticies): - vertex_index = int(v[1:]) - cur_face.add_vert(vertex_index) - self.facets_for_point[vertex_index].append(cur_face) - - nghbrs = [(facet_name, i) for i in neighboring_facets.split()] - facet_to_facets.extend(nghbrs) - - for rel in facet_to_facets: - if rel[1] in self.faces: - self.faces[rel[0]].add_neighbor(self.faces[rel[1]]) - - # for matcher in grid.point_re.finditer(qdelaunay_string): - # d = matcher.groupdict() - - # point = d['point'] - # neighboring_facets = d['neigh'] - - # self.facets_for_point[int(point[1:])] = [i for i in neighboring_facets.split() if i in self.faces] - - def for_qhull_generator(self): - """ - this returns a generator that should be fed into qdelaunay - """ - - yield '2'; - yield '%d' % len(self.points) - - for p in self.points: - yield "%f %f" % (p[0], p[1]) - - def for_qhull(self): - """ - this returns a single string that should be fed into qdelaunay - """ - r = '2\n' - r += '%d\n' % len(self.points) - for p in self.points: - r += "%f %f\n" % (p[0], p[1]) - return r - - 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.facets_for_point[c]]) - r += " faces: [%s]" % facet_str - r += "\n" - if self.faces: - for v in self.faces.itervalues(): - r += "%s\n" % v - return r - -class simple_rect_grid(grid): - def __init__(self, xres = 5, yres = 5): - xmin = -1.0 - xmax = 1.0 - xspan = xmax - xmin - xdel = xspan / float(xres - 1) - - ymin = -1.0 - ymay = 1.0 - yspan = ymay - ymin - ydel = yspan / float(yres - 1) - - - points = [] - q = [] - for x in xrange(xres): - cur_x = xmin + (x * xdel) - for y in xrange(yres): - cur_y = ymin + (y * ydel) - points.append([cur_x, cur_y]) - q.append(exact_func(cur_x, cur_y)) - grid.__init__(self, points, q) - self.construct_connectivity() - - - -class simple_random_grid(simple_rect_grid): - def __init__(self, num_points = 10): - points = [] - q = [] - - r = np.random - - for i in xrange(num_points): - cur_x = r.rand() - cur_y = r.rand() - - points.append([cur_x, cur_y]) - q.append(exact_func(cur_x, cur_y)) - grid.__init__(self, points, q) - - self.points = np.array(self.points) - self.q = np.array(self.q) - - -if __name__ == '__main__': - try: - resolution = int(sys.argv[1]) - except: - resolution = 10 - g = simple_rect_grid(resolution, resolution) - print g.for_qhull() +from grid import * diff --git a/lib/grid/grid.py b/lib/grid/grid.py new file mode 100755 index 0000000..6cbd413 --- /dev/null +++ b/lib/grid/grid.py @@ -0,0 +1,287 @@ +#!/usr/bin/python + +import sys +import re +from collections import defaultdict + +import numpy as np +import scipy.spatial + +from baker import run_baker, get_phis +from baker.tools import exact_func, smberror +from smcqdelaunay import * + +class face(object): + def __init__(self, name): + self.name = name + self.verts = [] + 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 contains(self, X, grid): + R = [grid.points[i] for i in self.verts] + + phis = get_phis(X, R) + + r = True + if [i for i in phis if i < 0.0]: + r = False + return r + + def __str__(self): + neighbors = [i.name for i in self.neighbors] + return '%s: verts: %s neighbors: [%s]' %\ + ( + self.name, + self.verts, + ", ".join(neighbors) + ) + + + + + +class grid(object): + facet_re = re.compile(r''' + -\s+(?Pf\d+).*? + vertices:\s(?P.*?)\n.*? + neighboring\s facets:\s+(?P[\sf\d]*) + ''', re.S|re.X) + + point_re = re.compile(r''' + -\s+(?Pp\d+).*? + neighbors:\s+(?P[\sf\d]*) + ''', re.S|re.X) + + vert_re = re.compile(r''' + (p\d+) + ''', re.S|re.X) + + + def __init__(self, points, q): + """ + this thing eats two pre-constructed arrays of stuff: + points = array of arrays (i will convert to numpy.array) + [[x0,y0], [x1,y1], ...] + q = array (1D) of important values + """ + + self.points = np.array(points) + self.q = np.array(q) + self.tree = scipy.spatial.KDTree(self.points) + 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. + + R is a grid object that is the (a) containing simplex around point X + S is S_j from baker's paper : some points from all point that are not the simplex + """ + (dist, indicies) = self.tree.query(X, 3 + extra_points) + + + # get the containing simplex + r_mesh = self.create_mesh(indicies[:simplex_size]) + # and some extra points + 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. + + this function differes from the get_simplex_and_nearest_points + function in that it builds up the extra points 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 + nearest points would all be colinear + + 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() + + # get closest point + (dist, indicies) = self.tree.query(X, 2) + + simplex = None + for i in self.facets_for_point[indicies[0]]: + if i.contains(X, self): + simplex = i + break + + if not simplex: + raise AssertionError('no containing simplex found') + + R = self.create_mesh(simplex.verts) + + + s = [] + for c,i in enumerate(simplex.neighbors): + s.extend([guy for guy in i.verts if not guy in simplex.verts]) + S = self.create_mesh(s) + + return R, S + + def run_baker(self, X): + answer = None + + try: + (R, S) = self.get_simplex_and_nearest_points(X) + answer = run_baker(X, R, S) + except smberror as e: + print "caught error: %s, trying with connectivity-based mesh" % e + (R, S) = self.get_points_conn(X) + answer = run_baker(X, R, S) + + return answer + + + + def construct_connectivity(self): + """ + a call to this method prepares the internal connectivity structure. + + this is part of the __init__ for a simple_rect_grid, but can be called from any grid object + """ + qdelaunay_string = get_qdelaunay_dump_str(self) + facet_to_facets = [] + for matcher in grid.facet_re.finditer(qdelaunay_string): + d = matcher.groupdict() + + facet_name = d['facet'] + verticies = d['verts'] + neighboring_facets = d['neigh'] + + cur_face = face(facet_name) + self.faces[facet_name] = cur_face + + for v in grid.vert_re.findall(verticies): + vertex_index = int(v[1:]) + cur_face.add_vert(vertex_index) + self.facets_for_point[vertex_index].append(cur_face) + + nghbrs = [(facet_name, i) for i in neighboring_facets.split()] + facet_to_facets.extend(nghbrs) + + for rel in facet_to_facets: + if rel[1] in self.faces: + self.faces[rel[0]].add_neighbor(self.faces[rel[1]]) + + # for matcher in grid.point_re.finditer(qdelaunay_string): + # d = matcher.groupdict() + + # point = d['point'] + # neighboring_facets = d['neigh'] + + # self.facets_for_point[int(point[1:])] = [i for i in neighboring_facets.split() if i in self.faces] + + def for_qhull_generator(self): + """ + this returns a generator that should be fed into qdelaunay + """ + + yield '2'; + yield '%d' % len(self.points) + + for p in self.points: + yield "%f %f" % (p[0], p[1]) + + def for_qhull(self): + """ + this returns a single string that should be fed into qdelaunay + """ + r = '2\n' + r += '%d\n' % len(self.points) + for p in self.points: + r += "%f %f\n" % (p[0], p[1]) + return r + + 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.facets_for_point[c]]) + r += " faces: [%s]" % facet_str + r += "\n" + if self.faces: + for v in self.faces.itervalues(): + r += "%s\n" % v + return r + +class simple_rect_grid(grid): + def __init__(self, xres = 5, yres = 5): + xmin = -1.0 + xmax = 1.0 + xspan = xmax - xmin + xdel = xspan / float(xres - 1) + + ymin = -1.0 + ymay = 1.0 + yspan = ymay - ymin + ydel = yspan / float(yres - 1) + + + points = [] + q = [] + for x in xrange(xres): + cur_x = xmin + (x * xdel) + for y in xrange(yres): + cur_y = ymin + (y * ydel) + points.append([cur_x, cur_y]) + q.append(exact_func(cur_x, cur_y)) + grid.__init__(self, points, q) + self.construct_connectivity() + + + +class simple_random_grid(simple_rect_grid): + def __init__(self, num_points = 10): + points = [] + q = [] + + r = np.random + + for i in xrange(num_points): + cur_x = r.rand() + cur_y = r.rand() + + points.append([cur_x, cur_y]) + q.append(exact_func(cur_x, cur_y)) + grid.__init__(self, points, q) + + self.points = np.array(self.points) + self.q = np.array(self.q) + + +if __name__ == '__main__': + try: + resolution = int(sys.argv[1]) + except: + resolution = 10 + g = simple_rect_grid(resolution, resolution) + print g.for_qhull()