Update To include RMA10 - 3D - and mid side node
Mathieu Deiber
5 years ago
parent
93705a7a49
commit
ff8f9ff299
@ -0,0 +1,373 @@
|
||||
|
||||
# coding: utf-8
|
||||
|
||||
# In[64]:
|
||||
|
||||
import struct
|
||||
import matplotlib.pyplot as plt
|
||||
import math
|
||||
from py_rmatools import rma
|
||||
|
||||
plt.rcParams.update({'figure.max_open_warning': 0})
|
||||
|
||||
|
||||
# In[65]:
|
||||
|
||||
meshFilename = 'HWQ010.rm1'
|
||||
channelWidth = 100
|
||||
RMAfilename = 'HWQ277_2010_WQ'
|
||||
|
||||
#If RMA11
|
||||
constNum = [1,2,3,5,6,7,8,10,11,13,14,15,18]
|
||||
|
||||
|
||||
# In[66]:
|
||||
|
||||
def isElementOneD(nodelist):
|
||||
if len(nodelist) == 2:
|
||||
return True
|
||||
return False
|
||||
|
||||
def isElementSquare(nodelist):
|
||||
if len(nodelist) == 4:
|
||||
return True
|
||||
return False
|
||||
|
||||
def square2Triangle(ElementNum):
|
||||
nodelist = ElementDict[ElementNum]
|
||||
if isElementSquare(nodelist):
|
||||
ElementDict[ElementNum] = [nodelist[0], nodelist[1], nodelist[2]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist[0], nodelist[2], nodelist[3]]
|
||||
|
||||
def oneD2triangle(ElementNum):
|
||||
if isElementOneD(ElementDict[ElementNum]):
|
||||
nAe = ElementDict[ElementNum][0] #nAe Node A existing
|
||||
nBe = ElementDict[ElementNum][1]
|
||||
|
||||
if not nAe in node1Dduplicate: node1Dduplicate[nAe] = []
|
||||
if not nBe in node1Dduplicate: node1Dduplicate[nBe] = []
|
||||
|
||||
xA = nodeDict[nAe][0]
|
||||
xB = nodeDict[nBe][0]
|
||||
yA = nodeDict[nAe][1]
|
||||
yB = nodeDict[nBe][1]
|
||||
|
||||
normalVec = [-(yB - yA),(xB - xA)]
|
||||
dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2)
|
||||
normalVec[0] = normalVec[0] / dist
|
||||
normalVec[1] = normalVec[1] / dist
|
||||
xA2 = xA + channelWidth * normalVec[0]
|
||||
xB2 = xB + channelWidth * normalVec[0]
|
||||
yA2 = yA + channelWidth * normalVec[1]
|
||||
yB2 = yB + channelWidth * normalVec[1]
|
||||
|
||||
|
||||
nA = max(NodeList) + 1
|
||||
nB = max(NodeList) + 2
|
||||
|
||||
node1Dduplicate[nAe].append(nA)
|
||||
node1Dduplicate[nBe].append(nB)
|
||||
|
||||
node2nodevalue[nA] = nAe
|
||||
node2nodevalue[nB] = nBe
|
||||
|
||||
|
||||
NodeList.append(nA)
|
||||
NodeList.append(nB)
|
||||
nodeDict[nA] = [xA2, yA2, -1.01]
|
||||
nodeDict[nB] = [xB2, yB2, -1.01]
|
||||
|
||||
newEle = max(ElementList) + 1
|
||||
ElementList .append(newEle)
|
||||
ElementDict[ElementNum] = [nAe, nA, nBe]
|
||||
ElementDict[newEle] = [nA, nB, nBe]
|
||||
|
||||
def RMA11toSerafin():
|
||||
f = open('{}.slf'.format(RMAfilename), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
str='{0: >80}'.format('SERAFIN ')
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(constNum)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for c in constName:
|
||||
f.write(struct.pack(">l",32))
|
||||
str='{0: <32}'.format(c)
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",R.year))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
while R.next():
|
||||
#for i in range(3):
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",R.time * 3600))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
for c in constNum:
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.constit[c][node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.constit[c][key]))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
#R.next()
|
||||
|
||||
f.close()
|
||||
|
||||
|
||||
|
||||
def RMA2toSerafin():
|
||||
f = open('{}.slf'.format(RMAfilename), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
str='{0: >80}'.format('SERAFIN ')
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(constName)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for c in constName:
|
||||
f.write(struct.pack(">l",32))
|
||||
str='{0: <32}'.format(c)
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",R.year))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
while R.next():
|
||||
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",R.time * 3600))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.xvel[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.xvel[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.yvel[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.yvel[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.depth[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.depth[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.elevation[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.elevation[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
R.next()
|
||||
|
||||
f.close()
|
||||
|
||||
|
||||
# In[67]:
|
||||
|
||||
#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation
|
||||
NodeList = []
|
||||
ElementList = []
|
||||
ElementDict = {}
|
||||
nodeDict = {}
|
||||
node1Dduplicate = {} #Original Number: List of Duplicates
|
||||
node2nodevalue = {} #link between the node number and the node value to use
|
||||
#(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1},
|
||||
#but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated
|
||||
#node will be the same as the original so we might have {2050: 10})
|
||||
|
||||
with open(meshFilename) as f:
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
|
||||
cpt = 1
|
||||
while line and line != ' 9999\n':
|
||||
temp = line.split()
|
||||
ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9,2) if temp[i] != '0' and int(temp[9]) < 100]
|
||||
ElementList.append(int(temp[0]))
|
||||
line = f.readline()
|
||||
|
||||
for key, value in ElementDict.items():
|
||||
NodeList.extend(value)
|
||||
|
||||
NodeList = list(set(NodeList))
|
||||
|
||||
line = f.readline()
|
||||
while line and line != ' 9999\n':
|
||||
temp = line.split()
|
||||
if int(temp[0]) in NodeList:
|
||||
nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])]
|
||||
line = f.readline()
|
||||
|
||||
|
||||
for e in ElementList:
|
||||
oneD2triangle(e)
|
||||
square2Triangle(e)
|
||||
|
||||
for key in list(ElementDict): #Remove Special Element 902.....
|
||||
if len(ElementDict[key]) != 3:
|
||||
print(key, ElementDict[key])
|
||||
ElementDict.pop(key)
|
||||
ElementList.remove(key)
|
||||
|
||||
nodeOrdered = {}
|
||||
cpt = 1
|
||||
for key, value in nodeDict.items():
|
||||
nodeOrdered[key] = cpt
|
||||
cpt +=1
|
||||
|
||||
|
||||
|
||||
# # Open and Read First Step of the RMA File and Save a Serafin
|
||||
|
||||
# In[72]:
|
||||
|
||||
R=rma()
|
||||
R.open(RMAfilename)
|
||||
R.next()
|
||||
|
||||
if R.type==b'RMA11 ':
|
||||
constName = []
|
||||
for c in constNum:
|
||||
constName.append(R.constit_name[c].decode("utf-8"))
|
||||
RMA11toSerafin()
|
||||
|
||||
if R.type==b'RMA2 ':
|
||||
constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE']
|
||||
RMA2toSerafin()
|
||||
|
||||
@ -0,0 +1,738 @@
|
||||
|
||||
# coding: utf-8
|
||||
|
||||
# In[64]:
|
||||
|
||||
import struct
|
||||
import matplotlib.pyplot as plt
|
||||
import math
|
||||
from py_rmatools_v03 import rma
|
||||
import re
|
||||
from datetime import datetime, timedelta
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
plt.rcParams.update({'figure.max_open_warning': 0})
|
||||
|
||||
|
||||
# In[65]:
|
||||
|
||||
meshFilename = 'pok005.rm1'
|
||||
channelWidth = 100
|
||||
RMAfilenames = ['POA034']
|
||||
|
||||
#If RMA11
|
||||
constNum = [1]
|
||||
|
||||
#If RMA11 3D
|
||||
|
||||
#Options = [4,8,1,'Max'] #Option available 'Max' , an integer correspond to the layer number to extract the results, if 2D use options = [1] - 'Percentile' (percentile only work for RMA11 FOR now)
|
||||
Options = [8]
|
||||
percentiles = [0.5,0.8,0.9,0.95,1]
|
||||
NodeLayersFile = 'NodeLayers005.txt' #the node layer file was created by getting the nodes number from the RMA Outfile (just run RMA11 for a few steps)
|
||||
Nnodes = 6674 #Number of nodes (including mid side node) - thisn umber also correspond to the last node number in the nodelayers file
|
||||
# In[66]:
|
||||
|
||||
def isElementOneD(nodelist1):
|
||||
if len(nodelist1) == 3:
|
||||
return True
|
||||
return False
|
||||
|
||||
def isElementSquare(nodelist1):
|
||||
if len(nodelist1) == 8:
|
||||
return True
|
||||
return False
|
||||
|
||||
def isElementTriangle(nodelist1):
|
||||
if len(nodelist1) == 6:
|
||||
return True
|
||||
return False
|
||||
|
||||
def square2Triangle(ElementNum):
|
||||
nodelist1 = ElementDict[ElementNum]
|
||||
if isElementSquare(nodelist1):
|
||||
ElementDict[ElementNum] = [nodelist1[0], nodelist1[1], nodelist1[7]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[2], nodelist1[3]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[3], nodelist1[7]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[4], nodelist1[5]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[5], nodelist1[7]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[5], nodelist1[6], nodelist1[7]]
|
||||
|
||||
def triangle2Triangles(ElementNum):
|
||||
nodelist1 = ElementDict[ElementNum]
|
||||
if isElementTriangle(nodelist1):
|
||||
ElementDict[ElementNum] = [nodelist1[0], nodelist1[1], nodelist1[5]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[2], nodelist1[3]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[3], nodelist1[5]]
|
||||
ElementList.append(max(ElementList) + 1)
|
||||
ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[4], nodelist1[5]]
|
||||
|
||||
def oneD2triangle(ElementNum):
|
||||
if isElementOneD(ElementDict[ElementNum]):
|
||||
nAe = ElementDict[ElementNum][0] #nAe Node A existing
|
||||
nBe = ElementDict[ElementNum][1]
|
||||
|
||||
if not nAe in node1Dduplicate: node1Dduplicate[nAe] = []
|
||||
if not nBe in node1Dduplicate: node1Dduplicate[nBe] = []
|
||||
|
||||
xA = nodeDict[nAe][0]
|
||||
xB = nodeDict[nBe][0]
|
||||
yA = nodeDict[nAe][1]
|
||||
yB = nodeDict[nBe][1]
|
||||
|
||||
normalVec = [-(yB - yA),(xB - xA)]
|
||||
dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2)
|
||||
normalVec[0] = normalVec[0] / dist
|
||||
normalVec[1] = normalVec[1] / dist
|
||||
xA2 = xA + channelWidth * normalVec[0]
|
||||
xB2 = xB + channelWidth * normalVec[0]
|
||||
yA2 = yA + channelWidth * normalVec[1]
|
||||
yB2 = yB + channelWidth * normalVec[1]
|
||||
|
||||
|
||||
nA = max(NodeList) + 1
|
||||
nB = max(NodeList) + 2
|
||||
|
||||
node1Dduplicate[nAe].append(nA)
|
||||
node1Dduplicate[nBe].append(nB)
|
||||
|
||||
node2nodevalue[nA] = nAe
|
||||
node2nodevalue[nB] = nBe
|
||||
|
||||
|
||||
NodeList.append(nA)
|
||||
NodeList.append(nB)
|
||||
nodeDict[nA] = [xA2, yA2, -1.01]
|
||||
nodeDict[nB] = [xB2, yB2, -1.01]
|
||||
|
||||
newEle = max(ElementList) + 1
|
||||
ElementList .append(newEle)
|
||||
ElementDict[ElementNum] = [nAe, nA, nBe]
|
||||
ElementDict[newEle] = [nA, nB, nBe]
|
||||
|
||||
def RMA2toSerafin():
|
||||
f = open('{}.slf'.format(RMAfilename), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
str='{0: >80}'.format('SERAFIN ')
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(constName)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for c in constName:
|
||||
f.write(struct.pack(">l",32))
|
||||
str='{0: <32}'.format(c)
|
||||
f.write(str.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",R.year))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
while R.next():
|
||||
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",R.time * 3600))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.xvel[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.xvel[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.yvel[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.yvel[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.depth[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.depth[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
f.write(struct.pack(">f",R.elevation[node2nodevalue[key]]))
|
||||
else:
|
||||
f.write(struct.pack(">f",R.elevation[key]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
f.close()
|
||||
|
||||
def RMA11toSerafin(option=1):
|
||||
f = open('{}_{}.slf'.format(run,option), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
strtemp='{0: >80}'.format('SERAFIN ')
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(constName)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for idx,c in enumerate(constName):
|
||||
f.write(struct.pack(">l",32))
|
||||
strtemp='{0: <32}'.format(str(idx)+c)
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",startDate.year))
|
||||
f.write(struct.pack(">l",startDate.month))
|
||||
f.write(struct.pack(">l",startDate.day))
|
||||
f.write(struct.pack(">l",startDate.hour))
|
||||
f.write(struct.pack(">l",startDate.minute))
|
||||
f.write(struct.pack(">l",startDate.second))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
while R.next():
|
||||
#for i in range(3):
|
||||
currentDate = datetime(R.year,1,1) + timedelta(hours = R.time)
|
||||
timeCurrentStep = currentDate - startDate
|
||||
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",timeCurrentStep.total_seconds()))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
writeConst('X-VEL',key,f)
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
writeConst('Y-VEL',key,f)
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
writeConst('DEPTH',key,f)
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
writeConst('FREE SURFACE',key,f)
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
for c in constNum:
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
writeConst(c,key,f)
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
f.close()
|
||||
|
||||
|
||||
def RMA10toSerafin(option=1):
|
||||
f = open('{}_{}.slf'.format(run,option), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
strtemp='{0: >80}'.format('SERAFIN ')
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(constName)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for idx,c in enumerate(constName):
|
||||
f.write(struct.pack(">l",32))
|
||||
strtemp='{0: <32}'.format(str(idx)+c)
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",startDate.year))
|
||||
f.write(struct.pack(">l",startDate.month))
|
||||
f.write(struct.pack(">l",startDate.day))
|
||||
f.write(struct.pack(">l",startDate.hour))
|
||||
f.write(struct.pack(">l",startDate.minute))
|
||||
f.write(struct.pack(">l",startDate.second))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
while R.next():
|
||||
#for i in range(3):
|
||||
currentDate = datetime(R.year,1,1) + timedelta(hours = R.time)
|
||||
timeCurrentStep = currentDate - startDate
|
||||
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",timeCurrentStep.total_seconds()))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
|
||||
|
||||
# f.write(struct.pack(">l",len(NodeList)*4))
|
||||
# for key, value in nodeDict.items():
|
||||
# writeConst('X-VEL',key,f)
|
||||
# f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
|
||||
|
||||
for c in constName:
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
data = [None] * len(nodeDict)
|
||||
idx = 0
|
||||
for key, value in nodeDict.items():
|
||||
data[idx] = writeConst(c,key)
|
||||
idx += 1
|
||||
f.write(b''.join(data))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
f.close()
|
||||
|
||||
def writeConst(param,key):
|
||||
#get the surface node number
|
||||
# ['X-VEL','Y-VEL','DEPTH','FREE SURFACE','Constituent 1'......]
|
||||
|
||||
if param == 'X-VEL':
|
||||
tempR = R.xvel
|
||||
elif param == 'Y-VEL':
|
||||
tempR = R.yvel
|
||||
elif param == 'Z-VEL':
|
||||
tempR = R.zvel
|
||||
elif param == 'DEPTH':
|
||||
tempR = R.depth
|
||||
elif param == 'FREE SURFACE':
|
||||
tempR = R.elevation
|
||||
elif param == 'TEMPERATURE':
|
||||
tempR = R.temperature
|
||||
elif param == 'SALINITY':
|
||||
tempR = R.salinity
|
||||
elif param == 'SSED':
|
||||
tempR = R.sussed
|
||||
else:
|
||||
tempR = R.constit[c]
|
||||
|
||||
if key in node2nodevalue.keys():
|
||||
key = node2nodevalue[key]
|
||||
#based on the surface node number and option selected calculate the value
|
||||
if option == 'Max':
|
||||
tempArr = []
|
||||
for n in nodelayer[key]:
|
||||
tempArr.append(tempR[n])
|
||||
tempVal1 = max(tempArr)
|
||||
elif isinstance(option, int):
|
||||
temppos = (option - 1 ) * 2
|
||||
temppos2 = min([temppos,len(nodelayer[key])-1])
|
||||
tempVal1 = tempR[nodelayer[key][temppos2]]
|
||||
else:
|
||||
tempVal1 = tempR[key]
|
||||
#save the value in the selaphin file
|
||||
#f.write(struct.pack(">f",tempVal1))
|
||||
return struct.pack(">f",tempVal1)
|
||||
|
||||
|
||||
def writePCTL(c,pctl):
|
||||
|
||||
f = open('{}_{}_{}.slf'.format(run,c,option), 'wb')
|
||||
|
||||
f.write(struct.pack(">l",80))
|
||||
strtemp='{0: >80}'.format('SERAFIN ')
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",80))
|
||||
|
||||
f.write(struct.pack(">l",8))
|
||||
f.write(struct.pack(">l",len(pctl)))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",8))
|
||||
|
||||
for idx,p in enumerate(pctl):
|
||||
f.write(struct.pack(">l",32))
|
||||
strtemp='{0: <32}'.format(str(idx)+ 'pct ' + str(p))
|
||||
f.write(strtemp.encode('ascii'))
|
||||
f.write(struct.pack(">l",32))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",40))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",40))
|
||||
|
||||
f.write(struct.pack(">l",24))
|
||||
f.write(struct.pack(">l",startDate.year))
|
||||
f.write(struct.pack(">l",startDate.month))
|
||||
f.write(struct.pack(">l",startDate.day))
|
||||
f.write(struct.pack(">l",startDate.hour))
|
||||
f.write(struct.pack(">l",startDate.minute))
|
||||
f.write(struct.pack(">l",startDate.second))
|
||||
f.write(struct.pack(">l",24))
|
||||
|
||||
f.write(struct.pack(">l",16))
|
||||
f.write(struct.pack(">l",len(ElementList)))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
f.write(struct.pack(">l",3))
|
||||
f.write(struct.pack(">l",1))
|
||||
f.write(struct.pack(">l",16))
|
||||
|
||||
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
for el in ElementList:
|
||||
for nd in ElementDict[el]:
|
||||
|
||||
f.write(struct.pack(">l",nodeOrdered[nd]))
|
||||
f.write(struct.pack(">l",len(ElementList)*3*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
for i in range(0,len(NodeList)):
|
||||
f.write(struct.pack(">l",0))
|
||||
f.write(struct.pack(">l",len(NodeList)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[0]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
for key, value in nodeDict.items():
|
||||
f.write(struct.pack(">f",value[1]))
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
columnName = [node2nodevalue[key] if key in node2nodevalue.keys() else key for key, value in nodeDict.items()]
|
||||
|
||||
global dfAll
|
||||
dfAll = {}
|
||||
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
key = node2nodevalue[key]
|
||||
|
||||
dfAll[key] = np.array([])
|
||||
|
||||
cpt = 0
|
||||
|
||||
|
||||
currentDate = datetime(R.year,1,1) + timedelta(hours = R.time)
|
||||
timeCurrentStep = currentDate - startDate
|
||||
|
||||
f.write(struct.pack(">l",4))
|
||||
f.write(struct.pack(">f",timeCurrentStep.total_seconds()))
|
||||
f.write(struct.pack(">l",4))
|
||||
|
||||
while R.next():
|
||||
|
||||
tempR = R.constit[c]
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
key = node2nodevalue[key]
|
||||
|
||||
tempArr = []
|
||||
for n in nodelayer[key]:
|
||||
tempArr.append(tempR[n])
|
||||
tempVal1 = max(tempArr)
|
||||
dfAll[key] = np.append(dfAll[key],tempVal1)
|
||||
|
||||
cpt = cpt + 1
|
||||
|
||||
for p in pctl:
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
for key, value in nodeDict.items():
|
||||
if key in node2nodevalue.keys():
|
||||
key = node2nodevalue[key]
|
||||
f.write(struct.pack(">f",np.quantile(dfAll[key],p)))
|
||||
|
||||
f.write(struct.pack(">l",len(NodeList)*4))
|
||||
|
||||
def readNodeLayers(fname,Nnodes):
|
||||
fn = open(fname)
|
||||
lines = fn.readlines()
|
||||
|
||||
surfaceNode = int(lines[0])
|
||||
nodelayer[surfaceNode] = []
|
||||
for l in lines:
|
||||
l = int(l)
|
||||
if (l != surfaceNode) & (l < Nnodes + 1):
|
||||
surfaceNode = l
|
||||
nodelayer[surfaceNode] = []
|
||||
nodelayer[surfaceNode].append(l)
|
||||
|
||||
# In[67]:
|
||||
|
||||
#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation
|
||||
nodelayer = {}
|
||||
NodeList = []
|
||||
ElementList = []
|
||||
ElementDict = {}
|
||||
nodeDict = {}
|
||||
node1Dduplicate = {} #Original Number: List of Duplicates
|
||||
node2nodevalue = {} #link between the node number and the node value to use
|
||||
#(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1},
|
||||
#but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated
|
||||
#node will be the same as the original so we might have {2050: 10})
|
||||
|
||||
with open(meshFilename) as f:
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
line = f.readline()
|
||||
|
||||
cpt = 1
|
||||
while line and line != ' 9999\n':
|
||||
|
||||
#temp = line.split()
|
||||
|
||||
temp = re.findall('.....',line)
|
||||
ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9) if int(temp[i]) != 0 and int(temp[9]) < 100]
|
||||
ElementList.append(int(temp[0]))
|
||||
line = f.readline()
|
||||
|
||||
for key, value in ElementDict.items():
|
||||
NodeList.extend(value)
|
||||
|
||||
NodeList = list(set(NodeList))
|
||||
|
||||
line = f.readline()
|
||||
while line and line != ' 9999\n':
|
||||
|
||||
formatFix = (10,16,20,14,10,10)
|
||||
temp = []
|
||||
for form in formatFix:
|
||||
temp.append(line[:form])
|
||||
line = line[form:]
|
||||
|
||||
#line.split()
|
||||
if int(temp[0]) in NodeList:
|
||||
nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])]
|
||||
line = f.readline()
|
||||
|
||||
print('Processing triangles')
|
||||
ElementList2 = ElementList[:]
|
||||
for e in ElementList2:
|
||||
print(e)
|
||||
if isElementTriangle(ElementDict[e]):
|
||||
triangle2Triangles(e)
|
||||
elif isElementOneD(ElementDict[e]):
|
||||
oneD2triangle(e)
|
||||
else:
|
||||
square2Triangle(e)
|
||||
|
||||
for key in list(ElementDict): #Remove Special Element 902.....
|
||||
if len(ElementDict[key]) != 3:
|
||||
print(key, ElementDict[key])
|
||||
ElementDict.pop(key)
|
||||
ElementList.remove(key)
|
||||
|
||||
nodeOrdered = {}
|
||||
cpt = 1
|
||||
for key, value in nodeDict.items():
|
||||
nodeOrdered[key] = cpt
|
||||
cpt +=1
|
||||
|
||||
# # Open and Read First Step of the RMA File and Save a Serafin
|
||||
|
||||
# In[72]:
|
||||
for run in RMAfilenames:
|
||||
for option in Options:
|
||||
RMAfilename = '{}'.format(run)
|
||||
R=rma()
|
||||
R.open(RMAfilename)
|
||||
R.next()
|
||||
startDate = datetime(R.year,1,1) + timedelta(hours = R.time)
|
||||
if R.type==b'RMA11 ':
|
||||
constName = []
|
||||
readNodeLayers(NodeLayersFile,Nnodes)
|
||||
constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE']
|
||||
for c in constNum:
|
||||
constName.append(R.constit_name[c].decode("utf-8"))
|
||||
print(c)
|
||||
|
||||
if option == 'Percentile':
|
||||
for c in constNum:
|
||||
writePCTL(c,percentiles)
|
||||
else:
|
||||
RMA11toSerafin(option)
|
||||
|
||||
if R.type=='RMA10 ':
|
||||
constName = []
|
||||
readNodeLayers(NodeLayersFile,Nnodes)
|
||||
constName = ['X-VEL','Y-VEL','Z-VEL','FREE SURFACE','SALINITY','TEMPERATURE']
|
||||
|
||||
|
||||
if option == 'Percentile':
|
||||
for c in constNum:
|
||||
writePCTL(c,percentiles)
|
||||
else:
|
||||
RMA10toSerafin(option)
|
||||
|
||||
if R.type==b'RMA2 ':
|
||||
constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE']
|
||||
RMA2toSerafin()
|
||||
|
||||
|
||||
|
||||
@ -0,0 +1,212 @@
|
||||
from struct import unpack
|
||||
#from pylab import *
|
||||
import sys
|
||||
from array import array
|
||||
|
||||
|
||||
# Updated 10/11/15 BMM : Added mesh.
|
||||
# Updated 10/11/15 MD : Added dictionnary initialisation in mesh.
|
||||
# Update 8/11/2019 MD : Improve performance (used array instead of unpack)
|
||||
|
||||
class rma:
|
||||
def __init__(self):
|
||||
self.file = 'not opened yet'
|
||||
|
||||
# May want to do some other things here later
|
||||
|
||||
def open(self, filename):
|
||||
self.file = open(('%s.rma' % (filename)), 'rb')
|
||||
self.header = self.file.read(1000).decode("utf-8")
|
||||
self.type = self.header[0:10]
|
||||
self.title = self.header[100:172]
|
||||
self.geometry = self.header[200:300]
|
||||
self.num_nodes = int(self.header[40:50])
|
||||
self.num_elements = int(self.header[50:60])
|
||||
if self.type == 'RMA11 ':
|
||||
self.num_constits = int(self.header[60:70])
|
||||
if len(self.header[80:90].strip()) == 0:
|
||||
self.num_sedlayers = 0
|
||||
else:
|
||||
self.num_sedlayers = int(self.header[80:90])
|
||||
self.constit_name = []
|
||||
self.constit_name.append("NULL")
|
||||
i = 1
|
||||
print(self.num_constits)
|
||||
print(self.header[300:1000])
|
||||
while i <= self.num_constits:
|
||||
# print self.header[300:400]
|
||||
self.constit_name.append(self.header[300 + (i - 1) * 8:308 + (i - 1) * 8])
|
||||
if self.header[300 + (i - 1) * 8:308 + (i - 1) * 8] == " SSED ":
|
||||
self.constit_name.append(" BSHEAR")
|
||||
self.constit_name.append("BedThick")
|
||||
j = 1
|
||||
print(self.num_sedlayers)
|
||||
while j <= self.num_sedlayers:
|
||||
self.constit_name.append(" L%dThick" % j)
|
||||
j = j + 1
|
||||
i = i + 1
|
||||
# print i
|
||||
|
||||
# print self.num_constits
|
||||
# print i
|
||||
# if self.num_sedlayers > 0:
|
||||
# self.num_constits=self.num_constits+2+self.num_sedlayers
|
||||
|
||||
# Header format in RMA2
|
||||
|
||||
# HEADER(1:10) ='RMA2 '
|
||||
# HEADER(11:20)=DATEC = Date of run
|
||||
# HEADER(21:30)=TIMEC = Time of run
|
||||
# HEADER(31:40)=ZONEC = Time zone
|
||||
# HEADER(41:50)=NP = Number of nodes
|
||||
# HEADER(51:60)=NE = Number of elements
|
||||
# HEADER(101:172)=TITLE = Title line
|
||||
# HEADER(201:300)=FNAME = File name of binary geometry
|
||||
|
||||
# Header format in RMA11
|
||||
# HEADER(1:10) - Model
|
||||
# HEADER(11:20) - DATEC of model runtime
|
||||
# HEADER(21:30) - TIMEC of model runtime
|
||||
# HEADER(31:40) - ZONEC of model runtime
|
||||
# HEADER(41:90) - '(5I10)' NP,NE,NQAL,IGS,LGS or '(5I10)' NP,NE,NQAL,ISEDS,NLAYT
|
||||
# HEADER(101:172) - TITLE
|
||||
# HEADER(301:540) - '(30A8)' CLABL
|
||||
# HEADER(201:248) - FNAMD directory of the geometry
|
||||
# HEADER(251:298) - FNAM2 filename of the geometry
|
||||
|
||||
def mesh(self, filename):
|
||||
self.meshfile = open(('%s.rm1' % (filename)), 'r')
|
||||
l = self.meshfile.readline()
|
||||
|
||||
while l[0:5] != ' 9999' and len(l) != 5:
|
||||
l = self.meshfile.readline()
|
||||
l = self.meshfile.readline()
|
||||
self.node_e = {}
|
||||
self.node_n = {}
|
||||
self.node_z = {}
|
||||
while l[0:10] != ' 9999' and len(l) != 10:
|
||||
ll = l.split()
|
||||
self.node_e[ll[0]] = ll[1]
|
||||
self.node_n[ll[0]] = ll[2]
|
||||
self.node_z[ll[0]] = ll[3]
|
||||
l = self.meshfile.readline
|
||||
|
||||
def next(self):
|
||||
# This reads the next timestep and populates the variables.
|
||||
# Reset all of the variables
|
||||
self.time = 0.0
|
||||
self.year = 0
|
||||
self.xvel = []
|
||||
self.yvel = []
|
||||
self.zvel = []
|
||||
self.depth = []
|
||||
self.elevation = []
|
||||
self.temperature = []
|
||||
self.salinity = []
|
||||
self.sussed = []
|
||||
|
||||
# Add in an entry to fill position 0. Important since RMA files start numbering nodes from 1.
|
||||
self.xvel.append(-999)
|
||||
self.yvel.append(-999)
|
||||
self.zvel.append(-999)
|
||||
self.depth.append(-999)
|
||||
self.elevation.append(-999)
|
||||
self.temperature.append(-999)
|
||||
self.salinity.append(-999)
|
||||
self.sussed.append(-999)
|
||||
|
||||
if self.type == 'RMA2 ':
|
||||
# WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP)
|
||||
t = self.file.read(12)
|
||||
if t:
|
||||
a = unpack('fii', t)
|
||||
self.time = a[0]
|
||||
np = int(a[1])
|
||||
self.year = a[2]
|
||||
if (np != self.num_nodes):
|
||||
print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
|
||||
#b = unpack('%df' % 5 * np, self.file.read(20 * np))
|
||||
b = array('f')
|
||||
b.fromfile(self.file, 5 * np)
|
||||
|
||||
# b = np.fromfile(self.file, dtype=float,count = 5 * np)
|
||||
|
||||
tempA = [(x - 1) * 3 for x in range(1, np + 1)]
|
||||
tempB = [np * 3 + (x - 1) for x in range(1, np + 1)]
|
||||
self.xvel.extend([b[i] for i in tempA])
|
||||
self.yvel.extend([b[i + 1] for i in tempA])
|
||||
self.depth.extend([b[i + 2] for i in tempA])
|
||||
self.elevation.extend([b[i] for i in tempB])
|
||||
|
||||
if self.type == 'RMA11 ':
|
||||
self.constit = []
|
||||
c = 0
|
||||
# Start at zero to leave an empty array space. Constits numbered 1 .. num_constits
|
||||
while c <= self.num_constits:
|
||||
self.constit.append([])
|
||||
# Put a null into the 0 position so that RMA node numbers are in the same numbered array position.
|
||||
self.constit[c].append(-999)
|
||||
c = c + 1
|
||||
# READ(file1,END=100) TETT1,NQAL,NP,IYRR, ((VEL(K,J),J=1,NP),K=1,3), (wd(j),j=1,np), (wsel(j),j=1,np), ((TCON1(K,J),J=1,NP),K=1,NQAL-5)
|
||||
t = self.file.read(16)
|
||||
if t:
|
||||
a = unpack('fiii', t)
|
||||
self.time = a[0]
|
||||
nqal = int(a[1])
|
||||
np = int(a[2])
|
||||
self.year = a[3]
|
||||
if ((nqal - 5) != (self.num_constits)):
|
||||
print("Warning - NQAL-5 (%d) on this timestep does not match header (%d)" % (
|
||||
nqal - 5, self.num_constits))
|
||||
if (np != self.num_nodes):
|
||||
print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
|
||||
|
||||
b = array('f')
|
||||
b.fromfile(self.file, nqal * np)
|
||||
#b = unpack('%df' % nqal * np, self.file.read(4 * nqal * np)) #4 because single precision float type have 4 bytes
|
||||
|
||||
tempA = [(x - 1) * 3 for x in range(1, np + 1)]
|
||||
tempB = [np * 3 + (x - 1) for x in range(1, np + 1)]
|
||||
tempC = [np * 4 + (x - 1) for x in range(1, np + 1)]
|
||||
self.xvel.extend([b[i] for i in tempA])
|
||||
self.yvel.extend([b[i + 1] for i in tempA])
|
||||
self.zvel.extend([b[i + 2] for i in tempA])
|
||||
self.depth.extend([b[i] for i in tempB])
|
||||
self.elevation.extend([b[i] for i in tempC])
|
||||
for c in range(1, self.num_constits + 1):
|
||||
(self.constit[c].extend([b[np * ((c - 1) + 5) + (x - 1)] for x in range(1, np + 1)]))
|
||||
if self.type == 'RMA10 ':
|
||||
# WRITE(IRMAFM) TETT,NP,NDF,NE,IYRR,((VSING(K,J),K=1,NDF),VVEL(J),WSLL(J),J=1,NP),(DFCT(J),J=1,NE),(VSING(7,J),J=1,NP)
|
||||
# WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP)
|
||||
t = self.file.read(20)
|
||||
if t:
|
||||
a = unpack('fiiii', t)
|
||||
self.time = a[0]
|
||||
np = a[1]
|
||||
ndf = 6
|
||||
ne = a[3]
|
||||
self.year = a[4]
|
||||
if (np != self.num_nodes):
|
||||
print("Warning - NP1 (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
|
||||
tempRead = np * (3 + ndf) + ne
|
||||
b = array('f')
|
||||
b.fromfile(self.file, tempRead)
|
||||
# b = numpy.fromfile(self.file.read(4 * tempRead), dtype=float)
|
||||
#b = unpack('%df' % tempRead, self.file.read(4 * tempRead))
|
||||
i = 1
|
||||
while i <= np:
|
||||
self.xvel.append(b[0 + (i - 1) * 8])
|
||||
self.yvel.append(b[1 + (i - 1) * 8])
|
||||
self.depth.append(b[2 + (i - 1) * 8])
|
||||
self.salinity.append(b[3 + (i - 1) * 8])
|
||||
self.temperature.append(b[4 + (i - 1) * 8])
|
||||
self.sussed.append(b[5 + (i - 1) * 8])
|
||||
self.zvel.append(b[6 + (i - 1) * 8])
|
||||
self.elevation.append(b[7 + (i - 1) * 8])
|
||||
i = i + 1
|
||||
|
||||
if len(self.xvel) == 1:
|
||||
# Note that this is a 1 now as we've filled the zero array position
|
||||
return False
|
||||
else:
|
||||
return True
|
||||
Loading…
Reference in New Issue