From a2f193a834b49c81d4d472e400a0f815d3da2e97 Mon Sep 17 00:00:00 2001 From: Mathieu Deiber Date: Mon, 21 Oct 2019 15:01:10 +1100 Subject: [PATCH] Create script for RMA11 using RMA10 result file --- RMA2SERAPHIN_3DRMA11.py | 337 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 337 insertions(+) create mode 100644 RMA2SERAPHIN_3DRMA11.py diff --git a/RMA2SERAPHIN_3DRMA11.py b/RMA2SERAPHIN_3DRMA11.py new file mode 100644 index 0000000..92afcc4 --- /dev/null +++ b/RMA2SERAPHIN_3DRMA11.py @@ -0,0 +1,337 @@ + +# coding: utf-8 + +# In[64]: + +import struct +import matplotlib.pyplot as plt +import math +from py_rmatools import rma +import re +from datetime import datetime, timedelta +plt.rcParams.update({'figure.max_open_warning': 0}) + + +# In[65]: + +meshFilename = 'bub005c.rm1' +channelWidth = 100 +RMAfilename = 'BUB023_WQ' + +#If RMA11 +constNum = [1] + +#If RMA11 3D + +Options = ['Max',1,4,8] #Option available 'Max' , an integer correspond to the layer number to extract the results, if 2D use options = [1] +NodeLayersFile = 'NodeLayers.txt' #the node layer file was created by getting the nodes number from the RMA Outfile (just run RMA11 for a few steps) +Nnodes = 5769 #Number of nodes (including mid side node) - thisn umber also correspond to the last node number in the nodelayers file +# 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(option=1): + f = open('{}_{}.slf'.format(RMAfilename,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 writeConst(param,key,f): +#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 == 'DEPTH': + tempR = R.depth + elif param == 'FREE SURFACE': + tempR = R.elevation + 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[key]) + tempVal1 = max(tempArr) + elif isinstance(option, int): + temppos = (option - 1 ) * 2 + temppos2 = min([temppos,len(nodelayer[key])]) + + tempVal1 = tempR[nodelayer[key][temppos2]] + else: + tempVal1 = tempR[key] + #save the value in the selaphin file + f.write(struct.pack(">f",tempVal1)) + + +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,2) 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() + + +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 + +Nnodes = max(NodeList) +# # Open and Read First Step of the RMA File and Save a Serafin + +# In[72]: +for option in Options: + 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) + + RMA11toSerafin(option) + + +