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

This commit is contained in:
Stephen Mardson McQuay 2010-03-05 08:58:07 -07:00
parent a25e6d03d1
commit 9a1b8d14b2
4 changed files with 50 additions and 33 deletions

View File

@ -6,30 +6,38 @@ import pickle
from grid import simple_rect_grid, simple_random_grid from grid import simple_rect_grid, simple_random_grid
from baker import run_baker from baker import run_baker
from tools import smberror
qfile = '/tmp/grid_regular.txt' qfile = '/tmp/grid_regular.txt'
if __name__ == '__main__': if __name__ == '__main__':
try: try:
rx = int(sys.argv[1]) rx = int(sys.argv[1])
ry = int(sys.argv[2]) ry = int(sys.argv[2])
except ValueError as e: except (IndexError, ValueError) as e:
print e print "problem with argv: %s" % e
rx = 4 rx = 4
ry = 4 * rx ry = 4 * rx
source_mesh = simple_rect_grid(rx, ry) source_mesh = simple_rect_grid(rx, ry)
print source_mesh
# print source_mesh
X = [0.1, 0.1] X = [0.1, 0.1]
try:
(R, S) = source_mesh.get_simplex_and_nearest_points(X, extra_points=4) (R, S) = source_mesh.get_simplex_and_nearest_points(X, extra_points=4)
print "R for nearest-neighbor:\n", R print "R for nearest-neighbor:\n", R
print "S for nearest-neighbor:\n", S print "S for nearest-neighbor:\n", S
print run_baker(X, R, S) print run_baker(X, R, S)
except smberror as e:
print "caught error: %s" % e
(R, S) = source_mesh.get_points_conn(X) (R, S) = source_mesh.get_points_conn(X)
print "R for connectivity:\n", R print "R for connectivity:\n", R
print "S for connectivity:\n", S print "S for connectivity:\n", S
print run_baker(X, R, 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()) open(qfile, 'w').write(source_mesh.for_qhull())

View File

@ -1,6 +1,8 @@
import numpy as np import numpy as np
import sys import sys
from tools import smberror
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. 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: try:
phi = np.linalg.solve(A,b) phi = np.linalg.solve(A,b)
except: 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) phi = np.dot(np.linalg.pinv(A), b)
return phi return phi
@ -58,13 +61,13 @@ def get_phis_3D(X, r):
try: try:
phi = np.linalg.solve(A,b) phi = np.linalg.solve(A,b)
except: 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) phi = np.dot(np.linalg.pinv(A), b)
return phi return phi
def qlinear(X, R, q): def qlinear(X, R):
""" """
this calculates the linear portion of q from X to 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 q = CFD quantities of interest at the simplex points
""" """
phis = get_phis(X, R) phis = get_phis(X, R.points)
qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) qlin = sum([q_i * phi_i for q_i, phi_i in zip(R.q, phis)])
return qlin return phis, qlin
def qlinear_3D(X, R, q): def qlinear_3D(X, R, q):
""" """
@ -90,7 +93,7 @@ def qlinear_3D(X, R, q):
phis = get_phis_3D(X, R) phis = get_phis_3D(X, R)
qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) 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): def run_baker(X, R, S):
""" """
@ -106,8 +109,7 @@ def run_baker(X, R, S):
""" """
# calculate values only for the triangle # calculate values only for the triangle
phi = get_phis(X, R.points) phi, qlin = qlinear (X, R)
qlin = qlinear (X, R.points, R.q)
if len(S.points) == 0: if len(S.points) == 0:
answer = { answer = {
@ -124,10 +126,11 @@ def run_baker(X, R, S):
w = [] # baker eq 11 w = [] # baker eq 11
for (s, q) in zip(S.points, S.q): for (s, q) in zip(S.points, S.q):
(phi1, phi2, phi3) = get_phis(s, R.points) cur_phi, cur_qlin = qlinear(s, R)
B.append([phi1 * phi2, phi2*phi3, phi3*phi1]) (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) B = np.array(B)
w = np.array(w) w = np.array(w)
@ -139,7 +142,7 @@ def run_baker(X, R, S):
try: try:
(a, b, c) = np.linalg.solve(A,b) (a, b, c) = np.linalg.solve(A,b)
except: 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) (a, b, c) = np.dot(np.linalg.pinv(A), b)
error_term = a * phi[0] * phi[1]\ error_term = a * phi[0] * phi[1]\

View File

@ -8,7 +8,7 @@ import numpy as np
import scipy.spatial import scipy.spatial
from baker import run_baker, get_phis from baker import run_baker, get_phis
from tools import exact_func from tools import exact_func, smberror
from smcqdelaunay import * from smcqdelaunay import *
class face(object): class face(object):
@ -146,13 +146,20 @@ class grid(object):
return R, S return R, S
def run_baker(self, X, connectivity_based = False): def run_baker(self, X):
print >>sys.stderr, "suxor" answer = None
if connectivity_based:
(R, S) = self.get_points_conn(X) try:
else:
(R, S) = self.get_simplex_and_nearest_points(X) (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): def construct_connectivity(self):
""" """

View File

@ -30,7 +30,6 @@ class TestSequenceFunctions(unittest.TestCase):
import scipy import scipy
import grid import grid
import baker import baker
import delaunay
def testGetPhis(self): def testGetPhis(self):
@ -62,7 +61,7 @@ class TestSequenceFunctions(unittest.TestCase):
r = [[0, 0], [1, 0], [1, 1]] r = [[0, 0], [1, 0], [1, 1]]
q = [1, 0, 0] q = [1, 0, 0]
result = baker.qlinear(X, r, q) phi, result = baker.qlinear(X, grid.grid(r,q))
right_answer = 0.5 right_answer = 0.5