a3d72ed200
I'm going to hurry and clean up this method, but for now it is verbose, and uses sets. i'm going to just check it in for future reference
274 lines
7.0 KiB
Python
Executable File
274 lines
7.0 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.
|
|
|
|
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()
|
|
|
|
(dist, indicies) = self.tree.query(X, 2)
|
|
smplx = None
|
|
for i in self.facets_for_point[indicies[0]]:
|
|
if i.contains(X, self):
|
|
smplx = i
|
|
break
|
|
|
|
if not smplx:
|
|
raise AssertionError('no containing simplex found')
|
|
print "containing verts:", smplx.verts
|
|
smplx_set = set(smplx.verts)
|
|
R = self.create_mesh(smplx.verts)
|
|
|
|
|
|
s = []
|
|
for c,i in enumerate(smplx.neighbors):
|
|
print "neighbor %d: %s: %s" % (c, i.name, i.verts)
|
|
st = set(i.verts)
|
|
s.append((st - smplx_set).pop())
|
|
# s.append( [guy for guy in i.verts if not guy in smplx.verts])
|
|
print s
|
|
S = self.create_mesh(s)
|
|
|
|
print S
|
|
return R, S
|
|
|
|
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.
|
|
|
|
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()
|