diff --git a/bin/driver.py b/bin/driver.py index 63fef5c..a6e8985 100755 --- a/bin/driver.py +++ b/bin/driver.py @@ -4,9 +4,7 @@ import sys from optparse import OptionParser import numpy as np -import scipy.spatial -import baker from tools import rms, exact_func import grid @@ -72,36 +70,21 @@ if __name__ == '__main__': help = "verbosity") (options, args) = parser.parse_args() - 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()) - print >> sys.stderr, "wrote source mesh output to %s" % options.output + if options.verbose: + print >> sys.stderr, "options: %s, args: %s" % (options, args) + print >> sys.stderr, "wrote source mesh output to %s" % options.output errors = [] success = 0 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.verbose) + answer = mesh_source.run_baker(X) if answer['a'] == None: errors.append(0) @@ -114,9 +97,9 @@ if __name__ == '__main__': if options.verbose: print "current point : %s" % X print "exact : %0.4f" % exact - print "qlin : %0.4f" % answer['lin'] + print "qlin : %0.4f" % answer['qlin'] print "q_final : %0.4f" % answer['final'] - print "qlinerr : %1.4f" % (exact - anser['lin'],) + print "qlinerr : %1.4f" % (exact - answer['qlin'],) print "q_final_err : %0.4f" % (exact - answer['final'],) cur_error = np.abs(answer['final'] - exact) errors.append(cur_error) diff --git a/lib/baker.py b/lib/baker.py index 7196fec..7595f99 100644 --- a/lib/baker.py +++ b/lib/baker.py @@ -1,4 +1,3 @@ -from grid import exact_func import numpy as np import sys @@ -93,16 +92,16 @@ def qlinear_3D(X, R, q): qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)]) return qlin -def run_baker(X, R, S, verbose = False): +def run_baker(X, R, S): """ - This is the main function to call to get an interpolation to X from the tree + 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] - g -- the grid object + R = Simplex - tree -- the kdtree search object (built from the g mesh) + S = extra points """ diff --git a/lib/grid.py b/lib/grid.py index 8450794..fbc4361 100755 --- a/lib/grid.py +++ b/lib/grid.py @@ -7,6 +7,7 @@ from collections import defaultdict import numpy as np import scipy.spatial +from baker import run_baker from tools import exact_func from smcqdelaunay import * @@ -58,9 +59,36 @@ class grid(object): 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 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, 3 + extra_points) + + + # get the containing simplex + R = [self.points[i] for i in indicies[:simplex_size] ] + Rq = [self.q[i] for i in indicies[:simplex_size] ] + r_mesh = grid(R, Rq) + + # and some extra points + S = [self.points[i] for i in indicies[simplex_size:] ] + Sq = [self.q[i] for i in indicies[simplex_size:] ] + s_mesh = grid(S, Sq) + + return (r_mesh, s_mesh) + + def run_baker(self, X): + (R, S) = self.get_simplex_and_nearest_points(X) + return run_baker(X, R, S) + def construct_connectivity(self): """ a call to this method prepares the internal connectivity structure.