mathieud
/
RMA2SERAPHIN
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

First commit

Browse Source
master
mathieud 7 years ago
parent dd481ca4c0
commit a15701ec53
5 changed files with 1205 additions and 0 deletions
Show Stats Download Patch File Download Diff File

40
RMA11_ConstituentRev2.py
Unescape Escape View File

@ -0,0 +1,40 @@
from py_rmatools import rma
startyear=2013
runNumber = 'HWQ041'
DictNodeConst = {
1: [94,204,2281,213,2466,209,212,166,2468,2257,2256,194,2259,2471,2261,2252,2253,2467,2260,2262,197,2251,2255,2258,195,93,205,201,2254,200,2464,202,206,51,199,207,2465,2469,193,211,2263,196,208,2470,210,198,203],
2: [606,89,168,2445,600,614,166,2283,2804,2460,172,2824,2830,2808,2450,176,2452,2451,82,2801,2454,186,616,2447,2820,611,2455,2821,593,167,83,87,2459,605,84,621,618,2802,625,596,2813,2289,2803,180,2290,185,2809,2443,622,187,2811,184,626,617,601,80,85,81,619,88,181,592,2456,2287,603,610,2823,2825,177,2818,2444,179,602,170,2462,2805,2814,2292,2812,2442,2286,2798,595,2463,2822,2819,624,183,2807,2284,613,609,2291,2458,608,175,2806,173,599,598,2827,2446,2448,174,90,2461,2810,2829,604,2285,2800,171,2817,2826,182,2288,2828,91,2799,2453,615,627,2831,169,2797,178,612,620,623,607,597,2816,2796,594,2449,2815,2457,2282,86],
3: [498,507,516,2974,518,499,2942,531,512,33,64,485,34,41,2957,38,46,533,2947,524,2936,51,2915,2963,522,514,39,2928,2945,2932,2955,2960,2980,2927,489,2943,506,496,50,2967,508,2916,2956,2925,2930,2889,495,2970,520,2979,2977,2964,40,2968,511,2952,591,523,500,528,32,2950,525,502,488,2924,2934,2962,43,45,35,505,509,63,530,36,492,521,503,2971,2921,2917,519,2944,2946,2926,48,2922,2961,2938,2923,497,487,493,47,2975,2937,513,2972,504,517,501,2931,2953,2929,486,2920,2958,494,2966,515,2918,2983,2933,510,42,2935,2982,2940,2981,2939,37,490,527,2954,49,2965,2978,44,2973,2951,526,2941,491,2949,2969,2976,2919,2948,532,2959,484,31,529,2912],
4: [76,2871,2700,2768,2787,2715,2778,580,2271,575,1397,2850,436,2776,379,2844,403,2712,2713,349,2767,555,2703,543,2719,2842,2710,424,369,67,589,2745,2688,346,2848,377,2743,2763,400,2758,2277,2788,2874,428,1386,2741,417,332,374,361,2689,341,426,551,548,384,432,2782,401,382,577,542,574,2834,563,2781,2747,1396,66,2885,2739,2840,427,2769,591,339,351,2838,2714,547,4246,2861,2270,4249,2273,2882,583,2731,2869,347,381,2754,2722,430,2694,390,1399,2854,4252,2865,343,380,581,73,364,435,2702,1400,4250,587,69,567,2863,2784,75,406,558,534,2760,2888,2766,421,2852,545,2868,414,536,71,2867,1390,2441,557,410,416,2272,429,371,2740,541,2693,70,438,2735,423,2771,366,569,2870,2711,2790,2780,1385,4242,362,2732,2761,572,2880,2785,2695,2279,2750,590,1401,383,576,2691,394,544,4247,407,2851,1394,419,2707,2886,333,378,404,2728,392,415,368,395,2841,358,340,571,2846,367,385,2720,165,2847,586,359,344,579,2756,4257,2864,564,2786,422,437,2697,556,584,2699,2872,434,570,2859,72,549,537,2278,386,338,2274,2836,2698,2752,2753,405,582,4258,2862,397,2751,356,409,2705,337,2727,1392,398,540,2887,2724,2793,2757,2835,553,552,2876,2884,4254,538,2837,331,2764,373,352,1402,2855,1393,2773,4245,2269,2729,1398,560,1388,389,425,2762,2775,376,2856,2770,2730,573,74,354,335,2853,336,2792,2744,350,345,2748,4244,357,562,375,4241,418,578,1389,2765,2873,2736,433,2832,2881,2783,393,2875,2794,2791,561,2726,2275,363,1391,2268,2721,2738,2774,2718,2833,370,550,402,2759,4243,559,4255,565,2734,2696,2742,68,2857,2789,413,365,2843,2879,431,2723,408,342,4248,2772,2883,2717,2795,1395,4256,588,1387,360,420,2746,2877,2866,372,566,2733,2755,2845,2737,2725,412,348,2878,539,2839,554,2860,77,535,2779,2276,399,387,4251,2692,2749,411,2708,2701,334,2849,2690,396,4253,568,2704,2716,546,2777,2858,2706,585,353,2709,391,355,388],
5: [2218,2244,24,451,2232,2212,2893,442,26,483,2897,444,2231,1,2248,466,2892,452,2,3,464,453,2215,2233,454,2210,2234,2203,2896,2216,478,2903,440,470,2245,30,4,21,459,4364,19,2236,2891,445,2206,2899,460,1373,17,463,25,2207,2213,480,2221,443,2222,476,2909,2209,2908,471,6,20,13,15,482,2900,458,2895,456,18,2238,2247,2904,2239,450,457,2902,2201,2223,16,1369,27,2217,2205,2243,5994,2220,2890,9,1374,2911,12,2913,469,2240,2905,441,2246,2901,2894,2225,449,2914,455,2214,2280,472,11,2211,462,2204,2202,2241,1372,2237,5,481,448,475,477,2907,22,2910,7,2228,465,446,2242,14,2230,5993,8,461,2906,2226,468,2250,28,447,2898,439,2249,10,23,2208,467,2219,474,29,2227,2224,479,473,484,2235]
}
for key in DictNodeConst:
time=[]
node=DictNodeConst[key]
I_const=key
filename1='{}_{}_WQ'.format(runNumber,startyear)
filename2='{}_{}.csv'.format(runNumber,I_const)
target = open(filename2, 'w')
target.write("Year, Hour")
for n in node:
target.write(", {}".format(n))
target.write("\n")
R=rma()
R.open(filename1)
while R.next():
target.write("{}, {}".format(startyear,R.time))
for ii in range(0,len(node)):
target.write(", {}".format(R.constit[I_const][node[ii]]))
target.write("\n")
target.close()

560
RMA2SERAPHIN.ipynb
Unescape Escape View File

@ -0,0 +1,560 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import struct\n",
"import matplotlib.pyplot as plt\n",
"import math\n",
"from py_rmatools import rma\n",
"\n",
"%matplotlib qt\n",
"plt.rcParams.update({'figure.max_open_warning': 0})"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"meshFilename = 'ck_034.rm1'\n",
"channelWidth = 100\n",
"RMAfilename = 'CK011_WQ'\n",
"\n",
"#If RMA11\n",
"constNum = [1]"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def isElementOneD(nodelist):\n",
" if len(nodelist) == 2:\n",
" return True\n",
" return False\n",
"\n",
"def isElementSquare(nodelist):\n",
" if len(nodelist) == 4:\n",
" return True\n",
" return False\n",
"\n",
"def square2Triangle(ElementNum):\n",
" nodelist = ElementDict[ElementNum]\n",
" if isElementSquare(nodelist):\n",
" ElementDict[ElementNum] = [nodelist[0], nodelist[1], nodelist[2]]\n",
" ElementList.append(max(ElementList) + 1)\n",
" ElementDict[ElementList[-1]]= [nodelist[0], nodelist[2], nodelist[3]]\n",
" \n",
"def oneD2triangle(ElementNum):\n",
" if isElementOneD(ElementDict[ElementNum]):\n",
" nAe = ElementDict[ElementNum][0] #nAe Node A existing\n",
" nBe = ElementDict[ElementNum][1]\n",
" \n",
" if not nAe in node1Dduplicate: node1Dduplicate[nAe] = []\n",
" if not nBe in node1Dduplicate: node1Dduplicate[nBe] = []\n",
" \n",
" xA = nodeDict[nAe][0]\n",
" xB = nodeDict[nBe][0]\n",
" yA = nodeDict[nAe][1]\n",
" yB = nodeDict[nBe][1]\n",
" \n",
" normalVec = [-(yB - yA),(xB - xA)]\n",
" dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2)\n",
" normalVec[0] = normalVec[0] / dist\n",
" normalVec[1] = normalVec[1] / dist\n",
" xA2 = xA + channelWidth * normalVec[0]\n",
" xB2 = xB + channelWidth * normalVec[0]\n",
" yA2 = yA + channelWidth * normalVec[1]\n",
" yB2 = yB + channelWidth * normalVec[1]\n",
" \n",
"\n",
" nA = max(NodeList) + 1\n",
" nB = max(NodeList) + 2\n",
" \n",
" node1Dduplicate[nAe].append(nA)\n",
" node1Dduplicate[nBe].append(nB) \n",
" \n",
" node2nodevalue[nA] = nAe\n",
" node2nodevalue[nB] = nBe\n",
" \n",
" \n",
" NodeList.append(nA)\n",
" NodeList.append(nB)\n",
" nodeDict[nA] = [xA2, yA2, -1.01]\n",
" nodeDict[nB] = [xB2, yB2, -1.01]\n",
" \n",
" newEle = max(ElementList) + 1\n",
" ElementList .append(newEle)\n",
" ElementDict[ElementNum] = [nAe, nA, nBe]\n",
" ElementDict[newEle] = [nA, nB, nBe]\n",
" \n",
"def RMA11toSerafin():\n",
" f = open('{}.slf'.format(RMAfilename), 'wb')\n",
"\n",
" f.write(struct.pack(\">l\",80))\n",
" str='{0: >80}'.format('SERAFIN ')\n",
" f.write(str.encode('ascii'))\n",
" f.write(struct.pack(\">l\",80))\n",
"\n",
" f.write(struct.pack(\">l\",8))\n",
" f.write(struct.pack(\">l\",len(constNum)))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",8))\n",
"\n",
" for c in constName:\n",
" f.write(struct.pack(\">l\",32))\n",
" str='{0: <32}'.format(c)\n",
" f.write(str.encode('ascii'))\n",
" f.write(struct.pack(\">l\",32))\n",
"\n",
"\n",
" f.write(struct.pack(\">l\",40))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",40))\n",
"\n",
" f.write(struct.pack(\">l\",24))\n",
" f.write(struct.pack(\">l\",R.year))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",24))\n",
"\n",
" f.write(struct.pack(\">l\",16))\n",
" f.write(struct.pack(\">l\",len(ElementList)))\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
" f.write(struct.pack(\">l\",3))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",16))\n",
"\n",
"\n",
" f.write(struct.pack(\">l\",len(ElementList)*3*4))\n",
" for el in ElementList:\n",
" for nd in ElementDict[el]:\n",
" f.write(struct.pack(\">l\",nodeOrdered[nd])) \n",
" f.write(struct.pack(\">l\",len(ElementList)*3*4))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
" for i in range(0,len(NodeList)):\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" f.write(struct.pack(\">f\",value[0]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" f.write(struct.pack(\">f\",value[1]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
"\n",
"\n",
" while R.next():\n",
" #for i in range(3):\n",
" f.write(struct.pack(\">l\",4))\n",
" f.write(struct.pack(\">f\",R.time * 3600))\n",
" f.write(struct.pack(\">l\",4))\n",
"\n",
" for c in constNum:\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" if key in node2nodevalue.keys():\n",
" f.write(struct.pack(\">f\",R.constit[c][node2nodevalue[key]]))\n",
" else:\n",
" f.write(struct.pack(\">f\",R.constit[c][key]))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
" R.next()\n",
"\n",
" f.close() \n",
"\n",
" \n",
" \n",
"def RMA2toSerafin():\n",
" f = open('{}.slf'.format(RMAfilename), 'wb')\n",
"\n",
" f.write(struct.pack(\">l\",80))\n",
" str='{0: >80}'.format('SERAFIN ')\n",
" f.write(str.encode('ascii'))\n",
" f.write(struct.pack(\">l\",80))\n",
"\n",
" f.write(struct.pack(\">l\",8))\n",
" f.write(struct.pack(\">l\",len(constName)))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",8))\n",
"\n",
" for c in constName:\n",
" f.write(struct.pack(\">l\",32))\n",
" str='{0: <32}'.format(c)\n",
" f.write(str.encode('ascii'))\n",
" f.write(struct.pack(\">l\",32))\n",
"\n",
"\n",
" f.write(struct.pack(\">l\",40))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",40))\n",
"\n",
" f.write(struct.pack(\">l\",24))\n",
" f.write(struct.pack(\">l\",R.year))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",24))\n",
"\n",
" f.write(struct.pack(\">l\",16))\n",
" f.write(struct.pack(\">l\",len(ElementList)))\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
" f.write(struct.pack(\">l\",3))\n",
" f.write(struct.pack(\">l\",1))\n",
" f.write(struct.pack(\">l\",16))\n",
"\n",
"\n",
" f.write(struct.pack(\">l\",len(ElementList)*3*4))\n",
" for el in ElementList:\n",
" for nd in ElementDict[el]:\n",
" f.write(struct.pack(\">l\",nodeOrdered[nd]))\n",
" f.write(struct.pack(\">l\",len(ElementList)*3*4))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
" for i in range(0,len(NodeList)):\n",
" f.write(struct.pack(\">l\",0))\n",
" f.write(struct.pack(\">l\",len(NodeList)))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" f.write(struct.pack(\">f\",value[0]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" f.write(struct.pack(\">f\",value[1]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
"\n",
"\n",
" while R.next():\n",
"\n",
" f.write(struct.pack(\">l\",4))\n",
" f.write(struct.pack(\">f\",R.time * 3600))\n",
" f.write(struct.pack(\">l\",4))\n",
"\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" if key in node2nodevalue.keys():\n",
" f.write(struct.pack(\">f\",R.xvel[node2nodevalue[key]]))\n",
" else:\n",
" f.write(struct.pack(\">f\",R.xvel[key]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
"\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" if key in node2nodevalue.keys():\n",
" f.write(struct.pack(\">f\",R.yvel[node2nodevalue[key]]))\n",
" else:\n",
" f.write(struct.pack(\">f\",R.yvel[key]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" \n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" if key in node2nodevalue.keys():\n",
" f.write(struct.pack(\">f\",R.depth[node2nodevalue[key]]))\n",
" else:\n",
" f.write(struct.pack(\">f\",R.depth[key]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" \n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" for key, value in nodeDict.items():\n",
" if key in node2nodevalue.keys():\n",
" f.write(struct.pack(\">f\",R.elevation[node2nodevalue[key]]))\n",
" else:\n",
" f.write(struct.pack(\">f\",R.elevation[key]))\n",
" f.write(struct.pack(\">l\",len(NodeList)*4))\n",
" \n",
" \n",
" R.next()\n",
"\n",
" f.close() "
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation\n",
"NodeList = [] \n",
"ElementList = []\n",
"ElementDict = {}\n",
"nodeDict = {}\n",
"node1Dduplicate = {} #Original Number: List of Duplicates\n",
"node2nodevalue = {} #link between the node number and the node value to use \n",
" #(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1}, \n",
" #but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated \n",
" #node will be the same as the original so we might have {2050: 10})\n",
" \n",
"with open(meshFilename) as f:\n",
" line = f.readline()\n",
" line = f.readline()\n",
" line = f.readline()\n",
" line = f.readline()\n",
" \n",
" cpt = 1\n",
" while line and line != ' 9999\\n':\n",
" temp = line.split()\n",
" ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9,2) if temp[i] != '0' and int(temp[9]) < 100]\n",
" ElementList.append(int(temp[0]))\n",
" line = f.readline()\n",
" \n",
" for key, value in ElementDict.items():\n",
" NodeList.extend(value)\n",
" \n",
" NodeList = list(set(NodeList))\n",
"\n",
" line = f.readline()\n",
" while line and line != ' 9999\\n':\n",
" temp = line.split()\n",
" if int(temp[0]) in NodeList:\n",
" nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])]\n",
" line = f.readline()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for e in ElementList:\n",
" oneD2triangle(e)\n",
" square2Triangle(e)\n",
" \n",
"for key in list(ElementDict): #Remove Special Element 902.....\n",
" if len(ElementDict[key]) != 3:\n",
" print(key, ElementDict[key])\n",
" ElementDict.pop(key)\n",
" ElementList.remove(key)\n",
" \n",
"nodeOrdered = {}\n",
"cpt = 1\n",
"for key, value in nodeDict.items():\n",
" nodeOrdered[key] = cpt\n",
" cpt +=1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Open and Read First Step of the RMA File and Save a Serafin"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"R=rma()\n",
"R.open(RMAfilename)\n",
"R.next()\n",
"\n",
"if R.type==b'RMA11 ':\n",
" constName = []\n",
" for c in constNum:\n",
" constName.append(R.constit_name[c].decode(\"utf-8\"))\n",
" RMA11toSerafin()\n",
" \n",
"if R.type==b'RMA2 ':\n",
" constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] \n",
" RMA2toSerafin()\n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

442
RMA2SERAPHIN.py
Unescape Escape View File

@ -0,0 +1,442 @@
# coding: utf-8
# In[64]:
import struct
import matplotlib.pyplot as plt
import math
from py_rmatools import rma
get_ipython().magic('matplotlib qt')
plt.rcParams.update({'figure.max_open_warning': 0})
# In[65]:
meshFilename = 'ck_034.rm1'
channelWidth = 100
RMAfilename = 'CK011_WQ'
#If RMA11
constNum = [1]
# 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()
# In[ ]:
# In[ ]:
# In[75]:
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
# In[ ]:
# # 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()
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:
# In[ ]:

BIN
Seraphin files.pdf
View File

Binary file not shown.

163
py_rmatools.py
Unescape Escape View File

@ -0,0 +1,163 @@
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)
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==b'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)
while i <= self.num_constits:
#print(i, self.num_constits, i <= self.num_constits)
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
while j<=self.num_sedlayers:
self.constit_name.append(" L%dThick" % j)
j=j+1
i=i+1
#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=[]
# 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)
if self.type==b'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))
i=1
while i<=np:
self.xvel.append(b[0+(i-1)*3])
self.yvel.append(b[1+(i-1)*3])
self.depth.append(b[2+(i-1)*3])
self.elevation.append(b[np*3+(i-1)])
i=i+1
if self.type==b'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))
i=1
while i<=np:
self.xvel.append(b[0+(i-1)*3])
self.yvel.append(b[1+(i-1)*3])
self.zvel.append(b[2+(i-1)*3])
self.depth.append(b[np*3+(i-1)])
self.elevation.append(b[np*4+(i-1)])
c=1
while c<=self.num_constits:
self.constit[c].append(b[np*((c-1)+5)+(i-1)])
c=c+1
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
Write Preview
Loading…
Cancel
Save