a25e6d03d1
I spoke with Dr. Gorrell about first trying to call the nearest-neighbor routine, releasing the trapped exception where the matrix ends up being singular, and trapping it in the grid object's run_baker, where I could then call the connectivity-based lookup routine. Those are my next steps, as well as starting to think about pushing the 3D front ahead this week.
281 lines
7.3 KiB
Python
Executable File
281 lines
7.3 KiB
Python
Executable File
#!/usr/bin/python
|
|
|
|
import sys
|
|
import re
|
|
from collections import defaultdict
|
|
|
|
import numpy as np
|
|
import scipy.spatial
|
|
|
|
from baker import run_baker, get_phis
|
|
from tools import exact_func
|
|
from smcqdelaunay import *
|
|
|
|
class face(object):
|
|
def __init__(self, name):
|
|
self.name = name
|
|
self.verts = []
|
|
self.neighbors = []
|
|
|
|
def add_vert(self, v):
|
|
"""
|
|
v should be an index into grid.points
|
|
"""
|
|
self.verts.append(v)
|
|
|
|
def add_neighbor(self, n):
|
|
"""
|
|
reference to another face object
|
|
"""
|
|
self.neighbors.append(n)
|
|
|
|
def contains(self, X, grid):
|
|
R = [grid.points[i] for i in self.verts]
|
|
|
|
phis = get_phis(X, R)
|
|
|
|
r = True
|
|
if [i for i in phis if i < 0.0]:
|
|
r = False
|
|
return r
|
|
|
|
def __str__(self):
|
|
neighbors = [i.name for i in self.neighbors]
|
|
return '%s: verts: %s neighbors: [%s]' %\
|
|
(
|
|
self.name,
|
|
self.verts,
|
|
", ".join(neighbors)
|
|
)
|
|
|
|
|
|
|
|
|
|
|
|
class grid(object):
|
|
facet_re = re.compile(r'''
|
|
-\s+(?P<facet>f\d+).*?
|
|
vertices:\s(?P<verts>.*?)\n.*?
|
|
neighboring\s facets:\s+(?P<neigh>[\sf\d]*)
|
|
''', re.S|re.X)
|
|
|
|
point_re = re.compile(r'''
|
|
-\s+(?P<point>p\d+).*?
|
|
neighbors:\s+(?P<neigh>[\sf\d]*)
|
|
''', re.S|re.X)
|
|
|
|
vert_re = re.compile(r'''
|
|
(p\d+)
|
|
''', re.S|re.X)
|
|
|
|
|
|
def __init__(self, points, q):
|
|
"""
|
|
this thing eats two pre-constructed arrays of stuff:
|
|
points = array of arrays (i will convert to numpy.array)
|
|
[[x0,y0], [x1,y1], ...]
|
|
q = array (1D) of important values
|
|
"""
|
|
|
|
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 create_mesh(self, indicies):
|
|
p = [self.points[i] for i in indicies]
|
|
q = [self.q[i] for i in indicies]
|
|
return grid(p, q)
|
|
|
|
|
|
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_mesh = self.create_mesh(indicies[:simplex_size])
|
|
# and some extra points
|
|
s_mesh = self.create_mesh(indicies[simplex_size:])
|
|
|
|
return (r_mesh, s_mesh)
|
|
|
|
def get_points_conn(self, X):
|
|
"""
|
|
this returns two grid objects: R and S.
|
|
|
|
this function differes from the get_simplex_and_nearest_points
|
|
function in that it builds up the extra points based on
|
|
connectivity information, not just nearest-neighbor.
|
|
in theory, this will work much better for situations like
|
|
points near a short edge in a boundary layer cell where the
|
|
nearest points would all be colinear
|
|
|
|
R is a grid object that is the (a) containing simplex around point X
|
|
S is a connectivity-based nearest-neighbor lookup, limited to 3 extra points
|
|
"""
|
|
if not self.faces:
|
|
self.construct_connectivity()
|
|
|
|
# get closest point
|
|
(dist, indicies) = self.tree.query(X, 2)
|
|
|
|
simplex = None
|
|
for i in self.facets_for_point[indicies[0]]:
|
|
if i.contains(X, self):
|
|
simplex = i
|
|
break
|
|
|
|
if not simplex:
|
|
raise AssertionError('no containing simplex found')
|
|
|
|
R = self.create_mesh(simplex.verts)
|
|
|
|
|
|
s = []
|
|
for c,i in enumerate(simplex.neighbors):
|
|
s.extend([guy for guy in i.verts if not guy in simplex.verts])
|
|
S = self.create_mesh(s)
|
|
|
|
return R, S
|
|
|
|
def run_baker(self, X, connectivity_based = False):
|
|
print >>sys.stderr, "suxor"
|
|
if connectivity_based:
|
|
(R, S) = self.get_points_conn(X)
|
|
else:
|
|
(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.
|
|
|
|
this is part of the __init__ for a simple_rect_grid, but can be called from any grid object
|
|
"""
|
|
qdelaunay_string = get_qdelaunay_dump_str(self)
|
|
facet_to_facets = []
|
|
for matcher in grid.facet_re.finditer(qdelaunay_string):
|
|
d = matcher.groupdict()
|
|
|
|
facet_name = d['facet']
|
|
verticies = d['verts']
|
|
neighboring_facets = d['neigh']
|
|
|
|
cur_face = face(facet_name)
|
|
self.faces[facet_name] = cur_face
|
|
|
|
for v in grid.vert_re.findall(verticies):
|
|
vertex_index = int(v[1:])
|
|
cur_face.add_vert(vertex_index)
|
|
self.facets_for_point[vertex_index].append(cur_face)
|
|
|
|
nghbrs = [(facet_name, i) for i in neighboring_facets.split()]
|
|
facet_to_facets.extend(nghbrs)
|
|
|
|
for rel in facet_to_facets:
|
|
if rel[1] in self.faces:
|
|
self.faces[rel[0]].add_neighbor(self.faces[rel[1]])
|
|
|
|
# for matcher in grid.point_re.finditer(qdelaunay_string):
|
|
# d = matcher.groupdict()
|
|
|
|
# point = d['point']
|
|
# neighboring_facets = d['neigh']
|
|
|
|
# self.facets_for_point[int(point[1:])] = [i for i in neighboring_facets.split() if i in self.faces]
|
|
|
|
def for_qhull_generator(self):
|
|
"""
|
|
this returns a generator that should be fed into qdelaunay
|
|
"""
|
|
|
|
yield '2';
|
|
yield '%d' % len(self.points)
|
|
|
|
for p in self.points:
|
|
yield "%f %f" % (p[0], p[1])
|
|
|
|
def for_qhull(self):
|
|
"""
|
|
this returns a single string that should be fed into qdelaunay
|
|
"""
|
|
r = '2\n'
|
|
r += '%d\n' % len(self.points)
|
|
for p in self.points:
|
|
r += "%f %f\n" % (p[0], p[1])
|
|
return r
|
|
|
|
def __str__(self):
|
|
r = ''
|
|
assert( len(self.points) == len(self.q) )
|
|
for c, i in enumerate(zip(self.points, self.q)):
|
|
r += "%d %r: %0.4f" % (c,i[0], i[1])
|
|
facet_str = ", ".join([f.name for f in self.facets_for_point[c]])
|
|
r += " faces: [%s]" % facet_str
|
|
r += "\n"
|
|
if self.faces:
|
|
for v in self.faces.itervalues():
|
|
r += "%s\n" % v
|
|
return r
|
|
|
|
class simple_rect_grid(grid):
|
|
def __init__(self, xres = 5, yres = 5):
|
|
xmin = -1.0
|
|
xmax = 1.0
|
|
xspan = xmax - xmin
|
|
xdel = xspan / float(xres - 1)
|
|
|
|
ymin = -1.0
|
|
ymay = 1.0
|
|
yspan = ymay - ymin
|
|
ydel = yspan / float(yres - 1)
|
|
|
|
|
|
points = []
|
|
q = []
|
|
for x in xrange(xres):
|
|
cur_x = xmin + (x * xdel)
|
|
for y in xrange(yres):
|
|
cur_y = ymin + (y * ydel)
|
|
points.append([cur_x, cur_y])
|
|
q.append(exact_func(cur_x, cur_y))
|
|
grid.__init__(self, points, q)
|
|
self.construct_connectivity()
|
|
|
|
|
|
|
|
class simple_random_grid(simple_rect_grid):
|
|
def __init__(self, num_points = 10):
|
|
points = []
|
|
q = []
|
|
|
|
r = np.random
|
|
|
|
for i in xrange(num_points):
|
|
cur_x = r.rand()
|
|
cur_y = r.rand()
|
|
|
|
points.append([cur_x, cur_y])
|
|
q.append(exact_func(cur_x, cur_y))
|
|
grid.__init__(self, points, q)
|
|
|
|
self.points = np.array(self.points)
|
|
self.q = np.array(self.q)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
try:
|
|
resolution = int(sys.argv[1])
|
|
except:
|
|
resolution = 10
|
|
g = simple_rect_grid(resolution, resolution)
|
|
print g.for_qhull()
|