From 9a1b8d14b24b42d594788b34ca07193b6f87c8eb Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Fri, 5 Mar 2010 08:58:07 -0700 Subject: [PATCH] finally added proper (in my opinon) exception handling to the grid's run_baker method. it should properly autocorrect, and switch the extra point lookup scheme to connectivity-based if solving the system of equations yields a singular solution --- bin/grid_regular.py | 30 +++++++++++++++++++----------- lib/baker.py | 29 ++++++++++++++++------------- lib/grid.py | 21 ++++++++++++++------- test/baker.test.py | 3 +-- 4 files changed, 50 insertions(+), 33 deletions(-) diff --git a/bin/grid_regular.py b/bin/grid_regular.py index d91dc90..545e483 100755 --- a/bin/grid_regular.py +++ b/bin/grid_regular.py @@ -6,30 +6,38 @@ import pickle from grid import simple_rect_grid, simple_random_grid from baker import run_baker +from tools import smberror + qfile = '/tmp/grid_regular.txt' if __name__ == '__main__': try: rx = int(sys.argv[1]) ry = int(sys.argv[2]) - except ValueError as e: - print e + except (IndexError, ValueError) as e: + print "problem with argv: %s" % e rx = 4 ry = 4 * rx source_mesh = simple_rect_grid(rx, ry) - print source_mesh + + # print source_mesh X = [0.1, 0.1] - (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 run_baker(X, R, S) + try: + (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 run_baker(X, R, S) + except smberror as e: + print "caught error: %s" % e + (R, S) = source_mesh.get_points_conn(X) + print "R for connectivity:\n", R + print "S for connectivity:\n", S + print run_baker(X, R, S) - (R, S) = source_mesh.get_points_conn(X) - print "R for connectivity:\n", R - print "S for connectivity:\n", S - print run_baker(X, R, S) + print "repeating the above just using the grid object:" + print source_mesh.run_baker(X) open(qfile, 'w').write(source_mesh.for_qhull()) diff --git a/lib/baker.py b/lib/baker.py index 89ec8e4..d2289fb 100644 --- a/lib/baker.py +++ b/lib/baker.py @@ -1,6 +1,8 @@ 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. @@ -26,7 +28,8 @@ def get_phis(X, R): try: phi = np.linalg.solve(A,b) except: - print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system" + 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 @@ -58,13 +61,13 @@ def get_phis_3D(X, r): try: phi = np.linalg.solve(A,b) except: - print >> sys.stderr, "warning: calculation of phis yielded a linearly dependant system" + 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, q): +def qlinear(X, R): """ this calculates the linear portion of q from X to r @@ -75,9 +78,9 @@ def qlinear(X, R, q): 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, phis)]) - return qlin + 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): """ @@ -90,7 +93,7 @@ def qlinear_3D(X, R, q): phis = get_phis_3D(X, R) qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) - return qlin + return phis, qlin def run_baker(X, R, S): """ @@ -106,8 +109,7 @@ def run_baker(X, R, S): """ # calculate values only for the triangle - phi = get_phis(X, R.points) - qlin = qlinear (X, R.points, R.q) + phi, qlin = qlinear (X, R) if len(S.points) == 0: answer = { @@ -124,10 +126,11 @@ def run_baker(X, R, S): w = [] # baker eq 11 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]) + cur_phi, cur_qlin = qlinear(s, R) + (phi1, phi2, phi3) = cur_phi - w.append(q - qlinear(s, R.points, R.q)) + B.append([phi1 * phi2, phi2 * phi3, phi3 * phi1]) + w.append(q - cur_qlin) B = np.array(B) w = np.array(w) @@ -139,7 +142,7 @@ def run_baker(X, R, S): try: (a, b, c) = np.linalg.solve(A,b) except: - print >> sys.stderr, "warning: linear calculation went bad, resorting to np.linalg.pinv" + 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]\ diff --git a/lib/grid.py b/lib/grid.py index 389e3ad..950245f 100755 --- a/lib/grid.py +++ b/lib/grid.py @@ -8,7 +8,7 @@ import numpy as np import scipy.spatial from baker import run_baker, get_phis -from tools import exact_func +from tools import exact_func, smberror from smcqdelaunay import * class face(object): @@ -146,13 +146,20 @@ class grid(object): return R, S - def run_baker(self, X, connectivity_based = False): - print >>sys.stderr, "suxor" - if connectivity_based: - (R, S) = self.get_points_conn(X) - else: + def run_baker(self, X): + answer = None + + try: (R, S) = self.get_simplex_and_nearest_points(X) - return run_baker(X, R, S) + 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): """ diff --git a/test/baker.test.py b/test/baker.test.py index 91936aa..95cff5b 100755 --- a/test/baker.test.py +++ b/test/baker.test.py @@ -30,7 +30,6 @@ class TestSequenceFunctions(unittest.TestCase): import scipy import grid import baker - import delaunay def testGetPhis(self): @@ -62,7 +61,7 @@ class TestSequenceFunctions(unittest.TestCase): r = [[0, 0], [1, 0], [1, 1]] q = [1, 0, 0] - result = baker.qlinear(X, r, q) + phi, result = baker.qlinear(X, grid.grid(r,q)) right_answer = 0.5