smbinterp/interp/grid/gmsh.py

102 lines
2.6 KiB
Python
Raw Normal View History

2010-10-23 15:57:44 -07:00
from itertools import combinations
import numpy as np
from scipy.spatial import KDTree
from interp.grid import grid
2011-03-22 06:15:34 -07:00
from interp.grid import cell
2011-03-22 15:06:08 -07:00
import logging
log = logging.getLogger('interp')
THREE_NODE_TRIANGLE = 2
FOUR_NODE_TET = 4
EDGES_FOR_FACE_CONNECTIVITY = 2
EDGES_FOR_VOLUME_CONNECTIVITY = 3
2010-10-23 15:57:44 -07:00
class ggrid(grid):
2010-10-23 15:57:44 -07:00
def __init__(self, filename, dimension = 3):
2010-10-23 15:57:44 -07:00
"""
construct an interp.grid.grid-compliant grid
object out of a {2,3}D gmsh file
2010-10-23 15:57:44 -07:00
"""
self.dim = dimension
log.debug("dimension: %d", self.dim)
2010-10-23 15:57:44 -07:00
gmsh_file = open(filename, 'r')
gmsh_file.readline() # $MeshFormat
gmsh_file.readline()
2010-10-23 15:57:44 -07:00
gmsh_file.readline() # $EndMeshFormat
gmsh_file.readline() # $Nodes
node_count = int(gmsh_file.readline())
self.verts = np.empty((node_count, dimension))
2010-10-23 15:57:44 -07:00
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)
if self.dim == 3:
self.verts[i][2] = float(z)
2010-10-23 15:57:44 -07:00
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
# 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 == 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)
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, self.dim)]
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]])