From 1fc988428d00fd7f099d5af73f668c8d26165157 Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Sat, 20 Mar 2010 17:10:02 -0600 Subject: [PATCH] putting together a simple test case so that I can test my quad/cubic interpolator --- bin/driver.py | 12 +- bin/generate_qhull.pl | 16 --- bin/grid_regular.py | 25 +++- bin/test.py | 27 ++-- lib/baker/__init__.py | 273 ++++++++++++++++++++++++++++++++++++ lib/baker/baker.py | 272 ----------------------------------- lib/baker/tools.py | 4 +- lib/grid/DD.py | 57 ++++---- lib/grid/DDD.py | 31 +++- lib/grid/__init__.py | 174 ++++++++++++++++++++++- lib/grid/grid.py | 186 ------------------------ lib/grid/simplex.py | 19 +-- test/quad.test.py | 54 +++++++ tools/blender/plot_cloud.py | 4 +- 14 files changed, 616 insertions(+), 538 deletions(-) delete mode 100755 bin/generate_qhull.pl delete mode 100644 lib/baker/baker.py mode change 100755 => 100644 lib/grid/__init__.py delete mode 100755 lib/grid/grid.py create mode 100755 test/quad.test.py diff --git a/bin/driver.py b/bin/driver.py index 66aabb0..6d2a35a 100755 --- a/bin/driver.py +++ b/bin/driver.py @@ -6,18 +6,18 @@ from optparse import OptionParser import numpy as np from baker.tools import rms, exact_func -from grid.DD import simple_random_grid, simple_rect_grid +from grid.DD import random_grid, rect_grid def get_mesh(source, destination, use_structured_grid = False): mesh_source = None mesh_dest = None if use_structured_grid: - mesh_source = simple_rect_grid(source, source) - mesh_dest = simple_rect_grid(destination, destination) + mesh_source = rect_grid(source, source) + mesh_dest = rect_grid(destination, destination) else: - mesh_source = simple_random_grid(source) - mesh_dest = simple_random_grid(destination) + mesh_source = random_grid(source) + mesh_dest = random_grid(destination) if not (mesh_dest and mesh_source): raise smberror('problem creating mesh objects') @@ -88,7 +88,7 @@ if __name__ == '__main__': errors.append(0) continue - exact = exact_func(X[0], X[1]) + exact = exact_func(X) if np.abs(exact - answer['final']) <= np.abs(exact - answer['qlin']): success += 1 diff --git a/bin/generate_qhull.pl b/bin/generate_qhull.pl deleted file mode 100755 index 437a6a9..0000000 --- a/bin/generate_qhull.pl +++ /dev/null @@ -1,16 +0,0 @@ -#!/usr/bin/perl - -use strict; -use warnings; - -my $output_filename = $ARGV[0] or die "no args"; -my $tmp_grid_file = "$output_filename.grid.txt"; - -system("grid.py > $tmp_grid_file"); - - -system("cat $tmp_grid_file | qdelaunay GD2 s > $output_filename.txt\n" ); -system("cat $tmp_grid_file | qdelaunay GD2 s Qt > ${output_filename}_qt.txt\n"); -system("cat $tmp_grid_file | qdelaunay GD2 s QJ > ${output_filename}_qj.txt\n"); - -unlink($tmp_grid_file); diff --git a/bin/grid_regular.py b/bin/grid_regular.py index f25087d..be72b66 100755 --- a/bin/grid_regular.py +++ b/bin/grid_regular.py @@ -3,9 +3,9 @@ import sys import pickle -from grid.DD import simple_rect_grid, simple_random_grid +from grid.DD import rect_grid, random_grid from baker import run_baker -from baker.tools import smberror +from baker.tools import exact_func, smberror qfile = '/tmp/grid_regular.txt' @@ -18,19 +18,23 @@ if __name__ == '__main__': rx = 4 ry = 4 * rx - source_mesh = simple_rect_grid(rx, ry) + source_mesh = rect_grid(rx, ry) # print source_mesh - X = [0.1, 0.1] + X = [0.1, 0.01] try: + print "trying to get simplex and nearest points using nearest points (kdtree)" (R, S) = source_mesh.get_simplex_and_nearest_points(X, extra_points=4) print "R for nearest-neighbor:\n", R print "S for nearest-neighbor:\n", S + print "trying to run baker" print run_baker(X, R, S) except smberror as e: print "caught error: %s" % e + + print "trying to get simplex and nearest points using connectivity scheme" (R, S) = source_mesh.get_points_conn(X) print "R for connectivity:\n", R print "S for connectivity:\n", S @@ -38,5 +42,16 @@ if __name__ == '__main__': print "repeating the above just using the grid object:" - print source_mesh.run_baker(X) + r = source_mesh.run_baker(X) + exact = exact_func(X) + print r + print 'exact', exact + print 'qlin' , r['qlin'] + print 'error', r['error'] + print 'final', r['final'] + + if abs(r['final'] - exact) <= abs(r['qlin'] - exact): + print "win" + else: + print "failure" open(qfile, 'w').write(source_mesh.for_qhull()) diff --git a/bin/test.py b/bin/test.py index 11786d0..0ff5e48 100755 --- a/bin/test.py +++ b/bin/test.py @@ -1,9 +1,9 @@ #!/usr/bin/python2.6 import sys -from grid.DDD import simple_random_grid +from grid.DDD import random_grid from baker import get_phis_3D, run_baker_3D -from baker.tools import exact_func_3D +from baker.tools import exact_func_3D, smberror try: total_points = int(sys.argv[1]) @@ -12,7 +12,7 @@ except: print total_points -g = simple_random_grid(total_points) +g = random_grid(total_points) open('/tmp/for_qhull.txt', 'w').write(g.for_qhull()) @@ -31,13 +31,16 @@ print "phi values (should all be positive): ", phis, sum(phis) if [i for i in phis if i < 0.0]: print "problems" -r = run_baker_3D(X, R, S) +try: + r = run_baker_3D(X, R, S) + print 'qlin' , r['qlin'] + print 'error', r['error'] + print 'final', r['final'] -print 'qlin' , r['qlin'] -print 'error', r['error'] -print 'final', r['final'] - -if abs(r['final'] - exact) <= abs(r['qlin'] - exact): - print "win" -else: - print "failure" + if abs(r['final'] - exact) <= abs(r['qlin'] - exact): + print "win" + else: + print "failure" +except smberror as e: + print e + print 'TAINT' diff --git a/lib/baker/__init__.py b/lib/baker/__init__.py index 1cca1de..b82ceac 100644 --- a/lib/baker/__init__.py +++ b/lib/baker/__init__.py @@ -1 +1,274 @@ from baker import * +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 containing 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 np.linalg.LinAlgError as e: + msg = "warning: get_phis: calculation of phis yielded a linearly dependant system (%s)" % e + # TODO: log this -- > print >> sys.stderr, msg + raise smberror(msg) + 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 tetrahedron. + + X -- the destination point (3D) + X = [0,0,0] + R -- the four points that make up the containing simplex, tetrahedron (3D) + R = [ + [0.0, 0.0, 1.0], + [0.94280904333606508, 0.0, -0.3333333283722672], + [-0.47140452166803232, 0.81649658244673617, -0.3333333283722672], + [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672], + ] + + this (should) will return [0.25, 0.25, 0.25, 0.25] + """ + + # 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 np.linalg.LinAlgError as e: + print >> sys.stderr, "warning: get_phis_3D: calculation of phis yielded a linearly dependant system", e + 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): + """ + 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.points) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.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 simplex triangle + phi, qlin = qlinear(X, R) + + if [i for i in phi if i <= 0.0]: + s = "this is not a containing simplex:\n" + s += " X: %s\n" % X + s += " R: %s\n" % R + s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi)) + print >> sys.stderr, s + raise smberror("simplex does not contain point") + + 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 np.linalg.LinAlgError as e: + print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e + (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 + +def run_baker_3D(X, R, S): + """ + This is the main function to call to get an interpolation to X from the input meshes + + X -- the destination point (3D) + X = [0,0,0] + + R = Simplex (4 points, contains X) + S = extra points (surrounding, in some manner, R and X, but not in R) + """ + + # calculate values only for the triangle + phi, qlin = qlinear_3D(X, R) + + if [i for i in phi if i <= 0.0]: + s = "this is not a containing simplex:\n" + s += " X: %s\n" % X + s += " R: %s\n" % R + s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi)) + print >> sys.stderr, s + raise smberror("not containing simplex") + + if len(S.points) == 0: + answer = { + 'a': None, + 'b': None, + 'c': None, + 'd': None, + 'e': None, + 'f': 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_3D(s, R) + (phi1, phi2, phi3, phi4) = cur_phi + + B.append( + [ + phi1 * phi2, + phi1 * phi3, + phi1 * phi4, + phi2 * phi3, + phi2 * phi4, + phi3 * phi4, + ] + ) + + 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, d, e, f) = np.linalg.solve(A,b) + except np.linalg.LinAlgError as e: + print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e + (a, b, c, d, e, f) = np.dot(np.linalg.pinv(A), b) + + error_term = a * phi[0] * phi[1]\ + + b * phi[0] * phi[2]\ + + c * phi[0] * phi[3]\ + + d * phi[1] * phi[2]\ + + e * phi[1] * phi[3]\ + + f * phi[2] * phi[3] + + q_final = qlin + error_term + + answer = { + 'a': a, + 'b': b, + 'c': c, + 'd': d, + 'e': e, + 'f': f, + 'qlin': qlin, + 'error': error_term, + 'final': q_final, + } + + return answer diff --git a/lib/baker/baker.py b/lib/baker/baker.py deleted file mode 100644 index 6f4c4cc..0000000 --- a/lib/baker/baker.py +++ /dev/null @@ -1,272 +0,0 @@ -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 np.linalg.LinAlgError as e: - print >> sys.stderr, "warning: get_phis: calculation of phis yielded a linearly dependant system", e - 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 = [ - [0.0, 0.0, 1.0], - [0.94280904333606508, 0.0, -0.3333333283722672], - [-0.47140452166803232, 0.81649658244673617, -0.3333333283722672], - [-0.47140452166803298, -0.81649658244673584, -0.3333333283722672], - ] - - this will return [0.25, 0.25, 0.25, 0.25] - """ - - # 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 np.linalg.LinAlgError as e: - print >> sys.stderr, "warning: get_phis_3D: calculation of phis yielded a linearly dependant system", e - 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): - """ - 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.points) - qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.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 simplex triangle - phi, qlin = qlinear(X, R) - - if [i for i in phi if i <= 0.0]: - s = "this is not a containing simplex:\n" - s += " X: %s\n" % X - s += " R: %s\n" % R - s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi)) - print >> sys.stderr, s - raise smberror("not containing simplex") - - 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 np.linalg.LinAlgError as e: - print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e - (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 - -def run_baker_3D(X, R, S): - """ - This is the main function to call to get an interpolation to X from the input meshes - - X -- the destination point (3D) - X = [0,0,0] - - R = Simplex (4 points, contains X) - S = extra points (surrounding, in some manner, R and X, but not in R) - """ - - # calculate values only for the triangle - phi, qlin = qlinear_3D(X, R) - - if [i for i in phi if i <= 0.0]: - s = "this is not a containing simplex:\n" - s += " X: %s\n" % X - s += " R: %s\n" % R - s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi)) - print >> sys.stderr, s - raise smberror("not containing simplex") - - if len(S.points) == 0: - answer = { - 'a': None, - 'b': None, - 'c': None, - 'd': None, - 'e': None, - 'f': 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_3D(s, R) - (phi1, phi2, phi3, phi4) = cur_phi - - B.append( - [ - phi1 * phi2, - phi1 * phi3, - phi1 * phi4, - phi2 * phi3, - phi2 * phi4, - phi3 * phi4, - ] - ) - - 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, d, e, f) = np.linalg.solve(A,b) - except np.linalg.LinAlgError as e: - print >> sys.stderr, "warning: run_baker: linear calculation went bad, resorting to np.linalg.pinv", e - (a, b, c, d, e, f) = np.dot(np.linalg.pinv(A), b) - - error_term = a * phi[0] * phi[1]\ - + b * phi[0] * phi[2]\ - + c * phi[0] * phi[3]\ - + d * phi[1] * phi[2]\ - + e * phi[1] * phi[3]\ - + f * phi[2] * phi[3] - - q_final = qlin + error_term - - answer = { - 'a': a, - 'b': b, - 'c': c, - 'd': d, - 'e': e, - 'f': f, - 'qlin': qlin, - 'error': error_term, - 'final': q_final, - } - - return answer diff --git a/lib/baker/tools.py b/lib/baker/tools.py index b7f2e34..de417fc 100644 --- a/lib/baker/tools.py +++ b/lib/baker/tools.py @@ -19,10 +19,12 @@ def rms(errors): r = np.sqrt(r / len(errors)) return r -def exact_func(x, y): +def exact_func(X): """ the exact function used from baker's article (for testing) """ + x = X[0] + y = X[0] return np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) def exact_func_3D(X): diff --git a/lib/grid/DD.py b/lib/grid/DD.py index 23d31c1..0959be2 100644 --- a/lib/grid/DD.py +++ b/lib/grid/DD.py @@ -1,32 +1,13 @@ -from grid import grid +from grid import grid as basegrid + from baker.tools import exact_func import numpy as np -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 grid(basegrid): + def __init__(self, points, q): + basegrid.__init__(self, points, q) def for_qhull_generator(self): """ @@ -49,8 +30,32 @@ class simple_rect_grid(grid): r += "%f %f\n" % (p[0], p[1]) return r +class rect_grid(grid): + def __init__(self, xres = 5, yres = 5): + xmin = -1.0 + xmax = 1.0 + xspan = xmax - xmin + xdel = xspan / float(xres - 1) -class simple_random_grid(simple_rect_grid): + 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 random_grid(rect_grid): def __init__(self, num_points = 10): points = [] q = [] @@ -62,7 +67,7 @@ class simple_random_grid(simple_rect_grid): cur_y = r.rand() points.append([cur_x, cur_y]) - q.append(exact_func(cur_x, cur_y)) + q.append( exact_func( (cur_x, cur_y) ) ) grid.__init__(self, points, q) self.points = np.array(self.points) diff --git a/lib/grid/DDD.py b/lib/grid/DDD.py index 01b3155..a829f12 100644 --- a/lib/grid/DDD.py +++ b/lib/grid/DDD.py @@ -1,10 +1,35 @@ -from grid import grid +from grid import grid as basegrid from baker.tools import exact_func_3D import numpy as np -class simple_rect_grid(grid): +class grid(basegrid): + def __init__(self, points, q): + basegrid.__init__(self, points, q) + + def for_qhull_generator(self): + """ + this returns a generator that should be fed into qdelaunay + """ + + yield '3'; + yield '%d' % len(self.points) + + for p in self.points: + yield "%f %f %f" % tuple(p) + + def for_qhull(self): + """ + this returns a single string that should be fed into qdelaunay + """ + r = '3\n' + r += '%d\n' % len(self.points) + for p in self.points: + r += "%f %f %f\n" % tuple(p) + return r + +class rect_grid(grid): def __init__(self, xres = 5, yres = 5, zres = 5): xmin = -1.0 xmax = 1.0 @@ -56,7 +81,7 @@ class simple_rect_grid(grid): r += "%f %f %f\n" % tuple(p) return r -class simple_random_grid(simple_rect_grid): +class random_grid(rect_grid): def __init__(self, num_points = 10): points = [] q = [] diff --git a/lib/grid/__init__.py b/lib/grid/__init__.py old mode 100755 new mode 100644 index 985eee0..ac0095f --- a/lib/grid/__init__.py +++ b/lib/grid/__init__.py @@ -1 +1,173 @@ -from grid import * +import sys +import re +from collections import defaultdict + +import numpy as np +import scipy.spatial + +from baker import run_baker +from baker.tools import exact_func, smberror +from simplex import face +from smcqdelaunay import * + + + +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, simplex_size + 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, e: + print >> sys.stderr, "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 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 __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 diff --git a/lib/grid/grid.py b/lib/grid/grid.py deleted file mode 100755 index 58a39be..0000000 --- a/lib/grid/grid.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/python - -import sys -import re -from collections import defaultdict - -import numpy as np -import scipy.spatial - -from baker import run_baker -from baker.tools import exact_func, smberror -from simplex import face -from smcqdelaunay import * - - - -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, simplex_size + 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, 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 __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 - - - - -if __name__ == '__main__': - try: - resolution = int(sys.argv[1]) - except: - resolution = 10 - g = simple_rect_grid(resolution, resolution) - print g.for_qhull() diff --git a/lib/grid/simplex.py b/lib/grid/simplex.py index b35d8a9..046496d 100644 --- a/lib/grid/simplex.py +++ b/lib/grid/simplex.py @@ -1,5 +1,15 @@ from baker import get_phis +TOL = 1e-3 + +def contains(X, R): + phis = get_phis(X, R) + r = True + if [i for i in phis if i < 0.0 - TOL]: + r = False + return r + + class face(object): def __init__(self, name): self.name = name @@ -19,14 +29,7 @@ class 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 + return contains(X, grid.points) def __str__(self): neighbors = [i.name for i in self.neighbors] diff --git a/test/quad.test.py b/test/quad.test.py new file mode 100755 index 0000000..df02a61 --- /dev/null +++ b/test/quad.test.py @@ -0,0 +1,54 @@ +#!/usr/bin/python + +import unittest + +from baker import run_baker +from baker.tools import exact_func + +from grid.DD import grid +from grid.simplex import contains + + +class TestSequenceFunctions(unittest.TestCase): + def setUp(self): + self.points = [ + [ 0.25, 0.40], # 0 + [ 0.60, 0.80], # 1 + [ 0.65, 0.28], # 2 + [ 0.28, 0.65], # 3 + [ 1.00, 0.75], # 4 + [ 0.30, 0.95], # 5 + [ 0.80, 0.50], # 6 + [ 0.35, 0.15], # 7 + ] + self.q = [exact_func(p) for p in self.points] + + self.X = [0.55, 0.45] + self.X = [0.25, 0.4001] + + self.g = grid(self.points, self.q) + self.g.construct_connectivity() + self.R = self.g.create_mesh(range(3)) + self.S = self.g.create_mesh(range(3,len(self.points))) + + self.exact = exact_func(self.X) + + self.answer = run_baker(self.X, self.R, self.S) + + self.accuracy = 8 + + def test_R_contains_X(self): + self.assertTrue(contains(self.X, self.R.points)) + + def test_RunBaker(self): + print + print "X\n", self.X + print "R\n", self.R + print "S\n", self.S + print "exact\n",self.exact + print "qlin\n",self.answer['qlin'] + self.assertTrue(self.answer) + +if __name__ == '__main__': + suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) + unittest.TextTestRunner(verbosity=3).run(suite) diff --git a/tools/blender/plot_cloud.py b/tools/blender/plot_cloud.py index 99ba17a..80ef195 100644 --- a/tools/blender/plot_cloud.py +++ b/tools/blender/plot_cloud.py @@ -1,9 +1,9 @@ from Blender import * import bpy -from grid.test3d import simple_rect_grid +from grid.test3d import rect_grid -g = simple_rect_grid(10, 10, 20) +g = rect_grid(10, 10, 20) me = bpy.data.meshes.new('points') me.verts.extend(g.points)