from struct import * from pylab import * import sys # Updated 10/11/15 BMM : Added mesh. # Updated 10/11/15 MD : Added dictionnary initialisation in mesh. 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)) 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 = unpack('%df' % nqal * np, self.file.read(4 * nqal * np)) 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 = 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