#!/usr/bin/env python
#
# Read a gmsh (http://www.geuz.org/gmsh/) mesh file (.msh extension
# typical but we do not check), assumed to supply an all-quad mesh,
# and from this, produce a semtex session file skeleton, containing
# NODES and ELEMENTS sections. Wrap this with other standard semtex
# feml file sections in order to generate a valid semtex session file.
#
# Print all this to standard output.
#
# Usage: gmsh2sem.py gmshfile.msh
# Usage: python gmsh2sem.py gmshfile.msh
#
# NB: use Mesh.SubdivisionAlgorithm=1; in gmsh in order to generate quads.
# Standard semtex utility meshpr must be found in your shell's PATH variable.
#
# $Id: gmsh2sem.py,v 6.2 2010/05/11 11:47:01 hmb Exp $
##############################################################################
import sys, os
try:
infilename = sys.argv[1]
except:
print "Usage: gmesh2sem infile.msh"; sys.exit(1)
ifile = open (infilename, 'r')
t1 = open (".t1", 'w')
# -- We will assume that the file is a valid gmsh mesh file with quad meshes.
# First print a generic semtex header.
s = """# Skeleton semtex session file generated from gmsh .msh file
N_P = 5
N_STEP = 20
D_T = 0.01
KINVIS = 1.
u v p
1 w wall
1 w 3
u = 0.
v = 0.
p
\n
"""
t1.write(s)
# -- Extract NODES and ELEMENTS from given gmsh file.
while 1:
line = ifile.readline()
if '$Nodes' in line:
line = ifile.readline()
Nnode = int (line)
s = '\n'
t1.write(s)
for i in range (1, Nnode+1):
line = ifile.readline()
t1.write(line),
t1.write('\n\n')
elif '$Elements' in line:
line = ifile.readline()
Nel = int (line)
s = '\n'
t1.write(s)
for i in range (1, Nel+1):
# -- If all elements were quads, each line has 10 numbers.
# The first one is a tag, and the last 4 are node numbers.
line = ifile.readline()
words = line.split()
s = words[0]+' '+words[6]+' '+words[7]+' '+words[8]+' '+words[9]+'
\n'
t1.write(s)
t1.write('\n\n')
break
ifile.close()
# -- Generate a SURFACES section that contains all the unfilled
# surface valencies (those element edges that do not mate to
# another element). You can then edit this (by hand, or other
# means).
t1.close()
s = "meshpr -s .t1 > .t2"
os.system(s)
t1 = open (".t1", 'r')
t2 = open (".t2", 'r')
# -- Print up all of t1 and t2
while 1:
line = t1.readline()
if not line: break
print line,
while 1:
line = t2.readline()
if not line: break
print line,
t1.close(); t2.close()
os.remove('.t1'); os.remove('.t2')