diff --git a/bin/dgrid_bench.py b/bin/dgrid_bench.py index ca7df5f..0fba098 100644 --- a/bin/dgrid_bench.py +++ b/bin/dgrid_bench.py @@ -4,13 +4,21 @@ import pickle import numpy as np +from optparse import OptionParser + import interp.bootstrap from interp.grid.delaunay import dgrid from interp.tools import improved_answer, friendly_exact_3D as exact -def test_success(rcount, count): - v = np.random.random((rcount, 3)) * 10 +if __name__ == '__main__': + if len(sys.argv) != 3: + print >> sys.stderr, "usage: %s " % sys.argv[0] + sys.exit(1) + + rcount, count = (int(i) for i in sys.argv[1:]) + + v = np.random.random((rcount, 3)) pickle.dump(v, open('/tmp/v.p', 'w')) q = np.array([exact(x) for x in v]) @@ -34,15 +42,8 @@ def test_success(rcount, count): results[improved_answer(a, e)] += 1 i += 1 except Exception as e: + print e outside += 1 print "total skipped points: %d" % outside - return results - -if __name__ == '__main__': - if len(sys.argv) != 3: - print >> sys.stderr, "usage: %s " % sys.argv[0] - sys.exit(1) - - rcount, count = (int(i) for i in sys.argv[1:]) - print test_success(rcount, int(count)) + print results diff --git a/bin/gmsh_bench.py b/bin/gmsh_bench.py index 7a47bac..8e3ec5c 100644 --- a/bin/gmsh_bench.py +++ b/bin/gmsh_bench.py @@ -7,9 +7,9 @@ import numpy as np from optparse import OptionParser import interp.bootstrap -from interp.grid.gmsh import gmsh_grid3D, gmsh_grid -from interp.tools import * - +from interp.grid.gmsh import ggrid +from interp.grid.delaunay import dgrid +from interp.tools import improved_answer, baker_exact_2D, baker_exact_3D, friendly_exact_2D, friendly_exact_3D def get_right_exact_func(options): f = None @@ -25,10 +25,9 @@ def get_right_exact_func(options): f = friendly_exact_2D return f - - if __name__ == '__main__': - parser = OptionParser(usage = "%prog [options] ") + + parser = OptionParser(usage = "%prog [options] ") parser.add_option("-v", "--verbose", action="store_true", dest="verbose", default=False, @@ -38,12 +37,19 @@ if __name__ == '__main__': action="store_true", dest="baker", default=False, help="use more impressive baker funcs") + meshtype_options = ('gmsh', 'delaunay') + parser.add_option('-t', '--type', + type="str", dest="meshtype", default='dgrid', + help="specify mesh type (default: %default, options gmsh, delaunay)") + parser.add_option("-d", "--dimension", type="int", dest="dimension", default=3, help="use 3D or 2D gmsh lib (default: %default)") + parser.add_option("-o", "--order", type="int", dest="order", default=3, help="order of interpolation (default: %default)") + parser.add_option("-e", "--extra-points", type="int", dest="extra", default=3, help="number of extra points (default: %default)") @@ -51,26 +57,38 @@ if __name__ == '__main__': (options, args) = parser.parse_args() if len(args) != 2: - # print >> sys.stderr, "usage: %s " % sys.argv[0] parser.print_usage() sys.exit(1) input_file, count = args - if options.dimension == 3: - g = gmsh_grid3D(input_file) + if options.meshtype == 'gmsh': + g = ggrid(input_file, options.dimension) + elif options.meshtype == 'delaunay': + point_count = int(input_file) + v = np.random.random((point_count, options.dimension)) + g = dgrid(v) else: - g = gmsh_grid(input_file) + raise Exception("improper mesh type specified. options: %s\n%s" % (str(meshtype_options), options.meshtype)) g.q = np.array([get_right_exact_func(options)(x) for x in g.verts]) results = {True:0, False:0} - for i in xrange(int(count)): - X = np.random.random((1,options.dimension))[0] + i = 0 + outside = 0 + while i < int(count): + try: + X = np.random.random((1,options.dimension))[0] - a = g.run_baker(X, order = options.order, extra_points = options.extra) - e = get_right_exact_func(options)(X) + a = g.run_baker(X, order = options.order, extra_points = options.extra) + e = get_right_exact_func(options)(X) - results[improved_answer(a, e)] += 1 + results[improved_answer(a, e)] += 1 + i += 1 + except Exception as e: + print e + print X, i, count + outside += 1 + print "total skipped points: %d" % outside print results diff --git a/interp/grid/__init__.py b/interp/grid/__init__.py index 19807ed..43fee4f 100644 --- a/interp/grid/__init__.py +++ b/interp/grid/__init__.py @@ -77,6 +77,7 @@ class grid(object): log.debug("simplex vert indicies: %s" % simplex.verts) R = self.create_mesh(simplex.verts) + log.debug("R:\n%s", R) log.debug('total attempts before finding simplex: %d' % attempts) return R @@ -92,7 +93,7 @@ class grid(object): 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): + def get_simplex_and_nearest_points(self, X, extra_points = 3): """ this returns two grid objects: R and S. @@ -101,6 +102,7 @@ class grid(object): S is S_j from baker's paper : some verts from all point that are not the simplex """ + simplex_size = self.dim + 1 log.debug("extra verts: %d" % extra_points) log.debug("simplex size: %d" % simplex_size) @@ -234,12 +236,13 @@ class cell(object): return interp.grid.simplex.contains(X, [G.verts[i] for i in self.verts]) def __str__(self): - neighbors = [str(i.name) for i in self.neighbors] - return '' %\ + # neighbors = [str(i.name) for i in self.neighbors] + return '' %\ ( self.name, self.verts, - ", ".join(neighbors) + len(self.neighbors), + # ", ".join(neighbors) ) __repr__ = __str__ diff --git a/interp/grid/delaunay.py b/interp/grid/delaunay.py index 74d76cc..339a352 100644 --- a/interp/grid/delaunay.py +++ b/interp/grid/delaunay.py @@ -44,6 +44,7 @@ class dgrid(basegrid): ''', re.S|re.X) def __init__(self, verts, q = None): + self.dim = len(verts[0]) basegrid.__init__(self, verts,q) self.construct_connectivity() diff --git a/interp/grid/gmsh.py b/interp/grid/gmsh.py index 19ddcdd..58c8f78 100644 --- a/interp/grid/gmsh.py +++ b/interp/grid/gmsh.py @@ -22,27 +22,28 @@ EDGES_FOR_VOLUME_CONNECTIVITY = 3 -class gmsh_grid(grid): +class ggrid(grid): - def __init__(self, filename): + def __init__(self, filename, dimension = 3): """ construct an interp.grid.grid-compliant grid - object out of a 2D gmsh file + object out of a {2,3}D gmsh file """ + self.dim = dimension + log.debug("dimension: %d", self.dim) gmsh_file = open(filename, 'r') gmsh_file.readline() # $MeshFormat - format = gmsh_file.readline() + fmat = gmsh_file.readline() gmsh_file.readline() # $EndMeshFormat gmsh_file.readline() # $Nodes node_count = int(gmsh_file.readline()) - # for dim = 2, see note in next for loop - self.verts = np.empty((node_count, 2)) + self.verts = np.empty((node_count, dimension)) self.q = np.empty(node_count) for i in xrange(node_count): @@ -53,13 +54,16 @@ class gmsh_grid(grid): self.verts[i][0] = float(x) self.verts[i][1] = float(y) - # for the general baker method to work, it must have 2 - # components if it it a 2D mesh, so I removed: - # self.verts[i][2] = float(z) + if self.dim == 3: + self.verts[i][2] = float(z) - grid.__init__(self) self.tree = KDTree(self.verts) + + # initialize rest of structures about to be populated (cells, + # cells_for_vert) + grid.__init__(self) + gmsh_file.readline() # $EndNodes gmsh_file.readline() # $Elements @@ -74,7 +78,8 @@ class gmsh_grid(grid): int(cur_line[1]), [int(j) for j in cur_line[2:]]) - if(node_type == THREE_NODE_TRIANGLE): + if (node_type == THREE_NODE_TRIANGLE and self.dim == 2) \ + or (node_type == FOUR_NODE_TET and self.dim == 3): points_for_cur_cell = [i-1 for i in rest[rest[0]+1:]] cur_cell = cell(cur_cell_index) @@ -85,81 +90,7 @@ class gmsh_grid(grid): cur_cell.verts = points_for_cur_cell self.cells[cur_cell_index] = cur_cell - edges = [tuple(sorted(i)) for i in combinations(points_for_cur_cell, EDGES_FOR_FACE_CONNECTIVITY)] - - # edge is two verts - for edge in edges: - if edge in neighbors: - neighbors[edge].append(cur_cell_index) - else: - neighbors[edge] = [cur_cell_index] - - for k,v in neighbors.iteritems(): - if len(v) > 1: - self.cells[v[0]].add_neighbor(self.cells[v[1]]) - self.cells[v[1]].add_neighbor(self.cells[v[0]]) - - -class gmsh_grid3D(grid): - - def __init__(self, filename): - """ - construct an interp.grid.grid-compliant grid - object out of a 3D gmsh file - """ - - gmsh_file = open(filename, 'r') - - - gmsh_file.readline() # $MeshFormat - format = gmsh_file.readline() - gmsh_file.readline() # $EndMeshFormat - - gmsh_file.readline() # $Nodes - - node_count = int(gmsh_file.readline()) - - self.verts = np.empty((node_count, 3)) - self.q = np.empty(node_count) - - for i in xrange(node_count): - cur_line = gmsh_file.readline() - (index, x,y,z) = cur_line.split() - index = int(index) - 1 - - self.verts[i][0] = float(x) - self.verts[i][1] = float(y) - self.verts[i][2] = float(z) - - - grid.__init__(self) - self.tree = KDTree(self.verts) - gmsh_file.readline() # $EndNodes - gmsh_file.readline() # $Elements - - # temporary dict used to compute cell connectivity - neighbors = {} - - element_count = int(gmsh_file.readline()) - for i in xrange(element_count): - cur_line = gmsh_file.readline() - cur_line = cur_line.split() - cur_cell_index, node_type, rest = (int(cur_line[0]), - int(cur_line[1]), - [int(j) for j in cur_line[2:]]) - - if(node_type == FOUR_NODE_TET): - points_for_cur_cell = [i-1 for i in rest[rest[0]+1:]] - - cur_cell = cell(cur_cell_index) - - for cur_point in points_for_cur_cell: - self.cells_for_vert[cur_point].append(cur_cell) - - cur_cell.verts = points_for_cur_cell - - self.cells[cur_cell_index] = cur_cell - edges = [tuple(sorted(i)) for i in combinations(points_for_cur_cell, EDGES_FOR_VOLUME_CONNECTIVITY)] + edges = [tuple(sorted(i)) for i in combinations(points_for_cur_cell, self.dim)] for edge in edges: if edge in neighbors: