merged the gmsh meshes (2/3D) into a single mesh, updated other files to support this

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):
+if __name__ == '__main__':
-  v = np.random.random((rcount, 3)) * 10
+  if len(sys.argv) != 3:
+    print >> sys.stderr, "usage: %s <number of random points> <number of attempted interps>" % 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
+  print results
-
-if __name__ == '__main__':
-  if len(sys.argv) != 3:
-    print >> sys.stderr, "usage: %s <number of random points> <number of attempted interps>" % sys.argv[0]
-    sys.exit(1)
-
-  rcount, count = (int(i) for i in sys.argv[1:])
-  print test_success(rcount, int(count))

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.grid.gmsh     import ggrid
-from   interp.tools     import *
+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] <input file> <number of attempts>")
+
+  parser = OptionParser(usage = "%prog [options] <input file (gmsh) or random grid count (delaunay)> <number of attempts>")
 
   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 <input file> <number of attempts>" % sys.argv[0]
     parser.print_usage()
     sys.exit(1)
 
   input_file, count = args
 
-  if options.dimension == 3:
+  if options.meshtype == 'gmsh':
-    g = gmsh_grid3D(input_file)
+    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)):
+  i = 0
-    X = np.random.random((1,options.dimension))[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)
+      a = g.run_baker(X, order = options.order, extra_points = options.extra)
-    e = get_right_exact_func(options)(X)
+      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

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]
+    # neighbors = [str(i.name) for i in self.neighbors]
-    return '<cell %s: verts: %s neighbors: [%s]>' %\
+    return '<cell %s: verts: %s neighbor count: %s>' %\
            (
              self.name,
              self.verts,
-             ", ".join(neighbors)
+             len(self.neighbors),
+             # ", ".join(neighbors)
            )
 
   __repr__ = __str__

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()
 

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, dimension))
-    self.verts = np.empty((node_count, 2))
     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
+      if self.dim == 3:
-      # components if it it a 2D mesh, so I removed:
+        self.verts[i][2] = float(z)
-      # 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)]
+        edges = [tuple(sorted(i)) for i in combinations(points_for_cur_cell, self.dim)]
-
-        # 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)]
 
         for edge in edges:
           if edge in neighbors: