From 2eaec8268c45c343ca30d06a74fa86f86c2ebb82 Mon Sep 17 00:00:00 2001 From: Stephen Mardson McQuay Date: Sat, 23 Oct 2010 16:57:44 -0600 Subject: [PATCH] created a distince gmsh 2D plugin --- bin/parse_gmsh.py | 98 +++------------------------------------------ interp/grid/gmsh.py | 89 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 92 deletions(-) create mode 100644 interp/grid/gmsh.py diff --git a/bin/parse_gmsh.py b/bin/parse_gmsh.py index b0e011c..633cece 100755 --- a/bin/parse_gmsh.py +++ b/bin/parse_gmsh.py @@ -1,106 +1,20 @@ #!/usr/bin/env python -import atexit -import os -import sys -import readline -import rlcompleter - -historyPath = os.path.expanduser("~/.pyhistory") - -def save_history(historyPath=historyPath): - import readline - readline.write_history_file(historyPath) - -if os.path.exists(historyPath): - readline.read_history_file(historyPath) - -atexit.register(save_history) -del os, atexit, readline, rlcompleter, save_history, historyPath - import sys -import pickle -from itertools import combinations -from collections import defaultdict -from scipy.spatial import KDTree - -from grid.simplex import face import numpy as np -THREE_NODE_TRIANGLE = 2 +from interp.grid.gmsh import gmsh_grid if __name__ == '__main__': if len(sys.argv) != 2: print >> sys.stderr, "usage: %s " % sys.argv[0] sys.exit(1) - gmsh_file = open(sys.argv[1], 'r') + g = gmsh_grid(sys.argv[1]) + g.dump_to_blender_files() + X = np.array([0.2, 0.2, 0.0]) + R = g.get_containing_simplex(X) - gmsh_file.readline() # $MeshFormat - format = gmsh_file.readline() - gmsh_file.readline() # $EndMeshFormat - - gmsh_file.readline() # $Nodes - - node_count = int(gmsh_file.readline()) - - verts = np.zeros((node_count, 3)) - for i in xrange(node_count): - cur_line = gmsh_file.readline() - (index, x,y,z) = cur_line.split() - index = int(index) - 1 - x = float(x) - y = float(y) - z = float(z) - - verts[i][0] = x - verts[i][1] = y - verts[i][2] = z - - tree = KDTree(verts) - - gmsh_file.readline() # $E ndNodes - gmsh_file.readline() # $E lements - - element_count = int(gmsh_file.readline()) - - faces = {} - neighbors = {} - faces_for_vert = defaultdict(list) - - for i in xrange(element_count): - cur_line = gmsh_file.readline() - cur_line = cur_line.split() - cur_face_index, node_type, rest = (int(cur_line[0]), - int(cur_line[1]), - [int(j) for j in cur_line[2:]]) - - if(node_type == THREE_NODE_TRIANGLE): - points_for_cur_face = [i-1 for i in rest[rest[0]+1:]] - - cur_face = face(cur_face_index) - - for cur_point in points_for_cur_face: - faces_for_vert[cur_point].append(cur_face_index) - - cur_face.verts = points_for_cur_face - - faces[cur_face_index] = cur_face - edges = [tuple(sorted(i)) for i in combinations(points_for_cur_face, 2)] - - # edge is two verts - for edge in edges: - if edge in neighbors: - neighbors[edge].append(cur_face_index) - else: - neighbors[edge] = [cur_face_index] - - for k,v in neighbors.iteritems(): - if len(v) > 1: - faces[v[0]].add_neighbor(faces[v[1]]) - faces[v[1]].add_neighbor(faces[v[0]]) - - pickle.dump([(p[0], p[1], p[2]) for p in verts], open('/tmp/points.p', 'w')) - pickle.dump([f.verts for f in faces.itervalues()], open('/tmp/faces.p', 'w')) + print g diff --git a/interp/grid/gmsh.py b/interp/grid/gmsh.py new file mode 100644 index 0000000..d4cb2aa --- /dev/null +++ b/interp/grid/gmsh.py @@ -0,0 +1,89 @@ +import pickle + +from itertools import combinations +from collections import defaultdict + +import numpy as np +from scipy.spatial import KDTree + +from interp.grid import grid +from interp.grid.simplex import face +from interp.tools import exact_func + +class gmsh_grid(grid): + THREE_NODE_TRIANGLE = 2 + + def __init__(self, filename): + """ + construct an interp.grid.grid-compliant grid + object out of a 2D 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) + self.q[i] = exact_func(self.verts[i]) + + + grid.__init__(self) + gmsh_file.readline() # $EndNodes + gmsh_file.readline() # $Elements + + # temporary dict used to compute face 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_face_index, node_type, rest = (int(cur_line[0]), + int(cur_line[1]), + [int(j) for j in cur_line[2:]]) + + if(node_type == gmsh_grid.THREE_NODE_TRIANGLE): + points_for_cur_face = [i-1 for i in rest[rest[0]+1:]] + + cur_face = face(cur_face_index) + + for cur_point in points_for_cur_face: + self.faces_for_vert[cur_point].append(cur_face) + + cur_face.verts = points_for_cur_face + + self.faces[cur_face_index] = cur_face + edges = [tuple(sorted(i)) for i in combinations(points_for_cur_face, 2)] + + # edge is two verts + for edge in edges: + if edge in neighbors: + neighbors[edge].append(cur_face_index) + else: + neighbors[edge] = [cur_face_index] + + for k,v in neighbors.iteritems(): + if len(v) > 1: + self.faces[v[0]].add_neighbor(self.faces[v[1]]) + self.faces[v[1]].add_neighbor(self.faces[v[0]]) + + def dump_to_blender_files(self, pfile = '/tmp/points.p', ffile = '/tmp/faces.p'): + pickle.dump([(p[0], p[1], p[2]) for p in self.verts], open(pfile, 'w')) + pickle.dump([f.verts for f in self.faces.itervalues()], open(ffile, 'w'))