changed the 2D gmsh to work with new baker method

bin/parse_gmsh.py

@@ -14,9 +14,10 @@ if __name__ == '__main__':
     sys.exit(1)
 
   g = gmsh_grid(sys.argv[1])
+  g.q = np.array([exact_func(x) for x in g.verts])
   # g.dump_to_blender_files()
 
-  X = np.array([0.2, 0.5, 0.0])
+  X = np.array([0.2, 0.5]) #TODO:, 0.0])
   R = g.get_containing_simplex(X)
   print R
   R, S = g.get_simplex_and_nearest_points(X, 10)
@@ -64,18 +65,17 @@ if __name__ == '__main__':
   e = exact_func(X)
   print e
   print improved_answer(a, e)
-  sys.exit(0)
 
   results = {True:0, False:0}
   for i in xrange(1000):
-    X = np.random.random((1,3))[0]
+    X = np.random.random((1,2))[0]
 
     try:
-      a = g.run_baker(X, order=2, extra_points = 10)
+      a = g.run_baker(X, order=2, extra_points = 3)
       e = exact_func(X)
       ia = improved_answer(a, e)
-      if not ia:
+      # if not ia:
-        print a['final'] - e, a, e
+      #   print a['final'] - e, a, e
       results[ia] += 1
     except Exception,e:
       print >>sys.stderr, e

interp/grid/gmsh.py

@@ -39,7 +39,8 @@ class gmsh_grid(grid):
 
     node_count = int(gmsh_file.readline())
 
-    self.verts = np.empty((node_count, 3))
+    # for dim = 2, see note in next for loop
+    self.verts = np.empty((node_count, 2))
     self.q     = np.empty(node_count)
 
     for i in xrange(node_count):
@@ -49,8 +50,10 @@ class gmsh_grid(grid):
 
       self.verts[i][0] = float(x)
       self.verts[i][1] = float(y)
-      self.verts[i][2] = float(z)
+
-      # self.q[i] = exact_func(self.verts[i])
+      # 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)
 
 
     grid.__init__(self)

interp/tools.py

@@ -46,9 +46,8 @@ def exact_func(X):
   """
     the exact function used from baker's article (for testing)
   """
-  x = X[0]
+  x ,y = X
-  y = X[1]
+  answer = 1.0 + x*x + y*y # np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2)
-  answer = np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2)
   log.debug(answer)
   return answer