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

working on the containing simplex problem. the get containing simplex function works, if the point is in the domain. I am trying to force the domain to contain the point (placing points along perimiter of the mesh). I am not doing that right, but have to run

Browse Source 
--HG--
rename : bin/qhull-029.txt => data/qhull-029.txt
rename : bin/qhull-029.txt.gv => data/qhull-029.txt.gv
This commit is contained in:
Stephen Mardson McQuay
2010-04-24 18:26:44 -06:00
parent cfcd1e15a2
commit 3dcc10ea0e
10 changed files with 604 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/baker/__init__.py
+17 -32
View File
@@ -1,4 +1,5 @@
from baker import *
from baker.tools import logging as log
import numpy as np
import sys
@@ -30,7 +31,7 @@ def get_phis(X, R):
phi = np.linalg.solve(A,b)
except np.linalg.LinAlgError as e:
msg = "warning: get_phis: calculation of phis yielded a linearly dependant system (%s)" % e
# TODO: log this -- > print >> sys.stderr, msg
log.error(msg)
raise smberror(msg)
phi = np.dot(np.linalg.pinv(A), b)
@@ -136,7 +137,7 @@ def get_error_quadratic(phi, R, S):
+ b * phi[1] * phi[2]\
+ c * phi[2] * phi[0]
return error_term, a, b, c
return error_term
def get_error_cubic(phi, R, S):
B = [] # baker eq 9
@@ -181,7 +182,7 @@ def get_error_cubic(phi, R, S):
+ f * phi[2] * phi[1] * phi[1]\
+ g * phi[0] * phi[1] * phi[2]\
return error_term, a, b, c
return error_term
def run_baker(X, R, S, order=2):
@@ -195,37 +196,29 @@ def run_baker(X, R, S, order=2):
S = extra points
"""
answer = {
'qlin': None,
'error': None,
'final': None,
}
# calculate values only for the simplex triangle
phi, qlin = qlinear(X, R)
if len(S.points) == 0:
answer = {
'a': None,
'b': None,
'c': None,
'qlin': qlin,
'error': None,
'final': None,
}
if order == 1:
answer['qlin'] = qlin
return answer
if order == 2:
error_term, a, b, c = get_error_quadratic(phi, R, S)
elif order == 2:
error_term = get_error_quadratic(phi, R, S)
elif order == 3:
error_term, a, b, c = get_error_cubic(phi, R, S)
error_term = get_error_cubic(phi, R, S)
else:
raise smberror('unacceptable order for baker method')
q_final = qlin + error_term
answer = {
'a': a,
'b': b,
'c': c,
'qlin': qlin,
'error': error_term,
'final': q_final,
}
answer['qlin' ] = qlin
answer['error'] = error_term
answer['final'] = q_final
return answer
@@ -243,14 +236,6 @@ def run_baker_3D(X, R, S):
# calculate values only for the triangle
phi, qlin = qlinear_3D(X, R)
if [i for i in phi if i <= 0.0]:
s = "this is not a containing simplex:\n"
s += " X: %s\n" % X
s += " R: %s\n" % R
s += " phi: %s, sum(%0.4e)\n" % (phi, sum(phi))
print >> sys.stderr, s
raise smberror("not containing simplex")
if len(S.points) == 0:
answer = {
'a': None,
lib/baker/tools.py
+10
View File
@@ -1,5 +1,15 @@
import os
import logging
logging.basicConfig(
level = logging.DEBUG,
format = '%(asctime)s %(levelname)s %(message)s',
filename = os.path.join(os.sep, 'tmp', 'baker.lol'),
)
import numpy as np
class smberror(Exception):
"""
this is a silly little exception subclass
lib/grid/DD.py
+9
View File
@@ -61,6 +61,15 @@ class random_grid(rect_grid):
q = []
r = np.random
appx_side_res = int(np.sqrt(num_points))
delta = 1.0 / float(appx_side_res)
for x in xrange(appx_side_res):
cur_x = x * delta
for y in xrange(appx_side_res):
cur_y = y * delta
new_point = [cur_x, cur_y]
points.append(new_point)
q.append(exact_func(new_point))
for i in xrange(num_points):
cur_x = r.rand()
lib/grid/__init__.py
+45 -7
View File
@@ -6,7 +6,7 @@ import numpy as np
import scipy.spatial
from baker import run_baker
from baker.tools import exact_func, smberror
from baker.tools import exact_func, smberror, logging
from simplex import face, contains
from smcqdelaunay import *
@@ -42,7 +42,7 @@ class grid(object):
self.tree = scipy.spatial.KDTree(self.points)
self.faces = {}
self.facets_for_point = defaultdict(list)
self.counter = 0
def create_mesh(self, indicies):
p = [self.points[i] for i in indicies]
@@ -50,7 +50,45 @@ class grid(object):
return grid(p, q)
def get_containing_simplex(self, X):
pass
if not self.faces:
logging.debug('get_containing_simplex: setting up connectivity')
self.construct_connectivity()
# get closest point
(dist, indicies) = self.tree.query(X, 2)
closest_point = indicies[0]
logging.debug('counter: %d' % self.counter)
self.counter += 1
logging.debug('X: %s' % X)
logging.debug('point index: %d' % closest_point)
logging.debug('actual point %s' % self.points[closest_point])
logging.debug('distance = %0.4f' % dist[0])
simplex = None
checked_facets = []
facets_to_check = self.facets_for_point[closest_point]
attempts = 0
while not simplex:
logging.debug('attempt: %d' % attempts)
attempts += 1
if attempts > 10:
raise smberror("probably recursing to many times")
cur_facet = facets_to_check.pop()
checked_facets.append(cur_facet)
facets_to_check.extend([i for i in cur_facet.neighbors if i not in checked_facets])
if cur_facet.contains(X, self):
simplex = cur_facet
if not simplex:
raise AssertionError('no containing simplex found')
R = self.create_mesh(simplex.verts)
logging.debug('this must be R: %s' % R)
return R
def get_simplex_and_nearest_points(self, X, extra_points = 3, simplex_size = 3):
@@ -64,9 +102,8 @@ class grid(object):
"""
(dist, indicies) = self.tree.query(X, simplex_size + extra_points)
# get the containing simplex
r_mesh = self.create_mesh(indicies[:simplex_size])
# r_mesh = self.create_mesh(indicies[:simplex_size])
r_mesh = self.get_containing_simplex(X)
# and some extra points
s_mesh = self.create_mesh(indicies[simplex_size:])
@@ -103,7 +140,7 @@ class grid(object):
if not simplex:
raise AssertionError('no containing simplex found')
R = self.create_mesh(simplex.verts)
R = get_containing_simplex(X)# self.create_mesh(simplex.verts)
s = []
@@ -136,6 +173,7 @@ class grid(object):
this is part of the __init__ for a rect_grid, but can be called from any grid object
"""
logging.debug('calling construct_connectivity')
qdelaunay_string = get_qdelaunay_dump_str(self)
facet_to_facets = []
for matcher in grid.facet_re.finditer(qdelaunay_string):