moved teh grid code into the grid class

by next friday i need to have implemented a way to select between normal kdtree point lookups, and nearest-neighbor lookups. Hopefully i'll have that done tonight.
This commit is contained in:
smcquay@cfdviz2 2010-02-20 13:18:24 -07:00
parent 2cf9a0b574
commit 6c4c6d4cfe
3 changed files with 39 additions and 29 deletions

View File

@ -4,9 +4,7 @@ import sys
from optparse import OptionParser from optparse import OptionParser
import numpy as np import numpy as np
import scipy.spatial
import baker
from tools import rms, exact_func from tools import rms, exact_func
import grid import grid
@ -72,36 +70,21 @@ if __name__ == '__main__':
help = "verbosity") help = "verbosity")
(options, args) = parser.parse_args() (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()) 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 = [] errors = []
success = 0 success = 0
for X in mesh_dest.points: for X in mesh_dest.points:
(dist, indicies) = tree.query(X, 3 + options.extra) answer = mesh_source.run_baker(X)
# 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)
if answer['a'] == None: if answer['a'] == None:
errors.append(0) errors.append(0)
@ -114,9 +97,9 @@ if __name__ == '__main__':
if options.verbose: if options.verbose:
print "current point : %s" % X print "current point : %s" % X
print "exact : %0.4f" % exact 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 "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'],) print "q_final_err : %0.4f" % (exact - answer['final'],)
cur_error = np.abs(answer['final'] - exact) cur_error = np.abs(answer['final'] - exact)
errors.append(cur_error) errors.append(cur_error)

View File

@ -1,4 +1,3 @@
from grid import exact_func
import numpy as np import numpy as np
import sys 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)]) qlin = sum([q_i * phi_i for q_i, phi_i in zip(q, phis)])
return qlin 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 -- the destination point (2D)
X = [0,0] X = [0,0]
g -- the grid object R = Simplex
tree -- the kdtree search object (built from the g mesh) S = extra points
""" """

View File

@ -7,6 +7,7 @@ from collections import defaultdict
import numpy as np import numpy as np
import scipy.spatial import scipy.spatial
from baker import run_baker
from tools import exact_func from tools import exact_func
from smcqdelaunay import * from smcqdelaunay import *
@ -58,9 +59,36 @@ class grid(object):
self.points = np.array(points) self.points = np.array(points)
self.q = np.array(q) self.q = np.array(q)
self.tree = scipy.spatial.KDTree(self.points)
self.faces = {} self.faces = {}
self.facets_for_point = defaultdict(list) 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): def construct_connectivity(self):
""" """
a call to this method prepares the internal connectivity structure. a call to this method prepares the internal connectivity structure.