gmsh2adcirc.py

#!/usr/bin/env python3
#
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#                                                                       #
#                                 gmsh2adcirc.py                        # 
#                                                                       #
#+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!
#
# Author: Pat Prodanovic, Ph.D., P.Eng.
# 
# Date: June 26, 2015
#
# Modified: Feb 21, 2016
# Made it work under python 2 or 3
#
# Modified: May 30, 2017
# Added a check that makes sure elements are orinted in a CCW fashion
# before being written to the ADCIRC format.
#
# Purpose: Script takes the file generated by gmsh mesh generator, and 
# converts it to an ADCIRC mesh format
#
# Uses: Python 2 or 3, Matplotlib v1.4.2, Numpy v1.8.2
#
# Example:
#
# python gmsh2adcirc.py -i out.msh -o out.grd
# where:
# -i input *.msh file generated by gmsh
# -o output *.grd adcirc geometry mesh file
# 
# note: boundary nodes are not written
# 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Global Imports
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import os,sys
import numpy as np

# this works for python 2 and 3
def CCW(x1,y1,x2,y2,x3,y3):
  return (y3-y1)*(x2-x1) > (y2-y1)*(x3-x1)	

curdir = os.getcwd()
#
# I/O
if len(sys.argv) != 5 :
  print('Wrong number of Arguments, stopping now...')
  print('Usage:')
  print('python gmsh2adcirc.py -i out.msh -o out.grd')
  sys.exit()
dummy1 =  sys.argv[1]
gmsh_file = sys.argv[2]
dummy2 =  sys.argv[3]
adcirc_file = sys.argv[4]

# to create the output file
fout = open(adcirc_file,"w")

# to read the entire *.msh file as string (and clip of the stuff that is
# not needed
# each line in the file is a list object
line = list()
with open(gmsh_file, "r") as f1:
  for i in f1:
    line.append(i)
    
# to write to a temp file nodes and 2d elements only
temp_nodes_file ="temp_nodes"
temp_elements_file ="temp_elements"

fout_nodes = open(temp_nodes_file,"w")
fout_elements = open(temp_elements_file,"w")

# read how many nodes in the file (it is in *.gmsh file as line 5)
n = line[4] # it is a string
n = int(n)  # now it becomes an integer

# read the next n lines (this is the nodes file)
for i in range(5,5+n):
  fout_nodes.write(line[i])
  
# the next two lines are 
# $EndNodes 
# $Elements

# then number of elements
e = line[5+n+2]
e = int(e)

for i in range(5+n+3,5+n+3+e):
  # write only the 2d elements (i.e., if there are 7 blanks in the line)
  if(line[i].count(' ') == 7):
    fout_elements.write(line[i])
    
fout_nodes.close()
fout_elements.close()

# now open the nodes and elements files using numpy arrays

# use numpy to read the file
# each column in the file is a row in data read by no.loadtxt method
nodes_data = np.genfromtxt(temp_nodes_file,unpack=True)
elements_data = np.genfromtxt(temp_elements_file,unpack=True)

# nodes 
node_id = nodes_data[0,:]
node_id = node_id.astype(np.int32)
x = nodes_data[1,:]
y = nodes_data[2,:]
z = nodes_data[3,:]

# elements
e1 = elements_data[5,:]
e1 = e1.astype(np.int32)
e2 = elements_data[6,:]
e2 = e2.astype(np.int32)
e3 = elements_data[7,:]
e3 = e3.astype(np.int32)

# now to write the adcirc mesh file
fout.write("ADCIRC" + "\n")
# writes the number of elements and number of nodes in the header file
fout.write(str(len(e1)) + " " + str(len(node_id)) + "\n")

# added on 2017.05.30
# #######################
# make sure the elements are oriented in CCW fashion

ikle = np.column_stack((e1,e2,e3))

# go through each element, and make sure it is oriented in CCW fashion

for i in range(len(ikle)):
  
  # if the element is not CCW then must change its orientation
  if not CCW( x[ikle[i,0]-1], y[ikle[i,0]-1], x[ikle[i,1]-1], y[ikle[i,1]-1], 
    x[ikle[i,2]-1], y[ikle[i,2]-1] ):
    
    t0 = ikle[i,0]
    t1 = ikle[i,1]
    t2 = ikle[i,2]
    
    # switch orientation
    ikle[i,0] = t2
    ikle[i,2] = t0
# #######################

# writes the nodes
for i in range(0,len(node_id)):
  fout.write(str(node_id[i]) + " " + str("{:.3f}".format(x[i])) + " " + 
    str("{:.3f}".format(y[i])) + " " + str("{:.3f}".format(z[i])) + "\n")

# writes the elements
for i in range(0,len(e1)):
  fout.write(str(i+1) + " 3 " + str(ikle[i,0]) + " " + str(ikle[i,1]) + " " + 
    str(ikle[i,2]) + "\n")  

# now we can delete the temp file
os.remove(temp_nodes_file)
os.remove(temp_elements_file)