sm/smbinterp
1
0
Fork 0
Code Issues Pull Requests Activity

MAJOR: updated the baker method. it's more generic, and should allow me to 3D-ifiy it more simply. a ton of other things

Browse Source
This commit is contained in:
Stephen Mardson McQuay 2010-01-31 20:47:03 -07:00
parent 98b13fb8c5
commit 80720c45fe
6 changed files with 223 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
bin/driver.py
View File

@ -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)

85
lib/baker.py
View File

@ -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

6
lib/grid.py
View File

@ -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()

18
lib/tools.py
View File

@ -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

106
test/baker.test.py
View File

@ -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)

76
test/baker3d.test.py
View File

@ -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)