smbinterp/bin/parse_gmsh.py

123 lines
3.0 KiB
Python
Executable File

#!/usr/bin/env python
import sys
import pickle
from pudb import set_trace
from itertools import combinations
THREE_NODE_TRIANGLE = 2
class point(object):
def __init__(self, index, x, y, z):
self.index = index
self.x = x
self.y = y
self.z = z
self.points = []
self.faces = []
def __getitem__(self, i):
if isinstance(i, int):
if i == 0:
return self.x
elif i == 1:
return self.y
elif i == 2:
return self.z
else:
raise IndexError("there are only (x,y,z)")
elif isinstance(i, str):
if i == 'x':
return self.x
elif i == 'y':
return self.y
elif i == 'z':
return self.z
else:
raise IndexError
else:
raise TypeError("only p.x, p[0], and p['x'] access allowed")
def __str__(self):
return '(%f, %f, %f)' % (self.x, self.y, self.z)
__repr__ = __str__
class face(object):
def __init__(self, index):
self.index = index
self.neighbors = []
if __name__ == '__main__':
if len(sys.argv) != 2:
print >> sys.stderr, "usage: %s <gmsh file>" % sys.argv[0]
sys.exit(1)
gmsh_file = open(sys.argv[1], 'r')
gmsh_file.readline() # $MeshFormat
format = gmsh_file.readline()
gmsh_file.readline() # $EndMeshFormat
gmsh_file.readline() # $Nodes
node_count = int(gmsh_file.readline())
points = []
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)
points.append(point(index, x,y,z))
gmsh_file.readline() # $EndNodes
gmsh_file.readline() # $Elements
element_count = int(gmsh_file.readline())
faces = {}
faceobjs = {}
faces_for_point = {}
neighbors = {}
for i in xrange(element_count):
cur_line = gmsh_file.readline()
cur_line = cur_line.split()
cur_face_index, type, rest = (int(cur_line[0]),
int(cur_line[1]),
[int(j) for j in cur_line[2:]])
if(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:
points[cur_point].faces.append(cur_face_index)
faces[cur_face_index] = points_for_cur_face
cur_face.points = points_for_cur_face
faceobjs[cur_face_index] = cur_face
edges = [tuple(sorted(i)) for i in combinations(points_for_cur_face, 2)]
# edge is two points
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:
faceobjs[v[0]].neighbors.append(v[1])
faceobjs[v[1]].neighbors.append(v[0])
pickle.dump([(p[0], p[1], p[2]) for p in points], open('/tmp/points.p', 'w'))
pickle.dump(faces , open('/tmp/faces.p', 'w'))