diff --git a/RMA11_ConstituentRev2.py b/RMA11_ConstituentRev2.py new file mode 100644 index 0000000..b7f1b08 --- /dev/null +++ b/RMA11_ConstituentRev2.py @@ -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() \ No newline at end of file diff --git a/RMA2SERAPHIN.ipynb b/RMA2SERAPHIN.ipynb new file mode 100644 index 0000000..522e356 --- /dev/null +++ b/RMA2SERAPHIN.ipynb @@ -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 +} diff --git a/RMA2SERAPHIN.py b/RMA2SERAPHIN.py new file mode 100644 index 0000000..703087a --- /dev/null +++ b/RMA2SERAPHIN.py @@ -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[ ]: + + + diff --git a/Seraphin files.pdf b/Seraphin files.pdf new file mode 100644 index 0000000..3845c17 Binary files /dev/null and b/Seraphin files.pdf differ diff --git a/py_rmatools.py b/py_rmatools.py new file mode 100644 index 0000000..0a350fd --- /dev/null +++ b/py_rmatools.py @@ -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 +