diff --git a/bin/parse_gmsh.py b/bin/parse_gmsh.py index 3b6d2d3..1805d0d 100755 --- a/bin/parse_gmsh.py +++ b/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: - print a['final'] - e, a, e + # if not ia: + # print a['final'] - e, a, e results[ia] += 1 except Exception,e: print >>sys.stderr, e diff --git a/interp/grid/gmsh.py b/interp/grid/gmsh.py index 319d8de..90266e8 100644 --- a/interp/grid/gmsh.py +++ b/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) diff --git a/interp/tools.py b/interp/tools.py index 6feebd2..f8619ce 100644 --- a/interp/tools.py +++ b/interp/tools.py @@ -46,9 +46,8 @@ def exact_func(X): """ the exact function used from baker's article (for testing) """ - x = X[0] - y = X[1] - answer = np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) + x ,y = X + answer = 1.0 + x*x + y*y # np.power((np.sin(x * np.pi) * np.cos(y * np.pi)), 2) log.debug(answer) return answer