diff --git a/bin/driver.py b/bin/driver.py index 0bbc835..41fd3ae 100755 --- a/bin/driver.py +++ b/bin/driver.py @@ -7,24 +7,77 @@ import numpy as np import scipy.spatial import baker -from tools import rms, get_mesh, exact_func +from tools import rms, exact_func +import grid + + +def get_mesh(source, destination, use_structured_grid = False): + mesh_source = None + mesh_dest = None + if use_structured_grid: + mesh_source = grid.simple_rect_grid(source, source) + mesh_dest = grid.simple_rect_grid(destination, destination) + else: + mesh_source = grid.simple_random_grid(source) + mesh_dest = grid.simple_random_grid(destination) + + if not (mesh_dest and mesh_source): + raise smberror('problem creating mesh objects') + else: + return mesh_source, mesh_dest if __name__ == '__main__': parser = OptionParser() - parser.add_option("-o", "--output-file", dest="output", type='str', default = '/tmp/for_qhull.txt', help = "qhull output file") - parser.add_option("-e", "--extra-points", dest="extra", type='int', default = 3, help = "how many extra points") + parser.add_option("-o", + "--output-file", + dest="output", + type='str', + default = '/tmp/for_qhull.txt', + help = "qhull output file") - parser.add_option("-s", "--source-total", dest="source_total", type='int', default = 100, help = "total number of source points for random, resolution for structured") - parser.add_option("-d", "--destination-total", dest="destination_total", type='int', default = 100, help = "total number of destination points, resolution for structured") + parser.add_option("-e", + "--extra-points", + dest="extra", + type='int', + default = 3, + help = "how many extra points") - parser.add_option("-r", "--structured", action = 'store_true', default = False, help = "use a structured grid instead of random point cloud") - parser.add_option("-v", "--verbose", action = 'store_true', default = False, help = "verbosity") + parser.add_option("-s", + "--source-total", + dest="source_total", + type='int', + default = 100, + help = "total number of source points for random,\ + resolution for structured") + + parser.add_option("-d", + "--destination-total", + dest="destination_total", + type='int', + default = 100, + help = "total number of destination points,\ + resolution for structured") + + parser.add_option("-r", + "--structured", + action = 'store_true', + default = False, + help = "use a structured grid instead of random point cloud") + + parser.add_option("-v", + "--verbose", + action = 'store_true', + default = False, + help = "verbosity") (options, args) = parser.parse_args() - print >> sys.stderr, options + print >> sys.stderr, "options: %s, args: %s" % (options, args) - mesh_source, mesh_dest = get_mesh(options.source_total, options.destination_total, options.structured) + + + + (mesh_source, mesh_dest) = get_mesh(options.source_total, options.destination_total, options.structured) tree = scipy.spatial.KDTree(mesh_source.points) open(options.output, 'w').write(mesh_source.for_qhull()) @@ -32,20 +85,36 @@ if __name__ == '__main__': errors = [] success = 0 - for x in mesh_dest.points: - lin, error, final = baker.run_baker(x, mesh_source, tree, options.extra, options.verbose) - exact = exact_func(x[0], x[1]) - if np.abs(exact - final) < np.abs(exact - lin): + + for X in mesh_dest.points: + + (dist, indicies) = tree.query(X, 3 + options.extra) + + + # get the containing simplex + R = [mesh_source.points[i] for i in indicies[:3] ] + Rq = [mesh_source.q[i] for i in indicies[:3] ] + r_mesh = grid.grid(R, Rq) + + # and some extra points + S = [mesh_source.points[i] for i in indicies[3:] ] + Sq = [mesh_source.q[i] for i in indicies[3:] ] + s_mesh = grid.grid(S, Sq) + + answer = baker.run_baker(X, r_mesh, s_mesh, options.extra, options.verbose) + + exact = exact_func(X[0], X[1]) + if np.abs(exact - answer['final']) < np.abs(exact - answer['qlin']): success += 1 if options.verbose: - print "current point : %s" % x + print "current point : %s" % X print "exact : %0.4f" % exact - print "qlin : %0.4f" % lin - print "q_final : %0.4f" % final - print "qlinerr : %0.4f" % (exact - lin,) - print "q_final_err : %0.4f" % (exact - final,) - cur_error = np.abs(final - exact) + print "qlin : %0.4f" % answer['lin'] + print "q_final : %0.4f" % answer['final'] + print "qlinerr : %1.4f" % (exact - anser['lin'],) + print "q_final_err : %0.4f" % (exact - answer['final'],) + cur_error = np.abs(answer['final'] - exact) errors.append(cur_error) print rms(errors) diff --git a/lib/baker.py b/lib/baker.py index bd7697b..34f2cec 100644 --- a/lib/baker.py +++ b/lib/baker.py @@ -8,7 +8,7 @@ def get_phis(X, r): X -- the destination point (2D) X = [0,0] - r -- the three points that make up the triangle (2D) + 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] @@ -65,33 +65,35 @@ def get_phis_3D(X, r): return phi -def qlinear(X, r, q): +def qlinear(X, R, q): """ this calculates the linear portion of q from X to r + + also, this is baker eq 3 X = destination point - r = simplex points + R = simplex points q = CFD quantities of interest at the simplex points """ - phis = get_phis(X, r) - qlin = sum([q_i * phi_i for q_i, phi_i in zip(q[:len(phis)], phis)]) + phis = get_phis(X, R) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) return qlin -def qlinear_3D(X, r, q): +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) + 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[:len(phis)], phis)]) + phis = get_phis_3D(X, R) + qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) return qlin -def run_baker(X, g, tree, extra_points = 3, verbose = False): +def run_baker(X, R, S, extra_points = 3, verbose = False): """ This is the main function to call to get an interpolation to X from the tree @@ -104,44 +106,47 @@ def run_baker(X, g, tree, extra_points = 3, verbose = False): """ - (dist, indicies) = tree.query(X, 3 + extra_points) - - nn = [g.points[i] for i in indicies] - nq = [g.q[i] for i in indicies] - # calculate values only for the triangle - phi = get_phis(X, nn[:3]) - qlin = qlinear(X, nn[:3], nq[:3])# nq[0] * phi[0] + nq[1] * phi[1] + nq[2] * phi[2] + phi = get_phis(X, S.points) + qlin = qlinear (X, S.points, S.q) - error_term = 0.0 + if extra_points == 0: return qlin - if extra_points != 0: - B = [] # baker eq 9 - w = [] # baker eq 11 + B = [] # baker eq 9 + w = [] # baker eq 11 - for index in indicies[3:]: - (phi1,phi2,phi3) = get_phis(g.points[index], nn) - B.append([phi1 * phi2, phi2*phi3, phi3*phi1]) + for (s, q) in zip(S.points, S.q): + (phi1, phi2, phi3) = get_phis(s, R.points) + B.append([phi1 * phi2, phi2*phi3, phi3*phi1]) - w.append(g.q[index] - qlinear(g.points[index], nn, nq)) + w.append(q - qlinear(s, R.points, R.q)) - B = np.array(B) - w = np.array(w) + B = np.array(B) + w = np.array(w) - A = np.dot(B.T, B) - b = np.dot(B.T, 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: linear calculation went bad, resorting to np.linalg.pinv" - (a, b, c) = np.dot(np.linalg.pinv(A), b) + # baker solve eq 10 + try: + (a, b, c) = np.linalg.solve(A,b) + except: + print >> sys.stderr, "warning: 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] + error_term = a * phi[0] * phi[1]\ + + b * phi[1] * phi[2]\ + + c * phi[2] * phi[0] q_final = qlin + error_term - return qlin, error_term, q_final + answer = { + 'a': a, + 'b': b, + 'c': c, + 'qlin': qlin, + 'error': error_term, + 'final': q_final, + } + + return answer diff --git a/lib/grid.py b/lib/grid.py index 066d64a..6257fd9 100755 --- a/lib/grid.py +++ b/lib/grid.py @@ -71,5 +71,9 @@ class simple_random_grid(simple_rect_grid): if __name__ == '__main__': - g = simple_random_grid(100) + try: + resolution = int(sys.argv[1]) + except: + resolution = 10 + g = simple_rect_grid(resolution, resolution) print g.for_qhull() diff --git a/lib/tools.py b/lib/tools.py index 8a99e52..77bcb11 100644 --- a/lib/tools.py +++ b/lib/tools.py @@ -1,6 +1,4 @@ import numpy as np -import grid - class smberror(Exception): def __init__(self, val): @@ -21,19 +19,3 @@ def rms(errors): def exact_func(x, y): return np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) return np.sin(x * np.pi) * np.cos(y * np.pi) - - -def get_mesh(source, destination, use_structured_grid = False): - mesh_source = None - mesh_dest = None - if use_structured_grid: - mesh_source = grid.simple_rect_grid(source, source) - mesh_dest = grid.simple_rect_grid(destination, destination) - else: - mesh_source = grid.simple_random_grid(source) - mesh_dest = grid.simple_random_grid(destination) - - if not (mesh_dest and mesh_source): - raise smberror('problem creating mesh objects') - else: - return mesh_source, mesh_dest diff --git a/test/baker.test.py b/test/baker.test.py index ecdc90f..bd4c3d1 100755 --- a/test/baker.test.py +++ b/test/baker.test.py @@ -11,6 +11,19 @@ class TestSequenceFunctions(unittest.TestCase): def setUp(self): self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]] self.approx_fmt = "%0.6f" + self.all_points = [ + [ 0, 0], # 0 + [ 1, 0], # 1 + [ 1, 1], # 2 + [ 0, 1], # 3 + [ 1,-1], # 4 + [ 0,-1], # 5 + [-1, 1], # 6 + [-1, 0], # 7 + [-1,-1], # 8 + ] + self.q = [1, 0, 0, 0, 0, 0, 0, 0, 0] + self.X = [1.5, 10.25] def testGetPhis(self): @@ -48,29 +61,80 @@ class TestSequenceFunctions(unittest.TestCase): self.assertEqual(result, right_answer) - def testRunBaker(self): - X = [0.5, 0.25] + def testRunBaker_1(self): + size_of_simplex = 3 + extra_points = 3 - all_points = [ - [ 0, 0], # 0 - [ 1, 0], # 1 - [ 1, 1], # 2 - [ 0, 1], # 3 - [ 1,-1], # 4 - [ 0,-1], # 5 - [-1, 1], # 6 - [-1, 0], # 7 - [-1,-1], # 8 - ] - q = [1, 0, 0, 0, 0, 0, 0, 0, 0] - mesh = grid.grid(all_points, q) - tree = scipy.spatial.KDTree(all_points) - (final, exact) = baker.run_baker(X, mesh, tree) - print final, exact - result = 3 - right_answer = 3 + R = grid.grid(self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) - self.assertEqual(result, right_answer) + S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + + answer = baker.run_baker(self.X, R, S) + a = self.approx_fmt % answer['a'] + b = self.approx_fmt % answer['b'] + c = self.approx_fmt % answer['c'] + + self.assertEqual(a, c) + self.assertEqual(c, self.approx_fmt % 0.00) + self.assertEqual(b, self.approx_fmt % (1/3.0)) + + def testRunBaker_2(self): + size_of_simplex = 3 + extra_points = 4 + + R = grid.grid(self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) + + S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + + answer = baker.run_baker(self.X, R, S) + a = self.approx_fmt % answer['a'] + b = self.approx_fmt % answer['b'] + c = self.approx_fmt % answer['c'] + + self.assertEqual(a, c) + self.assertEqual(c, self.approx_fmt % float(2/3.0)) + + def testRunBaker_3(self): + size_of_simplex = 3 + extra_points = 5 + + R = grid.grid(self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) + + S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + + answer = baker.run_baker(self.X, R, S) + a = self.approx_fmt % answer['a'] + b = self.approx_fmt % answer['b'] + c = self.approx_fmt % answer['c'] + + self.assertEqual(a, self.approx_fmt % float(13/14.0)) + self.assertEqual(b, self.approx_fmt % float(2 / 7.0)) + self.assertEqual(c, self.approx_fmt % float(15/14.0)) + + def testRunBaker_4(self): + size_of_simplex = 3 + extra_points = 6 + + R = grid.grid(self.all_points[:size_of_simplex], + self.q[:size_of_simplex]) + + S = grid.grid(self.all_points[size_of_simplex:size_of_simplex + extra_points], + self.q[size_of_simplex:size_of_simplex + extra_points]) + + answer = baker.run_baker(self.X, R, S) + a = self.approx_fmt % answer['a'] + b = self.approx_fmt % answer['b'] + c = self.approx_fmt % answer['c'] + + self.assertEqual(a, self.approx_fmt % float(48/53.0)) + self.assertEqual(b, self.approx_fmt % float(15/53.0)) + self.assertEqual(c, self.approx_fmt % float(54/53.0)) if __name__ == '__main__': suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) diff --git a/test/baker3d.test.py b/test/baker3d.test.py deleted file mode 100755 index ec8b6b6..0000000 --- a/test/baker3d.test.py +++ /dev/null @@ -1,76 +0,0 @@ -#!/usr/bin/python - -import unittest -import baker - -import numpy as np -import scipy.spatial - -class TestSequenceFunctions(unittest.TestCase): - def setUp(self): - self.l = [[-1, 1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]] - self.r = [[1,2,3], [2,2,3], [1,3,3], [1,2,9]] - self.approx_fmt = "%0.6f" - - def testGetPhis(self): - - X = [0,0] - r = [[-1, -1], [0, 2], [1, -1]] - - result = baker.get_phis(X, r) - result = [self.approx_fmt % i for i in result] - - right_answer = [self.approx_fmt % i for i in [1/3.0, 1/3.0, 1/3.0]] - - for a,b in zip(result, right_answer): - self.assertEqual(a,b) - - def testGetPhis2(self): - - X = [0.5,0.25] - r = [[0, 0], [1, 0], [1, 1]] - - result = baker.get_phis(X, r) - - right_answer = [0.5, 0.25, 0.25] - - for a,b in zip(result, right_answer): - self.assertEqual(a,b) - - def testQlinear(self): - X = [0.5, 0.25] - r = [[0, 0], [1, 0], [1, 1]] - q = [1, 0, 0] - - result = baker.qlinear(X, r, q) - - right_answer = 0.5 - - self.assertEqual(result, right_answer) - - def testRunBaker(self): - X = [0.5, 0.25] - - all_points = [ - [ 0, 0], # 0 - [ 1, 0], # 1 - [ 1, 1], # 2 - [ 0, 1], # 3 - [ 1,-1], # 4 - [ 0,-1], # 5 - [-1, 1], # 6 - [-1, 0], # 7 - [-1,-1], # 8 - ] - q = [1, 0, 0, 0, 0, 0, 0, 0, 0] - - tree = scipy.spatial.KDTree(all_points) - baker.run_baker(X, - result = 3 - right_answer = 5 - - self.assertEqual(result, right_answer) - -if __name__ == '__main__': - suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) - unittest.TextTestRunner(verbosity=2).run(suite)