--- /dev/null
+1731 1827 0 0 0
+195 -0.964976 0.485465 0.000000
+133 -0.970149 0.500000 0.000000
+67 -1.000000 0.500000 0.000000
+180 -0.978723 0.478723 0.000000
+196 -0.945882 0.477585 0.000000
+132 -0.940299 0.500000 0.000000
+179 -0.957447 0.457447 0.000000
+197 -0.919731 0.466897 0.000000
+131 -0.910448 0.500000 0.000000
+178 -0.936170 0.436170 0.000000
+181 0.150000 -0.500000 0.000000
+187 0.150000 -0.471429 0.000000
+198 0.119014 -0.470800 0.000000
+185 0.120000 -0.500000 0.000000
+188 0.150000 -0.442857 0.000000
+199 0.117550 -0.441479 0.000000
+189 0.150000 -0.414286 0.000000
+200 0.117242 -0.412206 0.000000
+190 0.150000 -0.385714 0.000000
+201 0.118129 -0.383064 0.000000
+191 0.150000 -0.357143 0.000000
+202 0.120210 -0.354337 0.000000
+192 0.150000 -0.328571 0.000000
+203 0.121793 -0.325840 0.000000
+186 0.150000 -0.300000 0.000000
+204 0.122717 -0.296869 0.000000
+206 0.153091 -0.267708 0.000000
+205 0.123759 -0.266752 0.000000
+193 0.183333 -0.300000 0.000000
+207 0.184560 -0.267950 0.000000
+194 0.216667 -0.300000 0.000000
+208 0.216606 -0.268111 0.000000
+1 0.250000 -0.300000 0.000000
+209 0.248350 -0.268410 0.000000
+211 0.279437 -0.298917 0.000000
+210 0.278782 -0.268564 0.000000
+3 0.250000 -0.328571 0.000000
+212 0.279638 -0.328111 0.000000
+4 0.250000 -0.357143 0.000000
+213 0.279710 -0.356801 0.000000
+5 0.250000 -0.385714 0.000000
+214 0.279677 -0.385395 0.000000
+6 0.250000 -0.414286 0.000000
+215 0.279641 -0.413938 0.000000
+7 0.250000 -0.442857 0.000000
+216 0.279666 -0.442534 0.000000
+8 0.250000 -0.471429 0.000000
+217 0.279811 -0.471266 0.000000
+2 0.250000 -0.500000 0.000000
+10 0.280000 -0.500000 0.000000
+34 1.000000 0.500000 0.000000
+68 0.970149 0.500000 0.000000
+218 0.970096 0.469780 0.000000
+66 1.000000 0.469697 0.000000
+69 0.940299 0.500000 0.000000
+219 0.940231 0.469840 0.000000
+70 0.910448 0.500000 0.000000
+220 0.910371 0.469856 0.000000
+71 0.880597 0.500000 0.000000
+221 0.880527 0.469909 0.000000
+72 0.850746 0.500000 0.000000
+222 0.850684 0.469926 0.000000
+73 0.820896 0.500000 0.000000
+223 0.820842 0.469978 0.000000
+74 0.791045 0.500000 0.000000
+224 0.791001 0.470016 0.000000
+75 0.761194 0.500000 0.000000
+225 0.761159 0.470055 0.000000
+76 0.731343 0.500000 0.000000
+226 0.731316 0.470091 0.000000
+77 0.701493 0.500000 0.000000
+227 0.701472 0.470118 0.000000
+78 0.671642 0.500000 0.000000
+228 0.671627 0.470134 0.000000
+79 0.641791 0.500000 0.000000
+229 0.641780 0.470142 0.000000
+80 0.611940 0.500000 0.000000
+230 0.611934 0.470065 0.000000
+81 0.582090 0.500000 0.000000
+231 0.582089 0.470002 0.000000
+82 0.552239 0.500000 0.000000
+232 0.552248 0.469853 0.000000
+83 0.522388 0.500000 0.000000
+233 0.522413 0.469718 0.000000
+84 0.492537 0.500000 0.000000
+234 0.492577 0.469656 0.000000
+85 0.462687 0.500000 0.000000
+235 0.462743 0.469571 0.000000
+86 0.432836 0.500000 0.000000
+236 0.432895 0.469630 0.000000
+87 0.402985 0.500000 0.000000
+237 0.403046 0.469669 0.000000
+88 0.373134 0.500000 0.000000
+238 0.373191 0.469745 0.000000
+89 0.343284 0.500000 0.000000
+239 0.343335 0.469723 0.000000
+90 0.313433 0.500000 0.000000
+240 0.313478 0.469741 0.000000
+91 0.283582 0.500000 0.000000
+241 0.283622 0.469783 0.000000
+92 0.253731 0.500000 0.000000
+242 0.253768 0.469742 0.000000
+93 0.223881 0.500000 0.000000
+243 0.223918 0.469737 0.000000
+94 0.194030 0.500000 0.000000
+244 0.194064 0.469810 0.000000
+95 0.164179 0.500000 0.000000
+245 0.164210 0.469881 0.000000
+96 0.134328 0.500000 0.000000
+246 0.134355 0.469940 0.000000
+97 0.104478 0.500000 0.000000
+247 0.104500 0.469997 0.000000
+98 0.074627 0.500000 0.000000
+248 0.074645 0.470051 0.000000
+99 0.044776 0.500000 0.000000
+249 0.044790 0.470093 0.000000
+100 0.014925 0.500000 0.000000
+250 0.014935 0.470120 0.000000
+101 -0.014925 0.500000 0.000000
+251 -0.014918 0.470135 0.000000
+102 -0.044776 0.500000 0.000000
+252 -0.044771 0.470143 0.000000
+103 -0.074627 0.500000 0.000000
+253 -0.074623 0.470146 0.000000
+104 -0.104478 0.500000 0.000000
+254 -0.104496 0.469979 0.000000
+105 -0.134328 0.500000 0.000000
+255 -0.134359 0.469935 0.000000
+106 -0.164179 0.500000 0.000000
+256 -0.164225 0.469839 0.000000
+107 -0.194030 0.500000 0.000000
+257 -0.194083 0.469797 0.000000
+108 -0.223881 0.500000 0.000000
+258 -0.223940 0.469753 0.000000
+109 -0.253731 0.500000 0.000000
+259 -0.253798 0.469819 0.000000
+110 -0.283582 0.500000 0.000000
+260 -0.283647 0.470019 0.000000
+111 -0.313433 0.500000 0.000000
+261 -0.313491 0.470151 0.000000
+112 -0.343284 0.500000 0.000000
+262 -0.343348 0.470381 0.000000
+113 -0.373134 0.500000 0.000000
+263 -0.373195 0.470395 0.000000
+114 -0.402985 0.500000 0.000000
+264 -0.403049 0.470404 0.000000
+115 -0.432836 0.500000 0.000000
+265 -0.432937 0.470297 0.000000
+116 -0.462687 0.500000 0.000000
+266 -0.462848 0.469995 0.000000
+117 -0.492537 0.500000 0.000000
+267 -0.492807 0.469765 0.000000
+118 -0.522388 0.500000 0.000000
+268 -0.522760 0.469580 0.000000
+119 -0.552239 0.500000 0.000000
+269 -0.552730 0.469652 0.000000
+120 -0.582090 0.500000 0.000000
+270 -0.582662 0.469731 0.000000
+121 -0.611940 0.500000 0.000000
+271 -0.612594 0.469999 0.000000
+122 -0.641791 0.500000 0.000000
+272 -0.642647 0.470419 0.000000
+123 -0.671642 0.500000 0.000000
+273 -0.672638 0.470543 0.000000
+124 -0.701493 0.500000 0.000000
+274 -0.702968 0.469800 0.000000
+125 -0.731343 0.500000 0.000000
+275 -0.733232 0.468149 0.000000
+126 -0.761194 0.500000 0.000000
+276 -0.765399 0.466448 0.000000
+127 -0.791045 0.500000 0.000000
+277 -0.795702 0.469563 0.000000
+128 -0.820896 0.500000 0.000000
+278 -0.822797 0.472070 0.000000
+129 -0.850746 0.500000 0.000000
+279 -0.850723 0.467038 0.000000
+130 -0.880597 0.500000 0.000000
+303 -0.887028 0.454038 0.000000
+11 0.310000 -0.500000 0.000000
+280 0.309694 -0.471123 0.000000
+12 0.340000 -0.500000 0.000000
+281 0.339686 -0.471033 0.000000
+13 0.370000 -0.500000 0.000000
+282 0.369721 -0.471055 0.000000
+14 0.400000 -0.500000 0.000000
+283 0.399756 -0.470957 0.000000
+15 0.430000 -0.500000 0.000000
+284 0.429805 -0.470822 0.000000
+16 0.460000 -0.500000 0.000000
+285 0.459835 -0.470712 0.000000
+17 0.490000 -0.500000 0.000000
+286 0.489873 -0.470554 0.000000
+18 0.520000 -0.500000 0.000000
+287 0.519906 -0.470379 0.000000
+19 0.550000 -0.500000 0.000000
+288 0.549930 -0.470248 0.000000
+20 0.580000 -0.500000 0.000000
+289 0.579948 -0.470155 0.000000
+21 0.610000 -0.500000 0.000000
+290 0.609962 -0.470091 0.000000
+22 0.640000 -0.500000 0.000000
+291 0.639972 -0.470050 0.000000
+23 0.670000 -0.500000 0.000000
+292 0.669980 -0.470026 0.000000
+24 0.700000 -0.500000 0.000000
+293 0.699986 -0.470013 0.000000
+25 0.730000 -0.500000 0.000000
+294 0.729990 -0.470006 0.000000
+26 0.760000 -0.500000 0.000000
+295 0.759994 -0.470003 0.000000
+27 0.790000 -0.500000 0.000000
+296 0.789980 -0.469916 0.000000
+28 0.820000 -0.500000 0.000000
+297 0.819956 -0.469877 0.000000
+29 0.850000 -0.500000 0.000000
+298 0.849929 -0.469833 0.000000
+30 0.880000 -0.500000 0.000000
+299 0.879924 -0.469805 0.000000
+31 0.910000 -0.500000 0.000000
+300 0.909913 -0.469792 0.000000
+32 0.940000 -0.500000 0.000000
+301 0.939922 -0.469775 0.000000
+33 0.970000 -0.500000 0.000000
+302 0.969956 -0.469740 0.000000
+9 1.000000 -0.500000 0.000000
+35 1.000000 -0.469697 0.000000
+177 -0.914894 0.414894 0.000000
+176 -0.893617 0.393617 0.000000
+304 -0.867787 0.420844 0.000000
+175 -0.872340 0.372340 0.000000
+305 -0.847479 0.396796 0.000000
+174 -0.851064 0.351064 0.000000
+306 -0.826567 0.374510 0.000000
+173 -0.829787 0.329787 0.000000
+307 -0.806205 0.352141 0.000000
+172 -0.808511 0.308511 0.000000
+308 -0.786305 0.329477 0.000000
+171 -0.787234 0.287234 0.000000
+309 -0.766015 0.307496 0.000000
+170 -0.765957 0.265957 0.000000
+310 -0.745195 0.286345 0.000000
+169 -0.744681 0.244681 0.000000
+311 -0.723868 0.265704 0.000000
+168 -0.723404 0.223404 0.000000
+312 -0.702394 0.245107 0.000000
+167 -0.702128 0.202128 0.000000
+313 -0.680976 0.224300 0.000000
+166 -0.680851 0.180851 0.000000
+314 -0.659597 0.203254 0.000000
+165 -0.659574 0.159574 0.000000
+315 -0.638233 0.182030 0.000000
+164 -0.638298 0.138298 0.000000
+316 -0.616915 0.160630 0.000000
+163 -0.617021 0.117021 0.000000
+317 -0.595665 0.139192 0.000000
+162 -0.595745 0.095745 0.000000
+318 -0.574431 0.117764 0.000000
+161 -0.574468 0.074468 0.000000
+319 -0.553204 0.096309 0.000000
+160 -0.553191 0.053191 0.000000
+320 -0.531840 0.075023 0.000000
+159 -0.531915 0.031915 0.000000
+321 -0.510492 0.053733 0.000000
+158 -0.510638 0.010638 0.000000
+322 -0.489127 0.032461 0.000000
+157 -0.489362 -0.010638 0.000000
+323 -0.467791 0.011178 0.000000
+156 -0.468085 -0.031915 0.000000
+324 -0.446536 -0.010172 0.000000
+155 -0.446809 -0.053191 0.000000
+325 -0.425240 -0.031466 0.000000
+154 -0.425532 -0.074468 0.000000
+326 -0.404012 -0.052824 0.000000
+153 -0.404255 -0.095745 0.000000
+327 -0.382783 -0.074175 0.000000
+152 -0.382979 -0.117021 0.000000
+328 -0.361545 -0.095511 0.000000
+151 -0.361702 -0.138298 0.000000
+329 -0.340288 -0.116831 0.000000
+150 -0.340426 -0.159574 0.000000
+330 -0.319034 -0.138159 0.000000
+149 -0.319149 -0.180851 0.000000
+331 -0.297766 -0.159484 0.000000
+148 -0.297872 -0.202128 0.000000
+332 -0.276489 -0.180797 0.000000
+147 -0.276596 -0.223404 0.000000
+333 -0.255087 -0.202105 0.000000
+146 -0.255319 -0.244681 0.000000
+334 -0.233674 -0.223466 0.000000
+145 -0.234043 -0.265957 0.000000
+335 -0.212119 -0.244916 0.000000
+144 -0.212766 -0.287234 0.000000
+336 -0.190657 -0.266392 0.000000
+143 -0.191489 -0.308511 0.000000
+337 -0.169261 -0.287934 0.000000
+142 -0.170213 -0.329787 0.000000
+338 -0.147854 -0.309624 0.000000
+141 -0.148936 -0.351064 0.000000
+339 -0.126431 -0.331364 0.000000
+140 -0.127660 -0.372340 0.000000
+340 -0.105256 -0.353368 0.000000
+139 -0.106383 -0.393617 0.000000
+341 -0.083924 -0.375263 0.000000
+138 -0.085106 -0.414894 0.000000
+342 -0.062097 -0.396824 0.000000
+137 -0.063830 -0.436170 0.000000
+343 -0.040132 -0.418076 0.000000
+136 -0.042553 -0.457447 0.000000
+344 -0.019646 -0.439891 0.000000
+135 -0.021277 -0.478723 0.000000
+345 -0.002475 -0.461393 0.000000
+134 0.000000 -0.500000 0.000000
+346 0.009599 -0.478226 0.000000
+182 0.030000 -0.500000 0.000000
+347 0.029922 -0.474701 0.000000
+183 0.060000 -0.500000 0.000000
+348 0.057397 -0.471689 0.000000
+184 0.090000 -0.500000 0.000000
+349 0.087960 -0.470463 0.000000
+36 1.000000 -0.439394 0.000000
+350 0.969923 -0.439465 0.000000
+37 1.000000 -0.409091 0.000000
+351 0.969894 -0.409159 0.000000
+38 1.000000 -0.378788 0.000000
+352 0.969889 -0.378860 0.000000
+39 1.000000 -0.348485 0.000000
+353 0.969858 -0.348548 0.000000
+40 1.000000 -0.318182 0.000000
+354 0.969843 -0.318237 0.000000
+41 1.000000 -0.287879 0.000000
+355 0.969809 -0.287925 0.000000
+42 1.000000 -0.257576 0.000000
+356 0.969784 -0.257613 0.000000
+43 1.000000 -0.227273 0.000000
+357 0.969757 -0.227302 0.000000
+44 1.000000 -0.196970 0.000000
+358 0.969734 -0.196993 0.000000
+45 1.000000 -0.166667 0.000000
+359 0.969717 -0.166684 0.000000
+46 1.000000 -0.136364 0.000000
+360 0.969707 -0.136377 0.000000
+47 1.000000 -0.106061 0.000000
+361 0.969701 -0.106070 0.000000
+48 1.000000 -0.075758 0.000000
+362 0.969699 -0.075764 0.000000
+49 1.000000 -0.045455 0.000000
+363 0.969698 -0.045459 0.000000
+50 1.000000 -0.015152 0.000000
+364 0.969697 -0.015155 0.000000
+51 1.000000 0.015152 0.000000
+365 0.969697 0.015149 0.000000
+52 1.000000 0.045455 0.000000
+366 0.969697 0.045453 0.000000
+53 1.000000 0.075758 0.000000
+367 0.969697 0.075757 0.000000
+54 1.000000 0.106061 0.000000
+368 0.969697 0.106060 0.000000
+55 1.000000 0.136364 0.000000
+369 0.969697 0.136363 0.000000
+56 1.000000 0.166667 0.000000
+370 0.969697 0.166666 0.000000
+57 1.000000 0.196970 0.000000
+371 0.969697 0.196970 0.000000
+58 1.000000 0.227273 0.000000
+372 0.969697 0.227273 0.000000
+59 1.000000 0.257576 0.000000
+373 0.969795 0.257643 0.000000
+60 1.000000 0.287879 0.000000
+374 0.969874 0.287995 0.000000
+61 1.000000 0.318182 0.000000
+375 0.969958 0.318338 0.000000
+62 1.000000 0.348485 0.000000
+376 0.969973 0.348658 0.000000
+63 1.000000 0.378788 0.000000
+377 0.970000 0.378980 0.000000
+64 1.000000 0.409091 0.000000
+378 0.970010 0.409285 0.000000
+65 1.000000 0.439394 0.000000
+379 0.970062 0.439531 0.000000
+503 -0.841822 0.439272 0.000000
+380 -0.823383 0.418799 0.000000
+381 -0.801366 0.398075 0.000000
+382 -0.781899 0.373291 0.000000
+383 -0.763952 0.349541 0.000000
+384 -0.745019 0.327062 0.000000
+385 -0.724690 0.306382 0.000000
+386 -0.703265 0.286716 0.000000
+387 -0.681523 0.267067 0.000000
+388 -0.660018 0.246574 0.000000
+389 -0.638443 0.225771 0.000000
+390 -0.616833 0.204679 0.000000
+391 -0.595551 0.183139 0.000000
+392 -0.574409 0.161470 0.000000
+393 -0.553222 0.139792 0.000000
+394 -0.531865 0.118362 0.000000
+395 -0.510475 0.096979 0.000000
+396 -0.489009 0.075696 0.000000
+397 -0.467590 0.054392 0.000000
+398 -0.446270 0.033017 0.000000
+399 -0.424906 0.011735 0.000000
+400 -0.403621 -0.009619 0.000000
+401 -0.382400 -0.031023 0.000000
+402 -0.361173 -0.052410 0.000000
+403 -0.339922 -0.073800 0.000000
+404 -0.318715 -0.095237 0.000000
+405 -0.297453 -0.116646 0.000000
+406 -0.276156 -0.138018 0.000000
+407 -0.254723 -0.159496 0.000000
+408 -0.233335 -0.180914 0.000000
+409 -0.211596 -0.202409 0.000000
+410 -0.189882 -0.223912 0.000000
+411 -0.168328 -0.245541 0.000000
+412 -0.146829 -0.267432 0.000000
+413 -0.125347 -0.289413 0.000000
+414 -0.103926 -0.311617 0.000000
+415 -0.082894 -0.334299 0.000000
+416 -0.061472 -0.356952 0.000000
+417 -0.038982 -0.378982 0.000000
+418 -0.016521 -0.400436 0.000000
+419 0.006313 -0.423591 0.000000
+421 0.051537 -0.441884 0.000000
+420 0.021625 -0.449713 0.000000
+422 0.085122 -0.440903 0.000000
+423 0.082947 -0.409815 0.000000
+424 0.086501 -0.379229 0.000000
+425 0.090943 -0.350788 0.000000
+426 0.094089 -0.322730 0.000000
+427 0.095265 -0.294581 0.000000
+428 0.094987 -0.265825 0.000000
+430 0.123468 -0.236299 0.000000
+429 0.093679 -0.236645 0.000000
+431 0.153728 -0.236070 0.000000
+432 0.184970 -0.235916 0.000000
+433 0.216576 -0.236122 0.000000
+434 0.247804 -0.236880 0.000000
+435 0.278352 -0.237921 0.000000
+437 0.308849 -0.268892 0.000000
+436 0.308561 -0.238863 0.000000
+438 0.309174 -0.298672 0.000000
+439 0.309365 -0.327764 0.000000
+440 0.309409 -0.356452 0.000000
+441 0.309369 -0.385004 0.000000
+442 0.309363 -0.413591 0.000000
+443 0.309512 -0.442324 0.000000
+444 0.940150 0.439639 0.000000
+445 0.910293 0.439719 0.000000
+446 0.880439 0.439764 0.000000
+447 0.850606 0.439851 0.000000
+448 0.820773 0.439905 0.000000
+449 0.790942 0.439997 0.000000
+450 0.761110 0.440072 0.000000
+451 0.731277 0.440144 0.000000
+452 0.701441 0.440206 0.000000
+453 0.671603 0.440248 0.000000
+454 0.641769 0.440131 0.000000
+455 0.611935 0.440030 0.000000
+456 0.582116 0.439803 0.000000
+457 0.552311 0.439544 0.000000
+458 0.522496 0.439372 0.000000
+459 0.492687 0.439151 0.000000
+460 0.462850 0.439112 0.000000
+461 0.433008 0.439157 0.000000
+462 0.403152 0.439190 0.000000
+463 0.373296 0.439261 0.000000
+464 0.343430 0.439326 0.000000
+465 0.313568 0.439388 0.000000
+466 0.283709 0.439386 0.000000
+467 0.253856 0.439367 0.000000
+468 0.223995 0.439416 0.000000
+469 0.194134 0.439557 0.000000
+470 0.164274 0.439658 0.000000
+471 0.134410 0.439807 0.000000
+472 0.104545 0.439947 0.000000
+473 0.074682 0.440045 0.000000
+474 0.044821 0.440112 0.000000
+475 0.014962 0.440171 0.000000
+476 -0.014897 0.440222 0.000000
+477 -0.044754 0.440258 0.000000
+478 -0.074665 0.440098 0.000000
+479 -0.104552 0.439944 0.000000
+480 -0.134459 0.439734 0.000000
+481 -0.164338 0.439555 0.000000
+482 -0.194223 0.439404 0.000000
+483 -0.224113 0.439388 0.000000
+484 -0.253960 0.439636 0.000000
+485 -0.283803 0.439999 0.000000
+486 -0.313653 0.440461 0.000000
+487 -0.343481 0.440757 0.000000
+488 -0.373319 0.440824 0.000000
+489 -0.403191 0.440764 0.000000
+490 -0.433152 0.440474 0.000000
+491 -0.463150 0.439907 0.000000
+492 -0.493171 0.439333 0.000000
+493 -0.523226 0.438998 0.000000
+494 -0.553219 0.439109 0.000000
+495 -0.583244 0.439382 0.000000
+496 -0.613191 0.439846 0.000000
+497 -0.642956 0.440951 0.000000
+498 -0.672668 0.441524 0.000000
+499 -0.702061 0.440023 0.000000
+500 -0.733531 0.435747 0.000000
+501 -0.771806 0.426884 0.000000
+553 -0.803350 0.441135 0.000000
+502 -0.822712 0.450765 0.000000
+504 0.939858 -0.439508 0.000000
+505 0.939812 -0.409239 0.000000
+506 0.939773 -0.378933 0.000000
+507 0.939746 -0.348626 0.000000
+508 0.939690 -0.318305 0.000000
+509 0.939650 -0.287984 0.000000
+510 0.939589 -0.257662 0.000000
+511 0.939539 -0.227343 0.000000
+512 0.939492 -0.197025 0.000000
+513 0.939452 -0.166709 0.000000
+514 0.939425 -0.136395 0.000000
+515 0.939409 -0.106084 0.000000
+516 0.939401 -0.075774 0.000000
+517 0.939397 -0.045467 0.000000
+518 0.939395 -0.015160 0.000000
+519 0.939394 0.015146 0.000000
+520 0.939394 0.045451 0.000000
+521 0.939394 0.075755 0.000000
+522 0.939394 0.106059 0.000000
+523 0.939394 0.136362 0.000000
+524 0.939394 0.166666 0.000000
+525 0.939394 0.196969 0.000000
+526 0.939509 0.227355 0.000000
+527 0.939655 0.257739 0.000000
+528 0.939814 0.288143 0.000000
+529 0.939913 0.318502 0.000000
+530 0.939975 0.348849 0.000000
+531 0.940009 0.379160 0.000000
+532 0.940091 0.409414 0.000000
+533 0.339432 -0.442176 0.000000
+534 0.369451 -0.442035 0.000000
+535 0.399523 -0.441897 0.000000
+536 0.429589 -0.441684 0.000000
+537 0.459678 -0.441399 0.000000
+538 0.489748 -0.441062 0.000000
+539 0.519800 -0.440764 0.000000
+540 0.549844 -0.440522 0.000000
+541 0.579882 -0.440330 0.000000
+542 0.609913 -0.440196 0.000000
+543 0.639937 -0.440110 0.000000
+544 0.669955 -0.440059 0.000000
+545 0.699968 -0.440030 0.000000
+546 0.729978 -0.440014 0.000000
+547 0.759947 -0.439909 0.000000
+548 0.789937 -0.439803 0.000000
+549 0.819882 -0.439709 0.000000
+550 0.849843 -0.439644 0.000000
+551 0.879819 -0.439597 0.000000
+552 0.909827 -0.439560 0.000000
+590 0.040432 -0.405830 0.000000
+670 -0.755455 0.393182 0.000000
+554 -0.740317 0.367025 0.000000
+555 -0.724443 0.344933 0.000000
+556 -0.704864 0.325740 0.000000
+557 -0.683231 0.307949 0.000000
+558 -0.661196 0.289226 0.000000
+559 -0.639331 0.269137 0.000000
+560 -0.617409 0.248586 0.000000
+561 -0.595465 0.227639 0.000000
+562 -0.574245 0.205848 0.000000
+563 -0.553198 0.183785 0.000000
+564 -0.531973 0.162012 0.000000
+565 -0.510563 0.140460 0.000000
+566 -0.489049 0.119070 0.000000
+567 -0.467535 0.097696 0.000000
+568 -0.446086 0.076355 0.000000
+569 -0.424583 0.055151 0.000000
+570 -0.403151 0.033875 0.000000
+571 -0.381891 0.012449 0.000000
+572 -0.360645 -0.008987 0.000000
+573 -0.339410 -0.030497 0.000000
+574 -0.318218 -0.052043 0.000000
+575 -0.296914 -0.073549 0.000000
+576 -0.275621 -0.095064 0.000000
+577 -0.254107 -0.116660 0.000000
+578 -0.232698 -0.138235 0.000000
+579 -0.211016 -0.159840 0.000000
+580 -0.189307 -0.181488 0.000000
+581 -0.167499 -0.203012 0.000000
+582 -0.145658 -0.224717 0.000000
+583 -0.124325 -0.246942 0.000000
+584 -0.102865 -0.269240 0.000000
+585 -0.081284 -0.291521 0.000000
+586 -0.060814 -0.314997 0.000000
+587 -0.039868 -0.338366 0.000000
+588 -0.017866 -0.359895 0.000000
+589 0.007087 -0.381124 0.000000
+591 0.053625 -0.372025 0.000000
+592 0.062247 -0.344115 0.000000
+593 0.067251 -0.318078 0.000000
+594 0.068547 -0.292412 0.000000
+595 0.066596 -0.265244 0.000000
+596 0.064013 -0.237360 0.000000
+598 0.092058 -0.207144 0.000000
+597 0.061821 -0.208936 0.000000
+599 0.122455 -0.205662 0.000000
+600 0.153304 -0.204437 0.000000
+601 0.184778 -0.203622 0.000000
+602 0.216464 -0.203889 0.000000
+603 0.247666 -0.205271 0.000000
+604 0.278273 -0.207077 0.000000
+605 0.308444 -0.208860 0.000000
+607 0.338640 -0.239619 0.000000
+606 0.338501 -0.210215 0.000000
+608 0.338849 -0.269123 0.000000
+609 0.339046 -0.298501 0.000000
+610 0.339149 -0.327425 0.000000
+611 0.339145 -0.356028 0.000000
+612 0.339123 -0.384605 0.000000
+613 0.339265 -0.413346 0.000000
+614 0.910184 0.409525 0.000000
+615 0.880341 0.409627 0.000000
+616 0.850499 0.409699 0.000000
+617 0.820678 0.409821 0.000000
+618 0.790857 0.409913 0.000000
+619 0.761038 0.410043 0.000000
+620 0.731217 0.410148 0.000000
+621 0.701394 0.410247 0.000000
+622 0.671582 0.410190 0.000000
+623 0.641769 0.410064 0.000000
+624 0.611986 0.409820 0.000000
+625 0.582223 0.409499 0.000000
+626 0.552442 0.409166 0.000000
+627 0.522678 0.408817 0.000000
+628 0.492873 0.408592 0.000000
+629 0.463046 0.408514 0.000000
+630 0.433212 0.408428 0.000000
+631 0.403368 0.408486 0.000000
+632 0.373497 0.408591 0.000000
+633 0.343636 0.408704 0.000000
+634 0.313765 0.408808 0.000000
+635 0.283900 0.408804 0.000000
+636 0.254031 0.408855 0.000000
+637 0.224157 0.409018 0.000000
+638 0.194295 0.409133 0.000000
+639 0.164411 0.409379 0.000000
+640 0.134527 0.409617 0.000000
+641 0.104649 0.409788 0.000000
+642 0.074780 0.409914 0.000000
+643 0.044908 0.410046 0.000000
+644 0.015035 0.410160 0.000000
+645 -0.014837 0.410241 0.000000
+646 -0.044786 0.410157 0.000000
+647 -0.074715 0.410004 0.000000
+648 -0.104689 0.409724 0.000000
+649 -0.134629 0.409413 0.000000
+650 -0.164573 0.409102 0.000000
+651 -0.194543 0.408821 0.000000
+652 -0.224425 0.408905 0.000000
+653 -0.254313 0.409349 0.000000
+654 -0.284136 0.410098 0.000000
+655 -0.313927 0.410835 0.000000
+656 -0.343709 0.411217 0.000000
+657 -0.373516 0.411287 0.000000
+658 -0.403429 0.411111 0.000000
+659 -0.433418 0.410580 0.000000
+660 -0.463454 0.409742 0.000000
+661 -0.493588 0.408650 0.000000
+662 -0.523716 0.408015 0.000000
+663 -0.553867 0.408390 0.000000
+664 -0.583894 0.408825 0.000000
+665 -0.613736 0.409400 0.000000
+666 -0.642923 0.411762 0.000000
+667 -0.671337 0.413661 0.000000
+668 -0.699264 0.412614 0.000000
+669 -0.726520 0.405751 0.000000
+671 0.369212 -0.413104 0.000000
+672 0.399261 -0.412881 0.000000
+673 0.429386 -0.412541 0.000000
+674 0.459495 -0.412121 0.000000
+675 0.489590 -0.411630 0.000000
+676 0.519678 -0.411233 0.000000
+677 0.549757 -0.410839 0.000000
+678 0.579820 -0.410560 0.000000
+679 0.609867 -0.410349 0.000000
+680 0.639903 -0.410201 0.000000
+681 0.669930 -0.410109 0.000000
+682 0.699950 -0.410056 0.000000
+683 0.729910 -0.409917 0.000000
+684 0.759884 -0.409777 0.000000
+685 0.789837 -0.409624 0.000000
+686 0.819781 -0.409515 0.000000
+687 0.849736 -0.409429 0.000000
+688 0.879719 -0.409378 0.000000
+689 0.909756 -0.409298 0.000000
+690 0.909695 -0.379025 0.000000
+691 0.909632 -0.348708 0.000000
+692 0.909581 -0.318392 0.000000
+693 0.909498 -0.288059 0.000000
+694 0.909434 -0.257728 0.000000
+695 0.909350 -0.227398 0.000000
+696 0.909280 -0.197070 0.000000
+697 0.909216 -0.166745 0.000000
+698 0.909165 -0.136423 0.000000
+699 0.909131 -0.106105 0.000000
+700 0.909111 -0.075791 0.000000
+701 0.909100 -0.045478 0.000000
+702 0.909095 -0.015168 0.000000
+703 0.909093 0.015140 0.000000
+704 0.909072 0.045446 0.000000
+705 0.909067 0.075752 0.000000
+706 0.909073 0.106056 0.000000
+707 0.909080 0.136361 0.000000
+708 0.909086 0.166664 0.000000
+709 0.909215 0.197066 0.000000
+710 0.909385 0.227464 0.000000
+711 0.909604 0.257898 0.000000
+712 0.909779 0.288313 0.000000
+713 0.909917 0.318693 0.000000
+714 0.909986 0.349031 0.000000
+715 0.910101 0.379296 0.000000
+716 0.024590 -0.357790 0.000000
+717 0.036820 -0.334922 0.000000
+718 0.044018 -0.313436 0.000000
+719 0.043162 -0.290917 0.000000
+720 0.038764 -0.265690 0.000000
+721 0.034906 -0.239103 0.000000
+722 0.031729 -0.211118 0.000000
+724 0.059321 -0.180222 0.000000
+723 0.028990 -0.183087 0.000000
+725 0.089882 -0.177435 0.000000
+726 0.120857 -0.174840 0.000000
+727 0.152193 -0.172458 0.000000
+728 0.184221 -0.171049 0.000000
+729 0.216369 -0.171146 0.000000
+730 0.247776 -0.173286 0.000000
+731 0.278417 -0.176312 0.000000
+732 0.308664 -0.179103 0.000000
+733 0.338547 -0.181101 0.000000
+735 0.368461 -0.210821 0.000000
+734 0.368350 -0.182042 0.000000
+736 0.368647 -0.239838 0.000000
+737 0.368823 -0.269040 0.000000
+738 0.368963 -0.298178 0.000000
+739 0.368962 -0.326977 0.000000
+740 0.368943 -0.355559 0.000000
+741 0.369058 -0.384300 0.000000
+742 0.880206 0.379422 0.000000
+743 0.850374 0.379547 0.000000
+744 0.820545 0.379654 0.000000
+745 0.790741 0.379815 0.000000
+746 0.760936 0.379938 0.000000
+747 0.731131 0.380103 0.000000
+748 0.701348 0.380156 0.000000
+749 0.671569 0.380011 0.000000
+750 0.641834 0.379835 0.000000
+751 0.612110 0.379504 0.000000
+752 0.582403 0.379062 0.000000
+753 0.552705 0.378593 0.000000
+754 0.522968 0.378130 0.000000
+755 0.493191 0.377826 0.000000
+756 0.463421 0.377567 0.000000
+757 0.433611 0.377484 0.000000
+758 0.403752 0.377581 0.000000
+759 0.373895 0.377725 0.000000
+760 0.344005 0.377927 0.000000
+761 0.314132 0.378023 0.000000
+762 0.284270 0.378127 0.000000
+763 0.254375 0.378254 0.000000
+764 0.224504 0.378365 0.000000
+765 0.194586 0.378654 0.000000
+766 0.164677 0.379020 0.000000
+767 0.134778 0.379321 0.000000
+768 0.104899 0.379529 0.000000
+769 0.074998 0.379743 0.000000
+770 0.045096 0.379901 0.000000
+771 0.015198 0.380037 0.000000
+772 -0.014766 0.380029 0.000000
+773 -0.044758 0.379983 0.000000
+774 -0.074799 0.379773 0.000000
+775 -0.104855 0.379434 0.000000
+776 -0.134915 0.378975 0.000000
+777 -0.165034 0.378407 0.000000
+778 -0.195055 0.378036 0.000000
+779 -0.225064 0.378155 0.000000
+780 -0.254915 0.379009 0.000000
+781 -0.284684 0.380269 0.000000
+782 -0.314377 0.381326 0.000000
+783 -0.344063 0.381861 0.000000
+784 -0.373781 0.381860 0.000000
+785 -0.403636 0.381479 0.000000
+786 -0.433541 0.380809 0.000000
+787 -0.463618 0.379545 0.000000
+788 -0.493878 0.377655 0.000000
+789 -0.524338 0.376897 0.000000
+790 -0.554662 0.377591 0.000000
+791 -0.584521 0.378512 0.000000
+792 -0.614901 0.378305 0.000000
+793 -0.644852 0.383150 0.000000
+794 -0.671012 0.387415 0.000000
+795 -0.694363 0.388117 0.000000
+796 -0.716721 0.380258 0.000000
+797 0.399067 -0.383941 0.000000
+798 0.429159 -0.383498 0.000000
+799 0.459314 -0.382900 0.000000
+800 0.489457 -0.382364 0.000000
+801 0.519590 -0.381760 0.000000
+802 0.549686 -0.381311 0.000000
+803 0.579760 -0.380889 0.000000
+804 0.609819 -0.380588 0.000000
+805 0.639864 -0.380360 0.000000
+806 0.669900 -0.380203 0.000000
+807 0.699859 -0.379982 0.000000
+808 0.729822 -0.379783 0.000000
+809 0.759760 -0.379578 0.000000
+810 0.789702 -0.379418 0.000000
+811 0.819647 -0.379286 0.000000
+812 0.849617 -0.379205 0.000000
+813 0.879641 -0.379103 0.000000
+814 -0.704962 0.359795 0.000000
+815 -0.687415 0.344916 0.000000
+816 -0.665253 0.330013 0.000000
+817 -0.640678 0.312888 0.000000
+818 -0.618867 0.291654 0.000000
+819 -0.596620 0.271672 0.000000
+820 -0.574345 0.250779 0.000000
+821 -0.552849 0.228703 0.000000
+822 -0.532003 0.206108 0.000000
+823 -0.510760 0.184088 0.000000
+824 -0.489162 0.162574 0.000000
+825 -0.467617 0.141083 0.000000
+826 -0.446103 0.119735 0.000000
+827 -0.424459 0.098611 0.000000
+828 -0.402806 0.077505 0.000000
+829 -0.381337 0.056223 0.000000
+830 -0.359980 0.034804 0.000000
+831 -0.338720 0.013213 0.000000
+832 -0.317525 -0.008472 0.000000
+833 -0.296233 -0.030164 0.000000
+834 -0.274944 -0.051850 0.000000
+835 -0.253374 -0.073554 0.000000
+836 -0.231843 -0.095322 0.000000
+837 -0.210150 -0.117070 0.000000
+838 -0.188532 -0.138903 0.000000
+839 -0.167016 -0.160779 0.000000
+840 -0.144837 -0.182173 0.000000
+841 -0.122682 -0.203880 0.000000
+842 -0.101823 -0.226595 0.000000
+843 -0.081098 -0.248884 0.000000
+844 -0.058277 -0.270876 0.000000
+845 -0.037972 -0.296496 0.000000
+1723 -0.019092 -0.320042 0.000000
+1725 0.000735 -0.340540 0.000000
+848 0.879560 -0.348818 0.000000
+849 0.879472 -0.318487 0.000000
+850 0.879399 -0.288159 0.000000
+851 0.879291 -0.257815 0.000000
+852 0.879205 -0.227473 0.000000
+853 0.879100 -0.197133 0.000000
+854 0.879013 -0.166796 0.000000
+855 0.878936 -0.136464 0.000000
+856 0.878856 -0.106138 0.000000
+857 0.878795 -0.075817 0.000000
+858 0.878746 -0.045502 0.000000
+859 0.878725 -0.015190 0.000000
+860 0.878647 0.015119 0.000000
+861 0.878640 0.045426 0.000000
+862 0.878639 0.075733 0.000000
+863 0.878671 0.106040 0.000000
+864 0.878713 0.136347 0.000000
+865 0.878883 0.166766 0.000000
+866 0.879094 0.197181 0.000000
+867 0.879362 0.227644 0.000000
+868 0.879600 0.258094 0.000000
+869 0.879805 0.288518 0.000000
+870 0.879937 0.318888 0.000000
+871 0.880089 0.349176 0.000000
+846 0.020727 -0.293356 0.000000
+847 0.026838 -0.309682 0.000000
+872 0.012547 -0.270085 0.000000
+873 0.004882 -0.242473 0.000000
+874 0.002208 -0.213381 0.000000
+875 -0.000861 -0.185685 0.000000
+877 0.026244 -0.155284 0.000000
+876 -0.002631 -0.158212 0.000000
+878 0.056250 -0.151549 0.000000
+879 0.087018 -0.147720 0.000000
+880 0.118301 -0.143854 0.000000
+881 0.150274 -0.140350 0.000000
+882 0.183089 -0.137656 0.000000
+883 0.216064 -0.138164 0.000000
+884 0.248380 -0.141045 0.000000
+885 0.279257 -0.145764 0.000000
+886 0.309290 -0.149826 0.000000
+887 0.338917 -0.152404 0.000000
+888 0.368368 -0.153538 0.000000
+890 0.398029 -0.182007 0.000000
+889 0.397803 -0.153563 0.000000
+891 0.398367 -0.210697 0.000000
+892 0.398637 -0.239587 0.000000
+893 0.398861 -0.268621 0.000000
+894 0.398942 -0.297627 0.000000
+895 0.398900 -0.326384 0.000000
+896 0.398951 -0.355111 0.000000
+897 -0.691273 0.370756 0.000000
+898 -0.674493 0.364355 0.000000
+1040 -0.650656 0.354888 0.000000
+976 -0.617833 0.341638 0.000000
+899 -0.598837 0.312875 0.000000
+900 -0.576909 0.294643 0.000000
+901 -0.552428 0.275224 0.000000
+902 -0.530818 0.250993 0.000000
+903 -0.510792 0.227821 0.000000
+904 -0.489498 0.205790 0.000000
+905 -0.467802 0.184485 0.000000
+906 -0.446323 0.163103 0.000000
+907 -0.424713 0.141919 0.000000
+908 -0.402894 0.120971 0.000000
+909 -0.381057 0.100028 0.000000
+910 -0.359274 0.078918 0.000000
+911 -0.337845 0.057399 0.000000
+912 -0.316626 0.035561 0.000000
+913 -0.295418 0.013637 0.000000
+914 -0.274138 -0.008245 0.000000
+915 -0.252589 -0.030100 0.000000
+916 -0.230965 -0.052030 0.000000
+917 -0.209177 -0.074014 0.000000
+918 -0.187347 -0.095982 0.000000
+919 -0.165790 -0.118137 0.000000
+920 -0.144877 -0.140439 0.000000
+921 -0.122656 -0.161118 0.000000
+922 -0.098962 -0.182854 0.000000
+923 -0.079724 -0.207401 0.000000
+924 -0.061150 -0.228771 0.000000
+925 -0.031281 -0.247098 0.000000
+926 -0.013937 -0.279764 0.000000
+927 0.850209 0.349326 0.000000
+928 0.820397 0.349486 0.000000
+929 0.790584 0.349621 0.000000
+930 0.760799 0.349820 0.000000
+931 0.731036 0.349946 0.000000
+932 0.701303 0.349866 0.000000
+933 0.671610 0.349779 0.000000
+934 0.641932 0.349508 0.000000
+935 0.612330 0.349074 0.000000
+936 0.582719 0.348493 0.000000
+937 0.553104 0.347828 0.000000
+938 0.523441 0.347216 0.000000
+939 0.493783 0.346743 0.000000
+940 0.464050 0.346432 0.000000
+941 0.434236 0.346369 0.000000
+942 0.404391 0.346489 0.000000
+943 0.374528 0.346684 0.000000
+944 0.344649 0.346945 0.000000
+945 0.314797 0.347163 0.000000
+946 0.284881 0.347290 0.000000
+947 0.254985 0.347409 0.000000
+948 0.225047 0.347654 0.000000
+949 0.195107 0.348125 0.000000
+950 0.165173 0.348601 0.000000
+951 0.135276 0.348956 0.000000
+952 0.105340 0.349258 0.000000
+953 0.075403 0.349451 0.000000
+954 0.045475 0.349683 0.000000
+955 0.015506 0.349758 0.000000
+956 -0.014542 0.349888 0.000000
+957 -0.044633 0.349863 0.000000
+958 -0.074856 0.349598 0.000000
+959 -0.105110 0.349093 0.000000
+960 -0.135408 0.348361 0.000000
+961 -0.165677 0.347487 0.000000
+962 -0.195905 0.346884 0.000000
+963 -0.226027 0.347092 0.000000
+964 -0.255914 0.348502 0.000000
+965 -0.285596 0.350414 0.000000
+966 -0.315123 0.352026 0.000000
+967 -0.344543 0.352795 0.000000
+968 -0.374040 0.352571 0.000000
+969 -0.403574 0.351874 0.000000
+970 -0.433297 0.351311 0.000000
+971 -0.463227 0.349634 0.000000
+972 -0.493249 0.346733 0.000000
+973 -0.525130 0.344192 0.000000
+974 -0.556069 0.347386 0.000000
+975 -0.584638 0.349299 0.000000
+977 0.429047 -0.354564 0.000000
+978 0.459194 -0.353904 0.000000
+979 0.489384 -0.353164 0.000000
+980 0.519514 -0.352512 0.000000
+981 0.549620 -0.351866 0.000000
+982 0.579698 -0.351370 0.000000
+983 0.609762 -0.350919 0.000000
+984 0.639814 -0.350586 0.000000
+985 0.669775 -0.350217 0.000000
+986 0.699737 -0.349897 0.000000
+987 0.729661 -0.349597 0.000000
+988 0.759593 -0.349367 0.000000
+989 0.789532 -0.349176 0.000000
+990 0.819498 -0.349051 0.000000
+991 0.849515 -0.348917 0.000000
+992 0.849408 -0.318619 0.000000
+993 0.849295 -0.288272 0.000000
+994 0.849198 -0.257931 0.000000
+995 0.849065 -0.227574 0.000000
+996 0.848962 -0.197220 0.000000
+997 0.848836 -0.166870 0.000000
+998 0.848656 -0.136533 0.000000
+999 0.848506 -0.106206 0.000000
+1000 0.848367 -0.075894 0.000000
+1001 0.848296 -0.045579 0.000000
+1002 0.848134 -0.015275 0.000000
+1003 0.848083 0.015034 0.000000
+1004 0.848033 0.045346 0.000000
+1005 0.848096 0.075663 0.000000
+1006 0.848193 0.105984 0.000000
+1007 0.848437 0.136431 0.000000
+1008 0.848728 0.166875 0.000000
+1009 0.849076 0.197378 0.000000
+1010 0.849382 0.227866 0.000000
+1011 0.849654 0.258328 0.000000
+1012 0.849850 0.288735 0.000000
+1013 0.850055 0.319057 0.000000
+1014 -0.026262 -0.214203 0.000000
+1015 -0.029383 -0.187423 0.000000
+1016 -0.032824 -0.158817 0.000000
+1018 -0.003084 -0.132754 0.000000
+1017 -0.029552 -0.130056 0.000000
+1019 0.022484 -0.128711 0.000000
+1020 0.051752 -0.123251 0.000000
+1021 0.082401 -0.118191 0.000000
+1022 0.113973 -0.113014 0.000000
+1023 0.146621 -0.107888 0.000000
+1024 0.180362 -0.103934 0.000000
+1025 0.216332 -0.102481 0.000000
+1026 0.250760 -0.108923 0.000000
+1027 0.282088 -0.115954 0.000000
+1028 0.310848 -0.121366 0.000000
+1029 0.339795 -0.124483 0.000000
+1030 0.368570 -0.125488 0.000000
+1031 0.397434 -0.125259 0.000000
+1033 0.427150 -0.152702 0.000000
+1032 0.426416 -0.124184 0.000000
+1034 0.427780 -0.181299 0.000000
+1035 0.428281 -0.210062 0.000000
+1036 0.428674 -0.238927 0.000000
+1037 0.428954 -0.267922 0.000000
+1038 0.429022 -0.296882 0.000000
+1039 0.429009 -0.325704 0.000000
+1041 0.820193 0.319232 0.000000
+1042 0.790405 0.319425 0.000000
+1043 0.760627 0.319590 0.000000
+1044 0.730930 0.319648 0.000000
+1045 0.701264 0.319642 0.000000
+1046 0.671658 0.319493 0.000000
+1047 0.642145 0.319171 0.000000
+1048 0.612661 0.318582 0.000000
+1049 0.583227 0.317841 0.000000
+1050 0.553746 0.316951 0.000000
+1051 0.524252 0.316219 0.000000
+1052 0.494666 0.315548 0.000000
+1053 0.464954 0.315177 0.000000
+1054 0.435153 0.315100 0.000000
+1055 0.405377 0.315185 0.000000
+1056 0.375483 0.315504 0.000000
+1057 0.345638 0.315877 0.000000
+1058 0.315734 0.316161 0.000000
+1059 0.285817 0.316366 0.000000
+1060 0.255889 0.316515 0.000000
+1061 0.225925 0.316910 0.000000
+1062 0.195940 0.317535 0.000000
+1063 0.166024 0.318176 0.000000
+1064 0.136095 0.318655 0.000000
+1065 0.106116 0.318949 0.000000
+1066 0.076161 0.319226 0.000000
+1067 0.046155 0.319410 0.000000
+1068 0.016077 0.319689 0.000000
+1069 -0.014080 0.319906 0.000000
+1070 -0.044435 0.319837 0.000000
+1071 -0.074887 0.319446 0.000000
+1072 -0.105364 0.318757 0.000000
+1073 -0.135875 0.317703 0.000000
+1074 -0.166440 0.316418 0.000000
+1075 -0.197098 0.315318 0.000000
+1076 -0.227516 0.315492 0.000000
+1077 -0.257562 0.317644 0.000000
+1078 -0.287116 0.320624 0.000000
+1079 -0.316202 0.323175 0.000000
+1080 -0.345132 0.324256 0.000000
+1081 -0.374089 0.323465 0.000000
+1082 -0.403254 0.322262 0.000000
+1083 -0.432294 0.322314 0.000000
+1084 -0.460563 0.320631 0.000000
+1085 -0.490322 0.315176 0.000000
+1086 -0.527019 0.305092 0.000000
+1087 -0.560434 0.318601 0.000000
+1088 -0.581094 0.326893 0.000000
+1089 0.459215 -0.324955 0.000000
+1090 0.489364 -0.324151 0.000000
+1091 0.519476 -0.323335 0.000000
+1092 0.549546 -0.322613 0.000000
+1093 0.579618 -0.321924 0.000000
+1094 0.609675 -0.321348 0.000000
+1095 0.639635 -0.320779 0.000000
+1096 0.669600 -0.320274 0.000000
+1097 0.699520 -0.319810 0.000000
+1098 0.729446 -0.319441 0.000000
+1099 0.759383 -0.319141 0.000000
+1100 0.789349 -0.318940 0.000000
+1101 0.819365 -0.318752 0.000000
+1213 -0.506862 0.272229 0.000000
+1102 -0.488578 0.247543 0.000000
+1103 -0.468378 0.226625 0.000000
+1104 -0.446529 0.206460 0.000000
+1105 -0.425174 0.185182 0.000000
+1106 -0.403551 0.164089 0.000000
+1107 -0.381461 0.143421 0.000000
+1108 -0.359047 0.122895 0.000000
+1109 -0.336992 0.101927 0.000000
+1110 -0.315505 0.080197 0.000000
+1111 -0.294359 0.058038 0.000000
+1112 -0.273215 0.035830 0.000000
+1113 -0.251740 0.013783 0.000000
+1114 -0.230105 -0.008250 0.000000
+1115 -0.208280 -0.030424 0.000000
+1116 -0.186340 -0.052658 0.000000
+1117 -0.164096 -0.074999 0.000000
+1118 -0.142249 -0.098075 0.000000
+1119 -0.124600 -0.120870 0.000000
+1120 -0.101689 -0.138969 0.000000
+1121 -0.070069 -0.159405 0.000000
+1143 -0.056720 -0.191602 0.000000
+1122 -0.048349 -0.211338 0.000000
+1123 0.819230 -0.288441 0.000000
+1124 0.819089 -0.258075 0.000000
+1125 0.818971 -0.227720 0.000000
+1126 0.818800 -0.197352 0.000000
+1127 0.818537 -0.167025 0.000000
+1128 0.818297 -0.136701 0.000000
+1129 0.818075 -0.106416 0.000000
+1130 0.817914 -0.076099 0.000000
+1131 0.817666 -0.045814 0.000000
+1132 0.817482 -0.015504 0.000000
+1133 0.817346 0.014821 0.000000
+1134 0.817349 0.045170 0.000000
+1135 0.817459 0.075531 0.000000
+1136 0.817790 0.106031 0.000000
+1137 0.818215 0.136528 0.000000
+1138 0.818711 0.167083 0.000000
+1139 0.819128 0.197617 0.000000
+1140 0.819470 0.228123 0.000000
+1141 0.819728 0.258571 0.000000
+1142 0.819993 0.288931 0.000000
+1244 -0.059543 -0.125203 0.000000
+1145 -0.022658 -0.104274 0.000000
+1144 -0.050604 -0.097371 0.000000
+1146 -0.003277 -0.113375 0.000000
+1147 0.015569 -0.103577 0.000000
+1148 0.043900 -0.096090 0.000000
+1149 0.074890 -0.089260 0.000000
+1150 0.106866 -0.082653 0.000000
+1151 0.139199 -0.075782 0.000000
+1152 0.174640 -0.068119 0.000000
+1153 0.218123 -0.058538 0.000000
+1154 0.257570 -0.075952 0.000000
+1155 0.287428 -0.087324 0.000000
+1156 0.314488 -0.094635 0.000000
+1157 0.341036 -0.097765 0.000000
+1158 0.368607 -0.097995 0.000000
+1159 0.396531 -0.096956 0.000000
+1160 0.424895 -0.095330 0.000000
+1162 0.455499 -0.122538 0.000000
+1161 0.453743 -0.093457 0.000000
+1163 0.456665 -0.151337 0.000000
+1164 0.457569 -0.180155 0.000000
+1165 0.458233 -0.209057 0.000000
+1166 0.458740 -0.238043 0.000000
+1167 0.459073 -0.267080 0.000000
+1168 0.459174 -0.296028 0.000000
+1169 0.790159 0.289139 0.000000
+1170 0.760439 0.289302 0.000000
+1171 0.730789 0.289384 0.000000
+1172 0.701251 0.289345 0.000000
+1173 0.671804 0.289210 0.000000
+1174 0.642476 0.288781 0.000000
+1175 0.613248 0.288165 0.000000
+1176 0.584018 0.287164 0.000000
+1177 0.554753 0.286148 0.000000
+1178 0.525405 0.285084 0.000000
+1179 0.495870 0.284318 0.000000
+1180 0.466191 0.283786 0.000000
+1181 0.436505 0.283622 0.000000
+1182 0.406661 0.283753 0.000000
+1183 0.376803 0.284156 0.000000
+1184 0.346949 0.284678 0.000000
+1185 0.317067 0.285146 0.000000
+1186 0.287175 0.285413 0.000000
+1187 0.257218 0.285557 0.000000
+1188 0.227216 0.286008 0.000000
+1189 0.197271 0.286938 0.000000
+1190 0.167376 0.287870 0.000000
+1191 0.137356 0.288464 0.000000
+1192 0.107379 0.288752 0.000000
+1193 0.077319 0.288996 0.000000
+1194 0.047198 0.289325 0.000000
+1195 0.016953 0.289780 0.000000
+1196 -0.013452 0.290044 0.000000
+1197 -0.044033 0.289941 0.000000
+1198 -0.074662 0.289458 0.000000
+1199 -0.105330 0.288557 0.000000
+1200 -0.136186 0.287158 0.000000
+1201 -0.167320 0.285224 0.000000
+1202 -0.198610 0.283265 0.000000
+1203 -0.229792 0.283738 0.000000
+1204 -0.260249 0.286423 0.000000
+1205 -0.289549 0.291136 0.000000
+1206 -0.317900 0.295204 0.000000
+1207 -0.345821 0.296788 0.000000
+1208 -0.373271 0.294906 0.000000
+1209 -0.402699 0.292329 0.000000
+1210 -0.431306 0.294539 0.000000
+1211 -0.456225 0.294743 0.000000
+1212 -0.480543 0.286589 0.000000
+1214 0.489324 -0.295158 0.000000
+1215 0.519386 -0.294275 0.000000
+1216 0.549442 -0.293393 0.000000
+1217 0.579473 -0.292561 0.000000
+1218 0.609427 -0.291788 0.000000
+1219 0.639393 -0.291062 0.000000
+1220 0.669310 -0.290383 0.000000
+1221 0.699237 -0.289798 0.000000
+1222 0.729178 -0.289297 0.000000
+1223 0.759153 -0.288937 0.000000
+1224 0.789181 -0.288640 0.000000
+1225 -0.468048 0.263946 0.000000
+1226 -0.449064 0.247078 0.000000
+1227 -0.425672 0.229002 0.000000
+1228 -0.403972 0.207099 0.000000
+1229 -0.382701 0.186041 0.000000
+1230 -0.360028 0.165890 0.000000
+1231 -0.336966 0.146039 0.000000
+1232 -0.314349 0.125401 0.000000
+1233 -0.292863 0.103310 0.000000
+1234 -0.271955 0.080617 0.000000
+1235 -0.250761 0.058067 0.000000
+1236 -0.229184 0.035846 0.000000
+1237 -0.207348 0.013734 0.000000
+1238 -0.185635 -0.008675 0.000000
+1239 -0.163704 -0.031228 0.000000
+1240 -0.140931 -0.053778 0.000000
+1241 -0.113524 -0.078759 0.000000
+1242 -0.108485 -0.104935 0.000000
+1243 -0.087393 -0.114472 0.000000
+1245 0.789012 -0.258320 0.000000
+1246 0.788816 -0.227953 0.000000
+1247 0.788528 -0.197675 0.000000
+1248 0.788213 -0.167365 0.000000
+1249 0.787924 -0.137121 0.000000
+1250 0.787641 -0.106821 0.000000
+1251 0.787323 -0.076569 0.000000
+1252 0.787001 -0.046263 0.000000
+1253 0.786721 -0.015923 0.000000
+1254 0.786547 0.014478 0.000000
+1255 0.786533 0.044918 0.000000
+1256 0.786809 0.075509 0.000000
+1257 0.787410 0.106098 0.000000
+1258 0.788167 0.136724 0.000000
+1259 0.788811 0.167320 0.000000
+1260 0.789306 0.197873 0.000000
+1261 0.789646 0.228368 0.000000
+1262 0.789928 0.258780 0.000000
+1263 -0.079885 -0.088544 0.000000
+1264 -0.068853 -0.063035 0.000000
+1354 -0.096410 -0.050332 0.000000
+1265 -0.040079 -0.071151 0.000000
+1266 -0.003904 -0.078351 0.000000
+1267 0.033276 -0.068635 0.000000
+1268 0.065450 -0.060908 0.000000
+1269 0.097012 -0.053165 0.000000
+1270 0.129040 -0.044661 0.000000
+1271 0.161078 -0.033951 0.000000
+1272 0.194763 -0.019351 0.000000
+1274 0.297376 -0.062036 0.000000
+1273 0.273270 -0.047050 0.000000
+1730 0.318766 -0.071448 0.000000
+1275 0.342006 -0.072798 0.000000
+1276 0.367455 -0.071083 0.000000
+1277 0.394160 -0.068575 0.000000
+1278 0.422042 -0.066164 0.000000
+1279 0.450643 -0.063564 0.000000
+1281 0.482961 -0.091446 0.000000
+1280 0.480092 -0.061554 0.000000
+1282 0.484919 -0.120729 0.000000
+1283 0.486360 -0.149771 0.000000
+1284 0.487432 -0.178843 0.000000
+1285 0.488230 -0.207943 0.000000
+1286 0.488820 -0.237072 0.000000
+1287 0.489125 -0.266135 0.000000
+1288 0.247345 -0.025834 0.000000
+1289 0.224714 0.002070 0.000000
+1290 0.760232 0.258961 0.000000
+1291 0.730702 0.259023 0.000000
+1292 0.701273 0.259041 0.000000
+1293 0.672070 0.258884 0.000000
+1294 0.643041 0.258580 0.000000
+1295 0.614100 0.257855 0.000000
+1296 0.585158 0.256709 0.000000
+1297 0.556141 0.255290 0.000000
+1298 0.526914 0.253967 0.000000
+1299 0.497488 0.252938 0.000000
+1300 0.467901 0.252297 0.000000
+1301 0.438185 0.251985 0.000000
+1302 0.408350 0.252105 0.000000
+1303 0.378486 0.252635 0.000000
+1304 0.348689 0.253392 0.000000
+1305 0.318849 0.254180 0.000000
+1306 0.288993 0.254506 0.000000
+1307 0.259104 0.254448 0.000000
+1308 0.229125 0.254911 0.000000
+1309 0.199160 0.256428 0.000000
+1310 0.169139 0.257831 0.000000
+1311 0.139113 0.258524 0.000000
+1312 0.109056 0.258739 0.000000
+1313 0.078905 0.258871 0.000000
+1314 0.048551 0.259413 0.000000
+1315 0.018071 0.260052 0.000000
+1316 -0.012555 0.260381 0.000000
+1317 -0.043290 0.260233 0.000000
+1318 -0.074004 0.259633 0.000000
+1319 -0.104865 0.258556 0.000000
+1320 -0.136045 0.256823 0.000000
+1321 -0.167724 0.254059 0.000000
+1322 -0.199349 0.251178 0.000000
+1323 -0.232977 0.249469 0.000000
+1324 -0.265381 0.255203 0.000000
+1325 -0.294998 0.262571 0.000000
+1326 -0.321525 0.268904 0.000000
+1327 -0.345370 0.271940 0.000000
+1328 -0.370494 0.266937 0.000000
+1329 -0.403103 0.256984 0.000000
+1330 -0.434068 0.268653 0.000000
+1331 0.519208 -0.265168 0.000000
+1332 0.549223 -0.264174 0.000000
+1333 0.579169 -0.263217 0.000000
+1334 0.609116 -0.262295 0.000000
+1335 0.639018 -0.261389 0.000000
+1336 0.668940 -0.260559 0.000000
+1337 0.698881 -0.259811 0.000000
+1338 0.728862 -0.259217 0.000000
+1339 0.758922 -0.258684 0.000000
+1340 -0.452781 0.275781 0.000000
+1341 -0.382340 0.228180 0.000000
+1342 -0.361422 0.207003 0.000000
+1343 -0.338786 0.188092 0.000000
+1344 -0.314784 0.169671 0.000000
+1345 -0.291764 0.149035 0.000000
+1346 -0.269929 0.126678 0.000000
+1347 -0.249214 0.103062 0.000000
+1348 -0.228144 0.080017 0.000000
+1349 -0.206339 0.057748 0.000000
+1350 -0.184396 0.035681 0.000000
+1351 -0.163152 0.013135 0.000000
+1352 -0.142212 -0.009361 0.000000
+1353 -0.120583 -0.030889 0.000000
+1355 0.758631 -0.228471 0.000000
+1356 0.758281 -0.198203 0.000000
+1357 0.757927 -0.168007 0.000000
+1358 0.757568 -0.137774 0.000000
+1359 0.757184 -0.107553 0.000000
+1360 0.756766 -0.077283 0.000000
+1361 0.756298 -0.046951 0.000000
+1362 0.755911 -0.016494 0.000000
+1363 0.755578 0.014071 0.000000
+1364 0.755530 0.044794 0.000000
+1365 0.756009 0.075565 0.000000
+1366 0.757080 0.106305 0.000000
+1367 0.758211 0.136989 0.000000
+1368 0.759083 0.167591 0.000000
+1369 0.759636 0.198105 0.000000
+1370 0.760021 0.228532 0.000000
+1371 -0.360079 0.243309 0.000000
+1372 -0.341982 0.226476 0.000000
+1373 -0.319699 0.210923 0.000000
+1374 -0.294047 0.194220 0.000000
+1375 -0.267776 0.174727 0.000000
+1376 -0.245955 0.149439 0.000000
+1377 -0.226125 0.124759 0.000000
+1378 -0.205333 0.101062 0.000000
+1379 -0.183133 0.079096 0.000000
+1380 -0.161064 0.056981 0.000000
+1381 -0.140320 0.034408 0.000000
+1382 -0.120789 0.011646 0.000000
+1383 -0.101307 -0.009581 0.000000
+1384 -0.080345 -0.027959 0.000000
+1385 -0.059504 -0.040679 0.000000
+1386 -0.035185 -0.045016 0.000000
+1387 -0.006250 -0.045396 0.000000
+1388 0.025205 -0.040029 0.000000
+1389 0.055491 -0.032869 0.000000
+1390 0.085590 -0.024698 0.000000
+1391 0.115671 -0.015090 0.000000
+1392 0.145238 -0.002575 0.000000
+1393 0.175423 0.012389 0.000000
+1394 0.204047 0.030821 0.000000
+1396 0.252433 0.023671 0.000000
+1395 0.231976 0.050390 0.000000
+1397 0.272784 -0.001727 0.000000
+1398 0.292997 -0.024140 0.000000
+1399 -0.344301 0.253627 0.000000
+1476 -0.327649 0.245851 0.000000
+1475 -0.303967 0.235261 0.000000
+1474 -0.275606 0.222900 0.000000
+1473 -0.238652 0.207257 0.000000
+1400 -0.218940 0.171385 0.000000
+1401 -0.201026 0.143814 0.000000
+1402 -0.182515 0.120273 0.000000
+1403 -0.160481 0.099745 0.000000
+1404 -0.136284 0.078764 0.000000
+1405 -0.115924 0.053849 0.000000
+1406 -0.097288 0.029482 0.000000
+1407 -0.080305 0.007789 0.000000
+1408 -0.064615 -0.010938 0.000000
+1409 -0.053222 -0.024526 0.000000
+1410 -0.035546 -0.021961 0.000000
+1411 -0.011462 -0.018372 0.000000
+1412 0.015936 -0.012857 0.000000
+1413 0.044720 -0.005800 0.000000
+1414 0.073013 0.002618 0.000000
+1415 0.101137 0.013052 0.000000
+1416 0.129049 0.026082 0.000000
+1417 0.156365 0.041756 0.000000
+1418 0.184159 0.058506 0.000000
+1419 0.210898 0.077145 0.000000
+1421 0.258587 0.071409 0.000000
+1420 0.235912 0.097537 0.000000
+1422 0.279603 0.044767 0.000000
+1423 0.299087 0.018547 0.000000
+1424 0.315969 -0.006303 0.000000
+1731 0.310464 -0.042060 0.000000
+1425 0.330327 -0.029348 0.000000
+1428 0.363244 -0.044857 0.000000
+1427 0.340084 -0.050081 0.000000
+1429 0.389469 -0.040337 0.000000
+1430 0.417204 -0.036504 0.000000
+1431 0.446200 -0.033641 0.000000
+1432 0.475860 -0.030939 0.000000
+1434 0.509909 -0.059392 0.000000
+1433 0.506170 -0.028724 0.000000
+1435 0.512650 -0.089356 0.000000
+1436 0.514653 -0.118826 0.000000
+1437 0.516198 -0.148172 0.000000
+1438 0.517372 -0.177496 0.000000
+1439 0.518251 -0.206797 0.000000
+1440 0.518809 -0.236000 0.000000
+1441 0.730529 0.228662 0.000000
+1442 0.701323 0.228754 0.000000
+1443 0.672423 0.228800 0.000000
+1444 0.643741 0.228640 0.000000
+1445 0.615194 0.227932 0.000000
+1446 0.586654 0.226401 0.000000
+1447 0.557941 0.224612 0.000000
+1448 0.528929 0.223020 0.000000
+1449 0.499611 0.221647 0.000000
+1450 0.470054 0.220728 0.000000
+1451 0.440302 0.220268 0.000000
+1452 0.410353 0.220340 0.000000
+1453 0.380500 0.220966 0.000000
+1454 0.350669 0.222138 0.000000
+1455 0.320983 0.223466 0.000000
+1456 0.291443 0.223868 0.000000
+1457 0.261694 0.223241 0.000000
+1458 0.231496 0.223560 0.000000
+1459 0.201231 0.226194 0.000000
+1460 0.171137 0.228421 0.000000
+1461 0.141193 0.229123 0.000000
+1462 0.111175 0.228865 0.000000
+1463 0.080816 0.228761 0.000000
+1464 0.050208 0.229564 0.000000
+1465 0.019472 0.230514 0.000000
+1466 -0.011343 0.230954 0.000000
+1467 -0.042117 0.230757 0.000000
+1468 -0.072814 0.229991 0.000000
+1469 -0.103589 0.228826 0.000000
+1470 -0.134529 0.227009 0.000000
+1471 -0.165046 0.223412 0.000000
+1472 -0.197706 0.217242 0.000000
+1477 0.548816 -0.234928 0.000000
+1478 0.578757 -0.233868 0.000000
+1479 0.608646 -0.232776 0.000000
+1480 0.638544 -0.231722 0.000000
+1481 0.668475 -0.230745 0.000000
+1482 0.698449 -0.229883 0.000000
+1483 0.728502 -0.229106 0.000000
+1484 0.728097 -0.198990 0.000000
+1485 0.727706 -0.168870 0.000000
+1486 0.727268 -0.138716 0.000000
+1487 0.726781 -0.108502 0.000000
+1488 0.726187 -0.078232 0.000000
+1489 0.725568 -0.047739 0.000000
+1490 0.724892 -0.017056 0.000000
+1491 0.724335 0.013822 0.000000
+1492 0.724164 0.044888 0.000000
+1493 0.724932 0.075908 0.000000
+1494 0.726778 0.106759 0.000000
+1495 0.728488 0.137429 0.000000
+1496 0.729633 0.167949 0.000000
+1497 0.730244 0.198323 0.000000
+1624 -0.176223 0.159072 0.000000
+1565 -0.188636 0.185770 0.000000
+1498 -0.161949 0.137398 0.000000
+1499 -0.140110 0.121576 0.000000
+1500 -0.106545 0.104454 0.000000
+1501 -0.087526 0.070737 0.000000
+1502 -0.072361 0.044417 0.000000
+1503 -0.057196 0.020338 0.000000
+1504 -0.044389 -0.001210 0.000000
+1505 -0.020536 0.005897 0.000000
+1506 0.005896 0.012768 0.000000
+1507 0.032832 0.019947 0.000000
+1508 0.059602 0.028562 0.000000
+1509 0.085773 0.039645 0.000000
+1510 0.111268 0.053907 0.000000
+1511 0.136949 0.069672 0.000000
+1512 0.163708 0.085560 0.000000
+1513 0.189111 0.102843 0.000000
+1514 0.212292 0.123719 0.000000
+1516 0.261207 0.118056 0.000000
+1515 0.235876 0.152123 0.000000
+1517 0.283993 0.091978 0.000000
+1518 0.308470 0.065373 0.000000
+1519 0.327160 0.035468 0.000000
+1520 0.342353 0.007735 0.000000
+1522 0.381592 -0.011628 0.000000
+1521 0.354627 -0.019214 0.000000
+1523 0.410351 -0.006674 0.000000
+1524 0.440111 -0.003133 0.000000
+1525 0.470506 -0.000324 0.000000
+1526 0.501426 0.002030 0.000000
+1528 0.536865 -0.026396 0.000000
+1527 0.532740 0.004491 0.000000
+1529 0.540168 -0.057063 0.000000
+1530 0.542678 -0.087141 0.000000
+1531 0.544636 -0.116906 0.000000
+1532 0.546173 -0.146556 0.000000
+1533 0.547355 -0.176132 0.000000
+1534 0.548177 -0.205590 0.000000
+1535 0.701238 0.198603 0.000000
+1536 0.672622 0.198890 0.000000
+1537 0.644404 0.199062 0.000000
+1538 0.616485 0.198422 0.000000
+1539 0.588583 0.196609 0.000000
+1540 0.560344 0.194273 0.000000
+1541 0.531588 0.192198 0.000000
+1542 0.502343 0.190448 0.000000
+1543 0.472715 0.189284 0.000000
+1544 0.442689 0.188639 0.000000
+1545 0.412525 0.188524 0.000000
+1546 0.382371 0.189211 0.000000
+1547 0.352532 0.190926 0.000000
+1548 0.323383 0.193301 0.000000
+1549 0.294480 0.194119 0.000000
+1550 0.265376 0.192145 0.000000
+1551 0.233926 0.191285 0.000000
+1552 0.201869 0.196460 0.000000
+1553 0.172943 0.200159 0.000000
+1554 0.143553 0.200482 0.000000
+1555 0.113782 0.199148 0.000000
+1556 0.083231 0.198393 0.000000
+1557 0.052135 0.199830 0.000000
+1558 0.021092 0.201266 0.000000
+1559 -0.009922 0.201866 0.000000
+1560 -0.040700 0.201511 0.000000
+1561 -0.071228 0.200508 0.000000
+1562 -0.101438 0.199427 0.000000
+1563 -0.131223 0.198305 0.000000
+1564 -0.159535 0.194395 0.000000
+1566 0.578134 -0.204375 0.000000
+1567 0.608049 -0.203173 0.000000
+1568 0.637979 -0.201995 0.000000
+1569 0.667952 -0.200914 0.000000
+1570 0.698007 -0.199916 0.000000
+1571 0.697532 -0.169905 0.000000
+1572 0.697003 -0.139791 0.000000
+1573 0.696380 -0.109617 0.000000
+1574 0.695597 -0.079155 0.000000
+1575 0.694688 -0.048469 0.000000
+1576 0.693667 -0.017488 0.000000
+1577 0.692727 0.013854 0.000000
+1578 0.692318 0.045364 0.000000
+1579 0.693432 0.076853 0.000000
+1580 0.696657 0.107710 0.000000
+1581 0.699224 0.138222 0.000000
+1582 0.700642 0.168450 0.000000
+1583 0.371214 0.017358 0.000000
+1584 0.401698 0.024012 0.000000
+1585 0.432909 0.027812 0.000000
+1586 0.464397 0.030558 0.000000
+1587 0.495928 0.033147 0.000000
+1588 0.527675 0.035870 0.000000
+1590 0.564245 0.007571 0.000000
+1589 0.559782 0.038972 0.000000
+1591 0.567852 -0.023731 0.000000
+1592 0.570714 -0.054539 0.000000
+1593 0.572988 -0.084893 0.000000
+1594 0.574798 -0.114991 0.000000
+1595 0.576239 -0.144944 0.000000
+1596 0.577286 -0.174743 0.000000
+1597 0.672269 0.169162 0.000000
+1598 0.644670 0.169920 0.000000
+1599 0.617885 0.169958 0.000000
+1600 0.591267 0.167167 0.000000
+1601 0.563698 0.164256 0.000000
+1602 0.535150 0.161558 0.000000
+1603 0.505769 0.159594 0.000000
+1604 0.475777 0.158082 0.000000
+1605 0.445365 0.156971 0.000000
+1606 0.414493 0.156592 0.000000
+1607 0.383621 0.157125 0.000000
+1608 0.353035 0.160313 0.000000
+1609 0.324219 0.164732 0.000000
+1610 0.298724 0.166543 0.000000
+1611 0.271789 0.160670 0.000000
+1651 0.199126 0.167293 0.000000
+1612 0.173017 0.174362 0.000000
+1613 0.147220 0.172881 0.000000
+1614 0.118370 0.169386 0.000000
+1615 0.086388 0.167114 0.000000
+1616 0.053408 0.170284 0.000000
+1617 0.022505 0.172511 0.000000
+1618 -0.008623 0.173287 0.000000
+1619 -0.039412 0.172568 0.000000
+1620 -0.069931 0.170982 0.000000
+1621 -0.100028 0.170519 0.000000
+1622 -0.128149 0.170934 0.000000
+1623 -0.152764 0.169059 0.000000
+1625 0.607304 -0.173381 0.000000
+1626 0.637335 -0.172134 0.000000
+1627 0.667402 -0.170981 0.000000
+1628 0.666783 -0.140930 0.000000
+1629 0.665958 -0.110595 0.000000
+1630 0.664919 -0.080023 0.000000
+1631 0.663654 -0.049157 0.000000
+1632 0.662211 -0.017858 0.000000
+1633 0.660734 0.013978 0.000000
+1634 0.659723 0.046364 0.000000
+1635 0.661116 0.078961 0.000000
+1636 0.667246 0.110337 0.000000
+1637 0.670745 0.139548 0.000000
+1638 -0.147799 0.150716 0.000000
+1674 -0.128684 0.145691 0.000000
+1673 -0.100884 0.140613 0.000000
+1640 -0.044570 0.054253 0.000000
+1639 -0.056446 0.080725 0.000000
+1641 -0.032234 0.029928 0.000000
+1642 -0.006024 0.037236 0.000000
+1643 0.020359 0.044478 0.000000
+1644 0.045874 0.052989 0.000000
+1645 0.069878 0.064328 0.000000
+1646 0.092192 0.079845 0.000000
+1647 0.116360 0.097485 0.000000
+1648 0.143869 0.111631 0.000000
+1649 0.168436 0.126174 0.000000
+1650 0.188283 0.142853 0.000000
+1652 0.358550 0.049279 0.000000
+1653 0.392847 0.055605 0.000000
+1654 0.425498 0.059383 0.000000
+1655 0.457980 0.061988 0.000000
+1656 0.489964 0.064583 0.000000
+1657 0.521756 0.067277 0.000000
+1658 0.553712 0.070546 0.000000
+1660 0.592058 0.043198 0.000000
+1659 0.586126 0.075499 0.000000
+1661 0.596024 0.010863 0.000000
+1662 0.599068 -0.020945 0.000000
+1663 0.601524 -0.052102 0.000000
+1664 0.603505 -0.082800 0.000000
+1665 0.605107 -0.113223 0.000000
+1666 0.606324 -0.143417 0.000000
+1667 0.171493 0.155016 0.000000
+1696 0.153180 0.147834 0.000000
+1695 0.126680 0.139244 0.000000
+1694 0.090766 0.129758 0.000000
+1668 0.053205 0.140969 0.000000
+1669 0.022718 0.144596 0.000000
+1670 -0.007526 0.145267 0.000000
+1671 -0.037468 0.143511 0.000000
+1672 -0.068877 0.141356 0.000000
+1675 0.643744 0.141837 0.000000
+1676 0.619340 0.143192 0.000000
+1677 0.595381 0.139507 0.000000
+1678 0.568637 0.134721 0.000000
+1679 0.539884 0.131489 0.000000
+1680 0.510058 0.128996 0.000000
+1681 0.479577 0.126939 0.000000
+1682 0.448267 0.125455 0.000000
+1683 0.416247 0.124554 0.000000
+1684 0.383738 0.125001 0.000000
+1685 0.350239 0.128717 0.000000
+1686 0.321325 0.138246 0.000000
+1687 0.301817 0.146183 0.000000
+1688 0.285401 0.133759 0.000000
+1689 0.636514 -0.142039 0.000000
+1719 0.635429 -0.111764 0.000000
+1718 0.634140 -0.081198 0.000000
+1717 0.632516 -0.050293 0.000000
+1716 0.630582 -0.018910 0.000000
+1715 0.628349 0.013302 0.000000
+1714 0.625872 0.047113 0.000000
+1713 0.623071 0.085699 0.000000
+1690 0.639907 0.114926 0.000000
+1707 0.307163 0.113588 0.000000
+1706 0.342214 0.089227 0.000000
+1711 -0.017903 0.061881 0.000000
+1710 0.008821 0.068462 0.000000
+1691 0.033317 0.075470 0.000000
+1692 0.053991 0.085182 0.000000
+1693 0.071183 0.102698 0.000000
+1697 0.620610 0.123584 0.000000
+1698 0.602483 0.113360 0.000000
+1699 0.575576 0.106296 0.000000
+1700 0.546046 0.101572 0.000000
+1701 0.515442 0.098403 0.000000
+1702 0.484200 0.095912 0.000000
+1703 0.452238 0.093868 0.000000
+1704 0.419736 0.091956 0.000000
+1705 0.384875 0.090284 0.000000
+1708 0.047562 0.116186 0.000000
+1709 0.022481 0.118825 0.000000
+1727 -0.004797 0.117886 0.000000
+1729 -0.034224 0.114874 0.000000
+1712 -0.065866 0.110365 0.000000
+1721 0.041866 0.099182 0.000000
+1720 0.024047 0.096178 0.000000
+1722 0.000778 -0.304198 0.000000
+1724 0.015932 -0.322856 0.000000
+1726 -0.000448 0.092548 0.000000
+1728 -0.027125 0.087460 0.000000
+1426 0.323105 -0.054530 0.000000
+1 1 quad 195 133 67 180
+2 1 quad 196 132 133 195
+3 1 quad 196 195 180 179
+4 1 quad 132 196 197 131
+5 1 quad 196 179 178 197
+6 1 quad 181 187 198 185
+7 1 quad 187 188 199 198
+8 1 quad 188 189 200 199
+9 1 quad 189 190 201 200
+10 1 quad 190 191 202 201
+11 1 quad 191 192 203 202
+12 1 quad 192 186 204 203
+13 1 quad 186 206 205 204
+14 1 quad 186 193 207 206
+15 1 quad 193 194 208 207
+16 1 quad 194 1 209 208
+17 1 quad 1 211 210 209
+18 1 quad 1 3 212 211
+19 1 quad 3 4 213 212
+20 1 quad 4 5 214 213
+21 1 quad 5 6 215 214
+22 1 quad 6 7 216 215
+23 1 quad 7 8 217 216
+24 1 quad 8 2 10 217
+25 1 quad 34 68 218 66
+26 1 quad 68 69 219 218
+27 1 quad 69 70 220 219
+28 1 quad 70 71 221 220
+29 1 quad 71 72 222 221
+30 1 quad 72 73 223 222
+31 1 quad 73 74 224 223
+32 1 quad 74 75 225 224
+33 1 quad 75 76 226 225
+34 1 quad 76 77 227 226
+35 1 quad 77 78 228 227
+36 1 quad 78 79 229 228
+37 1 quad 79 80 230 229
+38 1 quad 80 81 231 230
+39 1 quad 81 82 232 231
+40 1 quad 82 83 233 232
+41 1 quad 83 84 234 233
+42 1 quad 84 85 235 234
+43 1 quad 85 86 236 235
+44 1 quad 86 87 237 236
+45 1 quad 87 88 238 237
+46 1 quad 88 89 239 238
+47 1 quad 89 90 240 239
+48 1 quad 90 91 241 240
+49 1 quad 91 92 242 241
+50 1 quad 92 93 243 242
+51 1 quad 93 94 244 243
+52 1 quad 94 95 245 244
+53 1 quad 95 96 246 245
+54 1 quad 96 97 247 246
+55 1 quad 97 98 248 247
+56 1 quad 98 99 249 248
+57 1 quad 99 100 250 249
+58 1 quad 100 101 251 250
+59 1 quad 101 102 252 251
+60 1 quad 102 103 253 252
+61 1 quad 103 104 254 253
+62 1 quad 104 105 255 254
+63 1 quad 105 106 256 255
+64 1 quad 106 107 257 256
+65 1 quad 107 108 258 257
+66 1 quad 108 109 259 258
+67 1 quad 109 110 260 259
+68 1 quad 110 111 261 260
+69 1 quad 111 112 262 261
+70 1 quad 112 113 263 262
+71 1 quad 113 114 264 263
+72 1 quad 114 115 265 264
+73 1 quad 115 116 266 265
+74 1 quad 116 117 267 266
+75 1 quad 117 118 268 267
+76 1 quad 118 119 269 268
+77 1 quad 119 120 270 269
+78 1 quad 120 121 271 270
+79 1 quad 121 122 272 271
+80 1 quad 122 123 273 272
+81 1 quad 123 124 274 273
+82 1 quad 124 125 275 274
+83 1 quad 125 126 276 275
+84 1 quad 126 127 277 276
+85 1 quad 127 128 278 277
+86 1 quad 128 129 279 278
+87 1 quad 129 130 303 279
+88 1 quad 130 131 197 303
+89 1 quad 10 11 280 217
+90 1 quad 11 12 281 280
+91 1 quad 12 13 282 281
+92 1 quad 13 14 283 282
+93 1 quad 14 15 284 283
+94 1 quad 15 16 285 284
+95 1 quad 16 17 286 285
+96 1 quad 17 18 287 286
+97 1 quad 18 19 288 287
+98 1 quad 19 20 289 288
+99 1 quad 20 21 290 289
+100 1 quad 21 22 291 290
+101 1 quad 22 23 292 291
+102 1 quad 23 24 293 292
+103 1 quad 24 25 294 293
+104 1 quad 25 26 295 294
+105 1 quad 26 27 296 295
+106 1 quad 27 28 297 296
+107 1 quad 28 29 298 297
+108 1 quad 29 30 299 298
+109 1 quad 30 31 300 299
+110 1 quad 31 32 301 300
+111 1 quad 32 33 302 301
+112 1 quad 33 9 35 302
+113 1 quad 178 177 303 197
+114 1 quad 177 176 304 303
+115 1 quad 176 175 305 304
+116 1 quad 175 174 306 305
+117 1 quad 174 173 307 306
+118 1 quad 173 172 308 307
+119 1 quad 172 171 309 308
+120 1 quad 171 170 310 309
+121 1 quad 170 169 311 310
+122 1 quad 169 168 312 311
+123 1 quad 168 167 313 312
+124 1 quad 167 166 314 313
+125 1 quad 166 165 315 314
+126 1 quad 165 164 316 315
+127 1 quad 164 163 317 316
+128 1 quad 163 162 318 317
+129 1 quad 162 161 319 318
+130 1 quad 161 160 320 319
+131 1 quad 160 159 321 320
+132 1 quad 159 158 322 321
+133 1 quad 158 157 323 322
+134 1 quad 157 156 324 323
+135 1 quad 156 155 325 324
+136 1 quad 155 154 326 325
+137 1 quad 154 153 327 326
+138 1 quad 153 152 328 327
+139 1 quad 152 151 329 328
+140 1 quad 151 150 330 329
+141 1 quad 150 149 331 330
+142 1 quad 149 148 332 331
+143 1 quad 148 147 333 332
+144 1 quad 147 146 334 333
+145 1 quad 146 145 335 334
+146 1 quad 145 144 336 335
+147 1 quad 144 143 337 336
+148 1 quad 143 142 338 337
+149 1 quad 142 141 339 338
+150 1 quad 141 140 340 339
+151 1 quad 140 139 341 340
+152 1 quad 139 138 342 341
+153 1 quad 138 137 343 342
+154 1 quad 137 136 344 343
+155 1 quad 136 135 345 344
+156 1 quad 135 134 346 345
+157 1 quad 134 182 347 346
+158 1 quad 182 183 348 347
+159 1 quad 183 184 349 348
+160 1 quad 184 185 198 349
+161 1 quad 35 36 350 302
+162 1 quad 36 37 351 350
+163 1 quad 37 38 352 351
+164 1 quad 38 39 353 352
+165 1 quad 39 40 354 353
+166 1 quad 40 41 355 354
+167 1 quad 41 42 356 355
+168 1 quad 42 43 357 356
+169 1 quad 43 44 358 357
+170 1 quad 44 45 359 358
+171 1 quad 45 46 360 359
+172 1 quad 46 47 361 360
+173 1 quad 47 48 362 361
+174 1 quad 48 49 363 362
+175 1 quad 49 50 364 363
+176 1 quad 50 51 365 364
+177 1 quad 51 52 366 365
+178 1 quad 52 53 367 366
+179 1 quad 53 54 368 367
+180 1 quad 54 55 369 368
+181 1 quad 55 56 370 369
+182 1 quad 56 57 371 370
+183 1 quad 57 58 372 371
+184 1 quad 58 59 373 372
+185 1 quad 59 60 374 373
+186 1 quad 60 61 375 374
+187 1 quad 61 62 376 375
+188 1 quad 62 63 377 376
+189 1 quad 63 64 378 377
+190 1 quad 64 65 379 378
+191 1 quad 65 66 218 379
+192 1 quad 303 304 503 279
+193 1 quad 304 305 380 503
+194 1 quad 305 306 381 380
+195 1 quad 306 307 382 381
+196 1 quad 307 308 383 382
+197 1 quad 308 309 384 383
+198 1 quad 309 310 385 384
+199 1 quad 310 311 386 385
+200 1 quad 311 312 387 386
+201 1 quad 312 313 388 387
+202 1 quad 313 314 389 388
+203 1 quad 314 315 390 389
+204 1 quad 315 316 391 390
+205 1 quad 316 317 392 391
+206 1 quad 317 318 393 392
+207 1 quad 318 319 394 393
+208 1 quad 319 320 395 394
+209 1 quad 320 321 396 395
+210 1 quad 321 322 397 396
+211 1 quad 322 323 398 397
+212 1 quad 323 324 399 398
+213 1 quad 324 325 400 399
+214 1 quad 325 326 401 400
+215 1 quad 326 327 402 401
+216 1 quad 327 328 403 402
+217 1 quad 328 329 404 403
+218 1 quad 329 330 405 404
+219 1 quad 330 331 406 405
+220 1 quad 331 332 407 406
+221 1 quad 332 333 408 407
+222 1 quad 333 334 409 408
+223 1 quad 334 335 410 409
+224 1 quad 335 336 411 410
+225 1 quad 336 337 412 411
+226 1 quad 337 338 413 412
+227 1 quad 338 339 414 413
+228 1 quad 339 340 415 414
+229 1 quad 340 341 416 415
+230 1 quad 341 342 417 416
+231 1 quad 342 343 418 417
+232 1 quad 343 344 419 418
+233 1 quad 347 348 421 420
+234 1 quad 348 349 422 421
+235 1 quad 349 198 199 422
+236 1 quad 199 200 423 422
+237 1 quad 200 201 424 423
+238 1 quad 201 202 425 424
+239 1 quad 202 203 426 425
+240 1 quad 203 204 427 426
+241 1 quad 204 205 428 427
+242 1 quad 205 430 429 428
+243 1 quad 205 206 431 430
+244 1 quad 206 207 432 431
+245 1 quad 207 208 433 432
+246 1 quad 208 209 434 433
+247 1 quad 209 210 435 434
+248 1 quad 210 437 436 435
+249 1 quad 210 211 438 437
+250 1 quad 211 212 439 438
+251 1 quad 212 213 440 439
+252 1 quad 213 214 441 440
+253 1 quad 214 215 442 441
+254 1 quad 215 216 443 442
+255 1 quad 216 217 280 443
+256 1 quad 218 219 444 379
+257 1 quad 219 220 445 444
+258 1 quad 220 221 446 445
+259 1 quad 221 222 447 446
+260 1 quad 222 223 448 447
+261 1 quad 223 224 449 448
+262 1 quad 224 225 450 449
+263 1 quad 225 226 451 450
+264 1 quad 226 227 452 451
+265 1 quad 227 228 453 452
+266 1 quad 228 229 454 453
+267 1 quad 229 230 455 454
+268 1 quad 230 231 456 455
+269 1 quad 231 232 457 456
+270 1 quad 232 233 458 457
+271 1 quad 233 234 459 458
+272 1 quad 234 235 460 459
+273 1 quad 235 236 461 460
+274 1 quad 236 237 462 461
+275 1 quad 237 238 463 462
+276 1 quad 238 239 464 463
+277 1 quad 239 240 465 464
+278 1 quad 240 241 466 465
+279 1 quad 241 242 467 466
+280 1 quad 242 243 468 467
+281 1 quad 243 244 469 468
+282 1 quad 244 245 470 469
+283 1 quad 245 246 471 470
+284 1 quad 246 247 472 471
+285 1 quad 247 248 473 472
+286 1 quad 248 249 474 473
+287 1 quad 249 250 475 474
+288 1 quad 250 251 476 475
+289 1 quad 251 252 477 476
+290 1 quad 252 253 478 477
+291 1 quad 253 254 479 478
+292 1 quad 254 255 480 479
+293 1 quad 255 256 481 480
+294 1 quad 256 257 482 481
+295 1 quad 257 258 483 482
+296 1 quad 258 259 484 483
+297 1 quad 259 260 485 484
+298 1 quad 260 261 486 485
+299 1 quad 261 262 487 486
+300 1 quad 262 263 488 487
+301 1 quad 263 264 489 488
+302 1 quad 264 265 490 489
+303 1 quad 265 266 491 490
+304 1 quad 266 267 492 491
+305 1 quad 267 268 493 492
+306 1 quad 268 269 494 493
+307 1 quad 269 270 495 494
+308 1 quad 270 271 496 495
+309 1 quad 271 272 497 496
+310 1 quad 272 273 498 497
+311 1 quad 273 274 499 498
+312 1 quad 274 275 500 499
+313 1 quad 275 276 501 500
+314 1 quad 276 277 553 501
+315 1 quad 277 278 502 553
+316 1 quad 278 279 503 502
+317 1 quad 302 350 504 301
+318 1 quad 350 351 505 504
+319 1 quad 351 352 506 505
+320 1 quad 352 353 507 506
+321 1 quad 353 354 508 507
+322 1 quad 354 355 509 508
+323 1 quad 355 356 510 509
+324 1 quad 356 357 511 510
+325 1 quad 357 358 512 511
+326 1 quad 358 359 513 512
+327 1 quad 359 360 514 513
+328 1 quad 360 361 515 514
+329 1 quad 361 362 516 515
+330 1 quad 362 363 517 516
+331 1 quad 363 364 518 517
+332 1 quad 364 365 519 518
+333 1 quad 365 366 520 519
+334 1 quad 366 367 521 520
+335 1 quad 367 368 522 521
+336 1 quad 368 369 523 522
+337 1 quad 369 370 524 523
+338 1 quad 370 371 525 524
+339 1 quad 371 372 526 525
+340 1 quad 372 373 527 526
+341 1 quad 373 374 528 527
+342 1 quad 374 375 529 528
+343 1 quad 375 376 530 529
+344 1 quad 376 377 531 530
+345 1 quad 377 378 532 531
+346 1 quad 378 379 444 532
+347 1 quad 280 281 533 443
+348 1 quad 281 282 534 533
+349 1 quad 282 283 535 534
+350 1 quad 283 284 536 535
+351 1 quad 284 285 537 536
+352 1 quad 285 286 538 537
+353 1 quad 286 287 539 538
+354 1 quad 287 288 540 539
+355 1 quad 288 289 541 540
+356 1 quad 289 290 542 541
+357 1 quad 290 291 543 542
+358 1 quad 291 292 544 543
+359 1 quad 292 293 545 544
+360 1 quad 293 294 546 545
+361 1 quad 294 295 547 546
+362 1 quad 295 296 548 547
+363 1 quad 296 297 549 548
+364 1 quad 297 298 550 549
+365 1 quad 298 299 551 550
+366 1 quad 299 300 552 551
+367 1 quad 300 301 504 552
+368 1 quad 420 421 590 419
+369 1 quad 421 422 423 590
+370 1 quad 503 380 553 502
+371 1 quad 380 381 501 553
+372 1 quad 381 382 670 501
+373 1 quad 382 383 554 670
+374 1 quad 383 384 555 554
+375 1 quad 384 385 556 555
+376 1 quad 385 386 557 556
+377 1 quad 386 387 558 557
+378 1 quad 387 388 559 558
+379 1 quad 388 389 560 559
+380 1 quad 389 390 561 560
+381 1 quad 390 391 562 561
+382 1 quad 391 392 563 562
+383 1 quad 392 393 564 563
+384 1 quad 393 394 565 564
+385 1 quad 394 395 566 565
+386 1 quad 395 396 567 566
+387 1 quad 396 397 568 567
+388 1 quad 397 398 569 568
+389 1 quad 398 399 570 569
+390 1 quad 399 400 571 570
+391 1 quad 400 401 572 571
+392 1 quad 401 402 573 572
+393 1 quad 402 403 574 573
+394 1 quad 403 404 575 574
+395 1 quad 404 405 576 575
+396 1 quad 405 406 577 576
+397 1 quad 406 407 578 577
+398 1 quad 407 408 579 578
+399 1 quad 408 409 580 579
+400 1 quad 409 410 581 580
+401 1 quad 410 411 582 581
+402 1 quad 411 412 583 582
+403 1 quad 412 413 584 583
+404 1 quad 413 414 585 584
+405 1 quad 414 415 586 585
+406 1 quad 415 416 587 586
+407 1 quad 416 417 588 587
+408 1 quad 417 418 589 588
+409 1 quad 418 419 590 589
+410 1 quad 423 424 591 590
+411 1 quad 424 425 592 591
+412 1 quad 425 426 593 592
+413 1 quad 426 427 594 593
+414 1 quad 427 428 595 594
+415 1 quad 428 429 596 595
+416 1 quad 429 598 597 596
+417 1 quad 429 430 599 598
+418 1 quad 430 431 600 599
+419 1 quad 431 432 601 600
+420 1 quad 432 433 602 601
+421 1 quad 433 434 603 602
+422 1 quad 434 435 604 603
+423 1 quad 435 436 605 604
+424 1 quad 436 607 606 605
+425 1 quad 436 437 608 607
+426 1 quad 437 438 609 608
+427 1 quad 438 439 610 609
+428 1 quad 439 440 611 610
+429 1 quad 440 441 612 611
+430 1 quad 441 442 613 612
+431 1 quad 442 443 533 613
+432 1 quad 444 445 614 532
+433 1 quad 445 446 615 614
+434 1 quad 446 447 616 615
+435 1 quad 447 448 617 616
+436 1 quad 448 449 618 617
+437 1 quad 449 450 619 618
+438 1 quad 450 451 620 619
+439 1 quad 451 452 621 620
+440 1 quad 452 453 622 621
+441 1 quad 453 454 623 622
+442 1 quad 454 455 624 623
+443 1 quad 455 456 625 624
+444 1 quad 456 457 626 625
+445 1 quad 457 458 627 626
+446 1 quad 458 459 628 627
+447 1 quad 459 460 629 628
+448 1 quad 460 461 630 629
+449 1 quad 461 462 631 630
+450 1 quad 462 463 632 631
+451 1 quad 463 464 633 632
+452 1 quad 464 465 634 633
+453 1 quad 465 466 635 634
+454 1 quad 466 467 636 635
+455 1 quad 467 468 637 636
+456 1 quad 468 469 638 637
+457 1 quad 469 470 639 638
+458 1 quad 470 471 640 639
+459 1 quad 471 472 641 640
+460 1 quad 472 473 642 641
+461 1 quad 473 474 643 642
+462 1 quad 474 475 644 643
+463 1 quad 475 476 645 644
+464 1 quad 476 477 646 645
+465 1 quad 477 478 647 646
+466 1 quad 478 479 648 647
+467 1 quad 479 480 649 648
+468 1 quad 480 481 650 649
+469 1 quad 481 482 651 650
+470 1 quad 482 483 652 651
+471 1 quad 483 484 653 652
+472 1 quad 484 485 654 653
+473 1 quad 485 486 655 654
+474 1 quad 486 487 656 655
+475 1 quad 487 488 657 656
+476 1 quad 488 489 658 657
+477 1 quad 489 490 659 658
+478 1 quad 490 491 660 659
+479 1 quad 491 492 661 660
+480 1 quad 492 493 662 661
+481 1 quad 493 494 663 662
+482 1 quad 494 495 664 663
+483 1 quad 495 496 665 664
+484 1 quad 496 497 666 665
+485 1 quad 497 498 667 666
+486 1 quad 498 499 668 667
+487 1 quad 499 500 669 668
+488 1 quad 500 501 670 669
+489 1 quad 533 534 671 613
+490 1 quad 534 535 672 671
+491 1 quad 535 536 673 672
+492 1 quad 536 537 674 673
+493 1 quad 537 538 675 674
+494 1 quad 538 539 676 675
+495 1 quad 539 540 677 676
+496 1 quad 540 541 678 677
+497 1 quad 541 542 679 678
+498 1 quad 542 543 680 679
+499 1 quad 543 544 681 680
+500 1 quad 544 545 682 681
+501 1 quad 545 546 683 682
+502 1 quad 546 547 684 683
+503 1 quad 547 548 685 684
+504 1 quad 548 549 686 685
+505 1 quad 549 550 687 686
+506 1 quad 550 551 688 687
+507 1 quad 551 552 689 688
+508 1 quad 552 504 505 689
+509 1 quad 505 506 690 689
+510 1 quad 506 507 691 690
+511 1 quad 507 508 692 691
+512 1 quad 508 509 693 692
+513 1 quad 509 510 694 693
+514 1 quad 510 511 695 694
+515 1 quad 511 512 696 695
+516 1 quad 512 513 697 696
+517 1 quad 513 514 698 697
+518 1 quad 514 515 699 698
+519 1 quad 515 516 700 699
+520 1 quad 516 517 701 700
+521 1 quad 517 518 702 701
+522 1 quad 518 519 703 702
+523 1 quad 519 520 704 703
+524 1 quad 520 521 705 704
+525 1 quad 521 522 706 705
+526 1 quad 522 523 707 706
+527 1 quad 523 524 708 707
+528 1 quad 524 525 709 708
+529 1 quad 525 526 710 709
+530 1 quad 526 527 711 710
+531 1 quad 527 528 712 711
+532 1 quad 528 529 713 712
+533 1 quad 529 530 714 713
+534 1 quad 530 531 715 714
+535 1 quad 531 532 614 715
+536 1 quad 590 591 716 589
+537 1 quad 591 592 717 716
+538 1 quad 592 593 718 717
+539 1 quad 593 594 719 718
+540 1 quad 594 595 720 719
+541 1 quad 595 596 721 720
+542 1 quad 596 597 722 721
+543 1 quad 597 724 723 722
+544 1 quad 597 598 725 724
+545 1 quad 598 599 726 725
+546 1 quad 599 600 727 726
+547 1 quad 600 601 728 727
+548 1 quad 601 602 729 728
+549 1 quad 602 603 730 729
+550 1 quad 603 604 731 730
+551 1 quad 604 605 732 731
+552 1 quad 605 606 733 732
+553 1 quad 606 735 734 733
+554 1 quad 606 607 736 735
+555 1 quad 607 608 737 736
+556 1 quad 608 609 738 737
+557 1 quad 609 610 739 738
+558 1 quad 610 611 740 739
+559 1 quad 611 612 741 740
+560 1 quad 612 613 671 741
+561 1 quad 614 615 742 715
+562 1 quad 615 616 743 742
+563 1 quad 616 617 744 743
+564 1 quad 617 618 745 744
+565 1 quad 618 619 746 745
+566 1 quad 619 620 747 746
+567 1 quad 620 621 748 747
+568 1 quad 621 622 749 748
+569 1 quad 622 623 750 749
+570 1 quad 623 624 751 750
+571 1 quad 624 625 752 751
+572 1 quad 625 626 753 752
+573 1 quad 626 627 754 753
+574 1 quad 627 628 755 754
+575 1 quad 628 629 756 755
+576 1 quad 629 630 757 756
+577 1 quad 630 631 758 757
+578 1 quad 631 632 759 758
+579 1 quad 632 633 760 759
+580 1 quad 633 634 761 760
+581 1 quad 634 635 762 761
+582 1 quad 635 636 763 762
+583 1 quad 636 637 764 763
+584 1 quad 637 638 765 764
+585 1 quad 638 639 766 765
+586 1 quad 639 640 767 766
+587 1 quad 640 641 768 767
+588 1 quad 641 642 769 768
+589 1 quad 642 643 770 769
+590 1 quad 643 644 771 770
+591 1 quad 644 645 772 771
+592 1 quad 645 646 773 772
+593 1 quad 646 647 774 773
+594 1 quad 647 648 775 774
+595 1 quad 648 649 776 775
+596 1 quad 649 650 777 776
+597 1 quad 650 651 778 777
+598 1 quad 651 652 779 778
+599 1 quad 652 653 780 779
+600 1 quad 653 654 781 780
+601 1 quad 654 655 782 781
+602 1 quad 655 656 783 782
+603 1 quad 656 657 784 783
+604 1 quad 657 658 785 784
+605 1 quad 658 659 786 785
+606 1 quad 659 660 787 786
+607 1 quad 660 661 788 787
+608 1 quad 661 662 789 788
+609 1 quad 662 663 790 789
+610 1 quad 663 664 791 790
+611 1 quad 664 665 792 791
+612 1 quad 665 666 793 792
+613 1 quad 666 667 794 793
+614 1 quad 667 668 795 794
+615 1 quad 668 669 796 795
+616 1 quad 669 670 554 796
+617 1 quad 671 672 797 741
+618 1 quad 672 673 798 797
+619 1 quad 673 674 799 798
+620 1 quad 674 675 800 799
+621 1 quad 675 676 801 800
+622 1 quad 676 677 802 801
+623 1 quad 677 678 803 802
+624 1 quad 678 679 804 803
+625 1 quad 679 680 805 804
+626 1 quad 680 681 806 805
+627 1 quad 681 682 807 806
+628 1 quad 682 683 808 807
+629 1 quad 683 684 809 808
+630 1 quad 684 685 810 809
+631 1 quad 685 686 811 810
+632 1 quad 686 687 812 811
+633 1 quad 687 688 813 812
+634 1 quad 688 689 690 813
+635 1 quad 554 555 814 796
+636 1 quad 555 556 815 814
+637 1 quad 556 557 816 815
+638 1 quad 557 558 817 816
+639 1 quad 558 559 818 817
+640 1 quad 559 560 819 818
+641 1 quad 560 561 820 819
+642 1 quad 561 562 821 820
+643 1 quad 562 563 822 821
+644 1 quad 563 564 823 822
+645 1 quad 564 565 824 823
+646 1 quad 565 566 825 824
+647 1 quad 566 567 826 825
+648 1 quad 567 568 827 826
+649 1 quad 568 569 828 827
+650 1 quad 569 570 829 828
+651 1 quad 570 571 830 829
+652 1 quad 571 572 831 830
+653 1 quad 572 573 832 831
+654 1 quad 573 574 833 832
+655 1 quad 574 575 834 833
+656 1 quad 575 576 835 834
+657 1 quad 576 577 836 835
+658 1 quad 577 578 837 836
+659 1 quad 578 579 838 837
+660 1 quad 579 580 839 838
+661 1 quad 580 581 840 839
+662 1 quad 581 582 841 840
+663 1 quad 582 583 842 841
+664 1 quad 583 584 843 842
+665 1 quad 584 585 844 843
+666 1 quad 585 586 845 844
+667 1 quad 586 587 1723 845
+668 1 quad 587 588 1725 1723
+669 1 quad 588 589 716 1725
+670 1 quad 690 691 848 813
+671 1 quad 691 692 849 848
+672 1 quad 692 693 850 849
+673 1 quad 693 694 851 850
+674 1 quad 694 695 852 851
+675 1 quad 695 696 853 852
+676 1 quad 696 697 854 853
+677 1 quad 697 698 855 854
+678 1 quad 698 699 856 855
+679 1 quad 699 700 857 856
+680 1 quad 700 701 858 857
+681 1 quad 701 702 859 858
+682 1 quad 702 703 860 859
+683 1 quad 703 704 861 860
+684 1 quad 704 705 862 861
+685 1 quad 705 706 863 862
+686 1 quad 706 707 864 863
+687 1 quad 707 708 865 864
+688 1 quad 708 709 866 865
+689 1 quad 709 710 867 866
+690 1 quad 710 711 868 867
+691 1 quad 711 712 869 868
+692 1 quad 712 713 870 869
+693 1 quad 713 714 871 870
+694 1 quad 714 715 742 871
+695 1 quad 718 719 846 847
+696 1 quad 719 720 872 846
+697 1 quad 720 721 873 872
+698 1 quad 721 722 874 873
+699 1 quad 722 723 875 874
+700 1 quad 723 877 876 875
+701 1 quad 723 724 878 877
+702 1 quad 724 725 879 878
+703 1 quad 725 726 880 879
+704 1 quad 726 727 881 880
+705 1 quad 727 728 882 881
+706 1 quad 728 729 883 882
+707 1 quad 729 730 884 883
+708 1 quad 730 731 885 884
+709 1 quad 731 732 886 885
+710 1 quad 732 733 887 886
+711 1 quad 733 734 888 887
+712 1 quad 734 890 889 888
+713 1 quad 734 735 891 890
+714 1 quad 735 736 892 891
+715 1 quad 736 737 893 892
+716 1 quad 737 738 894 893
+717 1 quad 738 739 895 894
+718 1 quad 739 740 896 895
+719 1 quad 740 741 797 896
+720 1 quad 796 814 897 795
+721 1 quad 814 815 898 897
+722 1 quad 815 816 1040 898
+723 1 quad 816 817 976 1040
+724 1 quad 817 818 899 976
+725 1 quad 818 819 900 899
+726 1 quad 819 820 901 900
+727 1 quad 820 821 902 901
+728 1 quad 821 822 903 902
+729 1 quad 822 823 904 903
+730 1 quad 823 824 905 904
+731 1 quad 824 825 906 905
+732 1 quad 825 826 907 906
+733 1 quad 826 827 908 907
+734 1 quad 827 828 909 908
+735 1 quad 828 829 910 909
+736 1 quad 829 830 911 910
+737 1 quad 830 831 912 911
+738 1 quad 831 832 913 912
+739 1 quad 832 833 914 913
+740 1 quad 833 834 915 914
+741 1 quad 834 835 916 915
+742 1 quad 835 836 917 916
+743 1 quad 836 837 918 917
+744 1 quad 837 838 919 918
+745 1 quad 838 839 920 919
+746 1 quad 839 840 921 920
+747 1 quad 840 841 922 921
+748 1 quad 841 842 923 922
+749 1 quad 842 843 924 923
+750 1 quad 843 844 925 924
+751 1 quad 844 845 926 925
+752 1 quad 742 743 927 871
+753 1 quad 743 744 928 927
+754 1 quad 744 745 929 928
+755 1 quad 745 746 930 929
+756 1 quad 746 747 931 930
+757 1 quad 747 748 932 931
+758 1 quad 748 749 933 932
+759 1 quad 749 750 934 933
+760 1 quad 750 751 935 934
+761 1 quad 751 752 936 935
+762 1 quad 752 753 937 936
+763 1 quad 753 754 938 937
+764 1 quad 754 755 939 938
+765 1 quad 755 756 940 939
+766 1 quad 756 757 941 940
+767 1 quad 757 758 942 941
+768 1 quad 758 759 943 942
+769 1 quad 759 760 944 943
+770 1 quad 760 761 945 944
+771 1 quad 761 762 946 945
+772 1 quad 762 763 947 946
+773 1 quad 763 764 948 947
+774 1 quad 764 765 949 948
+775 1 quad 765 766 950 949
+776 1 quad 766 767 951 950
+777 1 quad 767 768 952 951
+778 1 quad 768 769 953 952
+779 1 quad 769 770 954 953
+780 1 quad 770 771 955 954
+781 1 quad 771 772 956 955
+782 1 quad 772 773 957 956
+783 1 quad 773 774 958 957
+784 1 quad 774 775 959 958
+785 1 quad 775 776 960 959
+786 1 quad 776 777 961 960
+787 1 quad 777 778 962 961
+788 1 quad 778 779 963 962
+789 1 quad 779 780 964 963
+790 1 quad 780 781 965 964
+791 1 quad 781 782 966 965
+792 1 quad 782 783 967 966
+793 1 quad 783 784 968 967
+794 1 quad 784 785 969 968
+795 1 quad 785 786 970 969
+796 1 quad 786 787 971 970
+797 1 quad 787 788 972 971
+798 1 quad 788 789 973 972
+799 1 quad 789 790 974 973
+800 1 quad 790 791 975 974
+801 1 quad 791 792 976 975
+802 1 quad 797 798 977 896
+803 1 quad 798 799 978 977
+804 1 quad 799 800 979 978
+805 1 quad 800 801 980 979
+806 1 quad 801 802 981 980
+807 1 quad 802 803 982 981
+808 1 quad 803 804 983 982
+809 1 quad 804 805 984 983
+810 1 quad 805 806 985 984
+811 1 quad 806 807 986 985
+812 1 quad 807 808 987 986
+813 1 quad 808 809 988 987
+814 1 quad 809 810 989 988
+815 1 quad 810 811 990 989
+816 1 quad 811 812 991 990
+817 1 quad 812 813 848 991
+818 1 quad 848 849 992 991
+819 1 quad 849 850 993 992
+820 1 quad 850 851 994 993
+821 1 quad 851 852 995 994
+822 1 quad 852 853 996 995
+823 1 quad 853 854 997 996
+824 1 quad 854 855 998 997
+825 1 quad 855 856 999 998
+826 1 quad 856 857 1000 999
+827 1 quad 857 858 1001 1000
+828 1 quad 858 859 1002 1001
+829 1 quad 859 860 1003 1002
+830 1 quad 860 861 1004 1003
+831 1 quad 861 862 1005 1004
+832 1 quad 862 863 1006 1005
+833 1 quad 863 864 1007 1006
+834 1 quad 864 865 1008 1007
+835 1 quad 865 866 1009 1008
+836 1 quad 866 867 1010 1009
+837 1 quad 867 868 1011 1010
+838 1 quad 868 869 1012 1011
+839 1 quad 869 870 1013 1012
+840 1 quad 870 871 927 1013
+841 1 quad 872 873 925 926
+842 1 quad 873 874 1014 925
+843 1 quad 874 875 1015 1014
+844 1 quad 875 876 1016 1015
+845 1 quad 876 1018 1017 1016
+846 1 quad 876 877 1019 1018
+847 1 quad 877 878 1020 1019
+848 1 quad 878 879 1021 1020
+849 1 quad 879 880 1022 1021
+850 1 quad 880 881 1023 1022
+851 1 quad 881 882 1024 1023
+852 1 quad 882 883 1025 1024
+853 1 quad 883 884 1026 1025
+854 1 quad 884 885 1027 1026
+855 1 quad 885 886 1028 1027
+856 1 quad 886 887 1029 1028
+857 1 quad 887 888 1030 1029
+858 1 quad 888 889 1031 1030
+859 1 quad 889 1033 1032 1031
+860 1 quad 889 890 1034 1033
+861 1 quad 890 891 1035 1034
+862 1 quad 891 892 1036 1035
+863 1 quad 892 893 1037 1036
+864 1 quad 893 894 1038 1037
+865 1 quad 894 895 1039 1038
+866 1 quad 895 896 977 1039
+867 1 quad 795 897 898 794
+868 1 quad 792 793 1040 976
+869 1 quad 1040 793 794 898
+870 1 quad 927 928 1041 1013
+871 1 quad 928 929 1042 1041
+872 1 quad 929 930 1043 1042
+873 1 quad 930 931 1044 1043
+874 1 quad 931 932 1045 1044
+875 1 quad 932 933 1046 1045
+876 1 quad 933 934 1047 1046
+877 1 quad 934 935 1048 1047
+878 1 quad 935 936 1049 1048
+879 1 quad 936 937 1050 1049
+880 1 quad 937 938 1051 1050
+881 1 quad 938 939 1052 1051
+882 1 quad 939 940 1053 1052
+883 1 quad 940 941 1054 1053
+884 1 quad 941 942 1055 1054
+885 1 quad 942 943 1056 1055
+886 1 quad 943 944 1057 1056
+887 1 quad 944 945 1058 1057
+888 1 quad 945 946 1059 1058
+889 1 quad 946 947 1060 1059
+890 1 quad 947 948 1061 1060
+891 1 quad 948 949 1062 1061
+892 1 quad 949 950 1063 1062
+893 1 quad 950 951 1064 1063
+894 1 quad 951 952 1065 1064
+895 1 quad 952 953 1066 1065
+896 1 quad 953 954 1067 1066
+897 1 quad 954 955 1068 1067
+898 1 quad 955 956 1069 1068
+899 1 quad 956 957 1070 1069
+900 1 quad 957 958 1071 1070
+901 1 quad 958 959 1072 1071
+902 1 quad 959 960 1073 1072
+903 1 quad 960 961 1074 1073
+904 1 quad 961 962 1075 1074
+905 1 quad 962 963 1076 1075
+906 1 quad 963 964 1077 1076
+907 1 quad 964 965 1078 1077
+908 1 quad 965 966 1079 1078
+909 1 quad 966 967 1080 1079
+910 1 quad 967 968 1081 1080
+911 1 quad 968 969 1082 1081
+912 1 quad 969 970 1083 1082
+913 1 quad 970 971 1084 1083
+914 1 quad 971 972 1085 1084
+915 1 quad 972 973 1086 1085
+916 1 quad 973 974 1087 1086
+917 1 quad 974 975 1088 1087
+918 1 quad 975 976 899 1088
+919 1 quad 977 978 1089 1039
+920 1 quad 978 979 1090 1089
+921 1 quad 979 980 1091 1090
+922 1 quad 980 981 1092 1091
+923 1 quad 981 982 1093 1092
+924 1 quad 982 983 1094 1093
+925 1 quad 983 984 1095 1094
+926 1 quad 984 985 1096 1095
+927 1 quad 985 986 1097 1096
+928 1 quad 986 987 1098 1097
+929 1 quad 987 988 1099 1098
+930 1 quad 988 989 1100 1099
+931 1 quad 989 990 1101 1100
+932 1 quad 990 991 992 1101
+933 1 quad 899 900 1087 1088
+934 1 quad 900 901 1086 1087
+935 1 quad 901 902 1213 1086
+936 1 quad 902 903 1102 1213
+937 1 quad 903 904 1103 1102
+938 1 quad 904 905 1104 1103
+939 1 quad 905 906 1105 1104
+940 1 quad 906 907 1106 1105
+941 1 quad 907 908 1107 1106
+942 1 quad 908 909 1108 1107
+943 1 quad 909 910 1109 1108
+944 1 quad 910 911 1110 1109
+945 1 quad 911 912 1111 1110
+946 1 quad 912 913 1112 1111
+947 1 quad 913 914 1113 1112
+948 1 quad 914 915 1114 1113
+949 1 quad 915 916 1115 1114
+950 1 quad 916 917 1116 1115
+951 1 quad 917 918 1117 1116
+952 1 quad 918 919 1118 1117
+953 1 quad 919 920 1119 1118
+954 1 quad 920 921 1120 1119
+955 1 quad 921 922 1121 1120
+956 1 quad 922 923 1143 1121
+957 1 quad 923 924 1122 1143
+958 1 quad 924 925 1014 1122
+959 1 quad 992 993 1123 1101
+960 1 quad 993 994 1124 1123
+961 1 quad 994 995 1125 1124
+962 1 quad 995 996 1126 1125
+963 1 quad 996 997 1127 1126
+964 1 quad 997 998 1128 1127
+965 1 quad 998 999 1129 1128
+966 1 quad 999 1000 1130 1129
+967 1 quad 1000 1001 1131 1130
+968 1 quad 1001 1002 1132 1131
+969 1 quad 1002 1003 1133 1132
+970 1 quad 1003 1004 1134 1133
+971 1 quad 1004 1005 1135 1134
+972 1 quad 1005 1006 1136 1135
+973 1 quad 1006 1007 1137 1136
+974 1 quad 1007 1008 1138 1137
+975 1 quad 1008 1009 1139 1138
+976 1 quad 1009 1010 1140 1139
+977 1 quad 1010 1011 1141 1140
+978 1 quad 1011 1012 1142 1141
+979 1 quad 1012 1013 1041 1142
+980 1 quad 1014 1015 1143 1122
+981 1 quad 1015 1016 1121 1143
+982 1 quad 1016 1017 1244 1121
+983 1 quad 1017 1145 1144 1244
+984 1 quad 1017 1018 1146 1145
+985 1 quad 1018 1019 1147 1146
+986 1 quad 1019 1020 1148 1147
+987 1 quad 1020 1021 1149 1148
+988 1 quad 1021 1022 1150 1149
+989 1 quad 1022 1023 1151 1150
+990 1 quad 1023 1024 1152 1151
+991 1 quad 1024 1025 1153 1152
+992 1 quad 1025 1026 1154 1153
+993 1 quad 1026 1027 1155 1154
+994 1 quad 1027 1028 1156 1155
+995 1 quad 1028 1029 1157 1156
+996 1 quad 1029 1030 1158 1157
+997 1 quad 1030 1031 1159 1158
+998 1 quad 1031 1032 1160 1159
+999 1 quad 1032 1162 1161 1160
+1000 1 quad 1032 1033 1163 1162
+1001 1 quad 1033 1034 1164 1163
+1002 1 quad 1034 1035 1165 1164
+1003 1 quad 1035 1036 1166 1165
+1004 1 quad 1036 1037 1167 1166
+1005 1 quad 1037 1038 1168 1167
+1006 1 quad 1038 1039 1089 1168
+1007 1 quad 1041 1042 1169 1142
+1008 1 quad 1042 1043 1170 1169
+1009 1 quad 1043 1044 1171 1170
+1010 1 quad 1044 1045 1172 1171
+1011 1 quad 1045 1046 1173 1172
+1012 1 quad 1046 1047 1174 1173
+1013 1 quad 1047 1048 1175 1174
+1014 1 quad 1048 1049 1176 1175
+1015 1 quad 1049 1050 1177 1176
+1016 1 quad 1050 1051 1178 1177
+1017 1 quad 1051 1052 1179 1178
+1018 1 quad 1052 1053 1180 1179
+1019 1 quad 1053 1054 1181 1180
+1020 1 quad 1054 1055 1182 1181
+1021 1 quad 1055 1056 1183 1182
+1022 1 quad 1056 1057 1184 1183
+1023 1 quad 1057 1058 1185 1184
+1024 1 quad 1058 1059 1186 1185
+1025 1 quad 1059 1060 1187 1186
+1026 1 quad 1060 1061 1188 1187
+1027 1 quad 1061 1062 1189 1188
+1028 1 quad 1062 1063 1190 1189
+1029 1 quad 1063 1064 1191 1190
+1030 1 quad 1064 1065 1192 1191
+1031 1 quad 1065 1066 1193 1192
+1032 1 quad 1066 1067 1194 1193
+1033 1 quad 1067 1068 1195 1194
+1034 1 quad 1068 1069 1196 1195
+1035 1 quad 1069 1070 1197 1196
+1036 1 quad 1070 1071 1198 1197
+1037 1 quad 1071 1072 1199 1198
+1038 1 quad 1072 1073 1200 1199
+1039 1 quad 1073 1074 1201 1200
+1040 1 quad 1074 1075 1202 1201
+1041 1 quad 1075 1076 1203 1202
+1042 1 quad 1076 1077 1204 1203
+1043 1 quad 1077 1078 1205 1204
+1044 1 quad 1078 1079 1206 1205
+1045 1 quad 1079 1080 1207 1206
+1046 1 quad 1080 1081 1208 1207
+1047 1 quad 1081 1082 1209 1208
+1048 1 quad 1082 1083 1210 1209
+1049 1 quad 1083 1084 1211 1210
+1050 1 quad 1084 1085 1212 1211
+1051 1 quad 1085 1086 1213 1212
+1052 1 quad 1089 1090 1214 1168
+1053 1 quad 1090 1091 1215 1214
+1054 1 quad 1091 1092 1216 1215
+1055 1 quad 1092 1093 1217 1216
+1056 1 quad 1093 1094 1218 1217
+1057 1 quad 1094 1095 1219 1218
+1058 1 quad 1095 1096 1220 1219
+1059 1 quad 1096 1097 1221 1220
+1060 1 quad 1097 1098 1222 1221
+1061 1 quad 1098 1099 1223 1222
+1062 1 quad 1099 1100 1224 1223
+1063 1 quad 1100 1101 1123 1224
+1064 1 quad 1213 1102 1225 1212
+1065 1 quad 1102 1103 1226 1225
+1066 1 quad 1103 1104 1227 1226
+1067 1 quad 1104 1105 1228 1227
+1068 1 quad 1105 1106 1229 1228
+1069 1 quad 1106 1107 1230 1229
+1070 1 quad 1107 1108 1231 1230
+1071 1 quad 1108 1109 1232 1231
+1072 1 quad 1109 1110 1233 1232
+1073 1 quad 1110 1111 1234 1233
+1074 1 quad 1111 1112 1235 1234
+1075 1 quad 1112 1113 1236 1235
+1076 1 quad 1113 1114 1237 1236
+1077 1 quad 1114 1115 1238 1237
+1078 1 quad 1115 1116 1239 1238
+1079 1 quad 1116 1117 1240 1239
+1080 1 quad 1117 1118 1241 1240
+1081 1 quad 1118 1119 1242 1241
+1082 1 quad 1119 1120 1243 1242
+1083 1 quad 1120 1121 1244 1243
+1084 1 quad 1123 1124 1245 1224
+1085 1 quad 1124 1125 1246 1245
+1086 1 quad 1125 1126 1247 1246
+1087 1 quad 1126 1127 1248 1247
+1088 1 quad 1127 1128 1249 1248
+1089 1 quad 1128 1129 1250 1249
+1090 1 quad 1129 1130 1251 1250
+1091 1 quad 1130 1131 1252 1251
+1092 1 quad 1131 1132 1253 1252
+1093 1 quad 1132 1133 1254 1253
+1094 1 quad 1133 1134 1255 1254
+1095 1 quad 1134 1135 1256 1255
+1096 1 quad 1135 1136 1257 1256
+1097 1 quad 1136 1137 1258 1257
+1098 1 quad 1137 1138 1259 1258
+1099 1 quad 1138 1139 1260 1259
+1100 1 quad 1139 1140 1261 1260
+1101 1 quad 1140 1141 1262 1261
+1102 1 quad 1141 1142 1169 1262
+1103 1 quad 1244 1144 1263 1243
+1104 1 quad 1243 1263 1241 1242
+1105 1 quad 1263 1264 1354 1241
+1106 1 quad 1263 1144 1265 1264
+1107 1 quad 1144 1145 1266 1265
+1108 1 quad 1146 1147 1266 1145
+1109 1 quad 1147 1148 1267 1266
+1110 1 quad 1148 1149 1268 1267
+1111 1 quad 1149 1150 1269 1268
+1112 1 quad 1150 1151 1270 1269
+1113 1 quad 1151 1152 1271 1270
+1114 1 quad 1152 1153 1272 1271
+1115 1 quad 1154 1155 1274 1273
+1116 1 quad 1155 1156 1730 1274
+1117 1 quad 1156 1157 1275 1730
+1118 1 quad 1157 1158 1276 1275
+1119 1 quad 1158 1159 1277 1276
+1120 1 quad 1159 1160 1278 1277
+1121 1 quad 1160 1161 1279 1278
+1122 1 quad 1161 1281 1280 1279
+1123 1 quad 1161 1162 1282 1281
+1124 1 quad 1162 1163 1283 1282
+1125 1 quad 1163 1164 1284 1283
+1126 1 quad 1164 1165 1285 1284
+1127 1 quad 1165 1166 1286 1285
+1128 1 quad 1166 1167 1287 1286
+1129 1 quad 1167 1168 1214 1287
+1130 1 quad 1153 1154 1273 1288
+1131 1 quad 1272 1153 1288 1289
+1132 1 quad 1169 1170 1290 1262
+1133 1 quad 1170 1171 1291 1290
+1134 1 quad 1171 1172 1292 1291
+1135 1 quad 1172 1173 1293 1292
+1136 1 quad 1173 1174 1294 1293
+1137 1 quad 1174 1175 1295 1294
+1138 1 quad 1175 1176 1296 1295
+1139 1 quad 1176 1177 1297 1296
+1140 1 quad 1177 1178 1298 1297
+1141 1 quad 1178 1179 1299 1298
+1142 1 quad 1179 1180 1300 1299
+1143 1 quad 1180 1181 1301 1300
+1144 1 quad 1181 1182 1302 1301
+1145 1 quad 1182 1183 1303 1302
+1146 1 quad 1183 1184 1304 1303
+1147 1 quad 1184 1185 1305 1304
+1148 1 quad 1185 1186 1306 1305
+1149 1 quad 1186 1187 1307 1306
+1150 1 quad 1187 1188 1308 1307
+1151 1 quad 1188 1189 1309 1308
+1152 1 quad 1189 1190 1310 1309
+1153 1 quad 1190 1191 1311 1310
+1154 1 quad 1191 1192 1312 1311
+1155 1 quad 1192 1193 1313 1312
+1156 1 quad 1193 1194 1314 1313
+1157 1 quad 1194 1195 1315 1314
+1158 1 quad 1195 1196 1316 1315
+1159 1 quad 1196 1197 1317 1316
+1160 1 quad 1197 1198 1318 1317
+1161 1 quad 1198 1199 1319 1318
+1162 1 quad 1199 1200 1320 1319
+1163 1 quad 1200 1201 1321 1320
+1164 1 quad 1201 1202 1322 1321
+1165 1 quad 1202 1203 1323 1322
+1166 1 quad 1203 1204 1324 1323
+1167 1 quad 1204 1205 1325 1324
+1168 1 quad 1205 1206 1326 1325
+1169 1 quad 1206 1207 1327 1326
+1170 1 quad 1207 1208 1328 1327
+1171 1 quad 1208 1209 1329 1328
+1172 1 quad 1209 1210 1330 1329
+1173 1 quad 1214 1215 1331 1287
+1174 1 quad 1215 1216 1332 1331
+1175 1 quad 1216 1217 1333 1332
+1176 1 quad 1217 1218 1334 1333
+1177 1 quad 1218 1219 1335 1334
+1178 1 quad 1219 1220 1336 1335
+1179 1 quad 1220 1221 1337 1336
+1180 1 quad 1221 1222 1338 1337
+1181 1 quad 1222 1223 1339 1338
+1182 1 quad 1223 1224 1245 1339
+1183 1 quad 1212 1225 1340 1211
+1184 1 quad 1211 1340 1330 1210
+1185 1 quad 1225 1226 1330 1340
+1186 1 quad 1226 1227 1329 1330
+1187 1 quad 1227 1228 1341 1329
+1188 1 quad 1228 1229 1342 1341
+1189 1 quad 1229 1230 1343 1342
+1190 1 quad 1230 1231 1344 1343
+1191 1 quad 1231 1232 1345 1344
+1192 1 quad 1232 1233 1346 1345
+1193 1 quad 1233 1234 1347 1346
+1194 1 quad 1234 1235 1348 1347
+1195 1 quad 1235 1236 1349 1348
+1196 1 quad 1236 1237 1350 1349
+1197 1 quad 1237 1238 1351 1350
+1198 1 quad 1238 1239 1352 1351
+1199 1 quad 1239 1240 1353 1352
+1200 1 quad 1240 1241 1354 1353
+1201 1 quad 1245 1246 1355 1339
+1202 1 quad 1246 1247 1356 1355
+1203 1 quad 1247 1248 1357 1356
+1204 1 quad 1248 1249 1358 1357
+1205 1 quad 1249 1250 1359 1358
+1206 1 quad 1250 1251 1360 1359
+1207 1 quad 1251 1252 1361 1360
+1208 1 quad 1252 1253 1362 1361
+1209 1 quad 1253 1254 1363 1362
+1210 1 quad 1254 1255 1364 1363
+1211 1 quad 1255 1256 1365 1364
+1212 1 quad 1256 1257 1366 1365
+1213 1 quad 1257 1258 1367 1366
+1214 1 quad 1258 1259 1368 1367
+1215 1 quad 1259 1260 1369 1368
+1216 1 quad 1260 1261 1370 1369
+1217 1 quad 1261 1262 1290 1370
+1218 1 quad 1329 1341 1371 1328
+1219 1 quad 1341 1342 1372 1371
+1220 1 quad 1342 1343 1373 1372
+1221 1 quad 1343 1344 1374 1373
+1222 1 quad 1344 1345 1375 1374
+1223 1 quad 1345 1346 1376 1375
+1224 1 quad 1346 1347 1377 1376
+1225 1 quad 1347 1348 1378 1377
+1226 1 quad 1348 1349 1379 1378
+1227 1 quad 1349 1350 1380 1379
+1228 1 quad 1350 1351 1381 1380
+1229 1 quad 1351 1352 1382 1381
+1230 1 quad 1352 1353 1383 1382
+1231 1 quad 1353 1354 1384 1383
+1232 1 quad 1354 1264 1385 1384
+1233 1 quad 1264 1265 1386 1385
+1234 1 quad 1265 1266 1387 1386
+1235 1 quad 1266 1267 1388 1387
+1236 1 quad 1267 1268 1389 1388
+1237 1 quad 1268 1269 1390 1389
+1238 1 quad 1269 1270 1391 1390
+1239 1 quad 1270 1271 1392 1391
+1240 1 quad 1271 1272 1393 1392
+1241 1 quad 1272 1289 1394 1393
+1242 1 quad 1289 1396 1395 1394
+1243 1 quad 1289 1288 1397 1396
+1244 1 quad 1288 1273 1398 1397
+1245 1 quad 1328 1371 1399 1327
+1246 1 quad 1371 1372 1476 1399
+1247 1 quad 1372 1373 1475 1476
+1248 1 quad 1373 1374 1474 1475
+1249 1 quad 1374 1375 1473 1474
+1250 1 quad 1375 1376 1400 1473
+1251 1 quad 1376 1377 1401 1400
+1252 1 quad 1377 1378 1402 1401
+1253 1 quad 1378 1379 1403 1402
+1254 1 quad 1379 1380 1404 1403
+1255 1 quad 1380 1381 1405 1404
+1256 1 quad 1381 1382 1406 1405
+1257 1 quad 1382 1383 1407 1406
+1258 1 quad 1383 1384 1408 1407
+1259 1 quad 1384 1385 1409 1408
+1260 1 quad 1385 1386 1410 1409
+1261 1 quad 1386 1387 1411 1410
+1262 1 quad 1387 1388 1412 1411
+1263 1 quad 1388 1389 1413 1412
+1264 1 quad 1389 1390 1414 1413
+1265 1 quad 1390 1391 1415 1414
+1266 1 quad 1391 1392 1416 1415
+1267 1 quad 1392 1393 1417 1416
+1268 1 quad 1393 1394 1418 1417
+1269 1 quad 1394 1395 1419 1418
+1270 1 quad 1395 1421 1420 1419
+1271 1 quad 1395 1396 1422 1421
+1272 1 quad 1396 1397 1423 1422
+1273 1 quad 1397 1398 1424 1423
+1274 1 quad 1398 1731 1425 1424
+1275 1 quad 1398 1273 1274 1731
+1276 1 quad 1275 1276 1428 1427
+1277 1 quad 1276 1277 1429 1428
+1278 1 quad 1277 1278 1430 1429
+1279 1 quad 1278 1279 1431 1430
+1280 1 quad 1279 1280 1432 1431
+1281 1 quad 1280 1434 1433 1432
+1282 1 quad 1280 1281 1435 1434
+1283 1 quad 1281 1282 1436 1435
+1284 1 quad 1282 1283 1437 1436
+1285 1 quad 1283 1284 1438 1437
+1286 1 quad 1284 1285 1439 1438
+1287 1 quad 1285 1286 1440 1439
+1288 1 quad 1286 1287 1331 1440
+1289 1 quad 1290 1291 1441 1370
+1290 1 quad 1291 1292 1442 1441
+1291 1 quad 1292 1293 1443 1442
+1292 1 quad 1293 1294 1444 1443
+1293 1 quad 1294 1295 1445 1444
+1294 1 quad 1295 1296 1446 1445
+1295 1 quad 1296 1297 1447 1446
+1296 1 quad 1297 1298 1448 1447
+1297 1 quad 1298 1299 1449 1448
+1298 1 quad 1299 1300 1450 1449
+1299 1 quad 1300 1301 1451 1450
+1300 1 quad 1301 1302 1452 1451
+1301 1 quad 1302 1303 1453 1452
+1302 1 quad 1303 1304 1454 1453
+1303 1 quad 1304 1305 1455 1454
+1304 1 quad 1305 1306 1456 1455
+1305 1 quad 1306 1307 1457 1456
+1306 1 quad 1307 1308 1458 1457
+1307 1 quad 1308 1309 1459 1458
+1308 1 quad 1309 1310 1460 1459
+1309 1 quad 1310 1311 1461 1460
+1310 1 quad 1311 1312 1462 1461
+1311 1 quad 1312 1313 1463 1462
+1312 1 quad 1313 1314 1464 1463
+1313 1 quad 1314 1315 1465 1464
+1314 1 quad 1315 1316 1466 1465
+1315 1 quad 1316 1317 1467 1466
+1316 1 quad 1317 1318 1468 1467
+1317 1 quad 1318 1319 1469 1468
+1318 1 quad 1319 1320 1470 1469
+1319 1 quad 1320 1321 1471 1470
+1320 1 quad 1321 1322 1472 1471
+1321 1 quad 1322 1323 1473 1472
+1322 1 quad 1323 1324 1474 1473
+1323 1 quad 1324 1325 1475 1474
+1324 1 quad 1325 1326 1476 1475
+1325 1 quad 1476 1326 1327 1399
+1326 1 quad 1331 1332 1477 1440
+1327 1 quad 1332 1333 1478 1477
+1328 1 quad 1333 1334 1479 1478
+1329 1 quad 1334 1335 1480 1479
+1330 1 quad 1335 1336 1481 1480
+1331 1 quad 1336 1337 1482 1481
+1332 1 quad 1337 1338 1483 1482
+1333 1 quad 1338 1339 1355 1483
+1334 1 quad 1355 1356 1484 1483
+1335 1 quad 1356 1357 1485 1484
+1336 1 quad 1357 1358 1486 1485
+1337 1 quad 1358 1359 1487 1486
+1338 1 quad 1359 1360 1488 1487
+1339 1 quad 1360 1361 1489 1488
+1340 1 quad 1361 1362 1490 1489
+1341 1 quad 1362 1363 1491 1490
+1342 1 quad 1363 1364 1492 1491
+1343 1 quad 1364 1365 1493 1492
+1344 1 quad 1365 1366 1494 1493
+1345 1 quad 1366 1367 1495 1494
+1346 1 quad 1367 1368 1496 1495
+1347 1 quad 1368 1369 1497 1496
+1348 1 quad 1369 1370 1441 1497
+1349 1 quad 1400 1401 1624 1565
+1350 1 quad 1401 1402 1498 1624
+1351 1 quad 1402 1403 1499 1498
+1352 1 quad 1403 1404 1500 1499
+1353 1 quad 1404 1405 1501 1500
+1354 1 quad 1405 1406 1502 1501
+1355 1 quad 1406 1407 1503 1502
+1356 1 quad 1407 1408 1504 1503
+1357 1 quad 1410 1411 1505 1504
+1358 1 quad 1411 1412 1506 1505
+1359 1 quad 1412 1413 1507 1506
+1360 1 quad 1413 1414 1508 1507
+1361 1 quad 1414 1415 1509 1508
+1362 1 quad 1415 1416 1510 1509
+1363 1 quad 1416 1417 1511 1510
+1364 1 quad 1417 1418 1512 1511
+1365 1 quad 1418 1419 1513 1512
+1366 1 quad 1419 1420 1514 1513
+1367 1 quad 1420 1516 1515 1514
+1368 1 quad 1420 1421 1517 1516
+1369 1 quad 1421 1422 1518 1517
+1370 1 quad 1422 1423 1519 1518
+1371 1 quad 1423 1424 1520 1519
+1372 1 quad 1428 1429 1522 1521
+1373 1 quad 1429 1430 1523 1522
+1374 1 quad 1430 1431 1524 1523
+1375 1 quad 1431 1432 1525 1524
+1376 1 quad 1432 1433 1526 1525
+1377 1 quad 1433 1528 1527 1526
+1378 1 quad 1433 1434 1529 1528
+1379 1 quad 1434 1435 1530 1529
+1380 1 quad 1435 1436 1531 1530
+1381 1 quad 1436 1437 1532 1531
+1382 1 quad 1437 1438 1533 1532
+1383 1 quad 1438 1439 1534 1533
+1384 1 quad 1439 1440 1477 1534
+1385 1 quad 1441 1442 1535 1497
+1386 1 quad 1442 1443 1536 1535
+1387 1 quad 1443 1444 1537 1536
+1388 1 quad 1444 1445 1538 1537
+1389 1 quad 1445 1446 1539 1538
+1390 1 quad 1446 1447 1540 1539
+1391 1 quad 1447 1448 1541 1540
+1392 1 quad 1448 1449 1542 1541
+1393 1 quad 1449 1450 1543 1542
+1394 1 quad 1450 1451 1544 1543
+1395 1 quad 1451 1452 1545 1544
+1396 1 quad 1452 1453 1546 1545
+1397 1 quad 1453 1454 1547 1546
+1398 1 quad 1454 1455 1548 1547
+1399 1 quad 1455 1456 1549 1548
+1400 1 quad 1456 1457 1550 1549
+1401 1 quad 1457 1458 1551 1550
+1402 1 quad 1458 1459 1552 1551
+1403 1 quad 1459 1460 1553 1552
+1404 1 quad 1460 1461 1554 1553
+1405 1 quad 1461 1462 1555 1554
+1406 1 quad 1462 1463 1556 1555
+1407 1 quad 1463 1464 1557 1556
+1408 1 quad 1464 1465 1558 1557
+1409 1 quad 1465 1466 1559 1558
+1410 1 quad 1466 1467 1560 1559
+1411 1 quad 1467 1468 1561 1560
+1412 1 quad 1468 1469 1562 1561
+1413 1 quad 1469 1470 1563 1562
+1414 1 quad 1470 1471 1564 1563
+1415 1 quad 1471 1472 1565 1564
+1416 1 quad 1477 1478 1566 1534
+1417 1 quad 1478 1479 1567 1566
+1418 1 quad 1479 1480 1568 1567
+1419 1 quad 1480 1481 1569 1568
+1420 1 quad 1481 1482 1570 1569
+1421 1 quad 1482 1483 1484 1570
+1422 1 quad 1484 1485 1571 1570
+1423 1 quad 1485 1486 1572 1571
+1424 1 quad 1486 1487 1573 1572
+1425 1 quad 1487 1488 1574 1573
+1426 1 quad 1488 1489 1575 1574
+1427 1 quad 1489 1490 1576 1575
+1428 1 quad 1490 1491 1577 1576
+1429 1 quad 1491 1492 1578 1577
+1430 1 quad 1492 1493 1579 1578
+1431 1 quad 1493 1494 1580 1579
+1432 1 quad 1494 1495 1581 1580
+1433 1 quad 1495 1496 1582 1581
+1434 1 quad 1496 1497 1535 1582
+1435 1 quad 1521 1522 1583 1520
+1436 1 quad 1522 1523 1584 1583
+1437 1 quad 1523 1524 1585 1584
+1438 1 quad 1524 1525 1586 1585
+1439 1 quad 1525 1526 1587 1586
+1440 1 quad 1526 1527 1588 1587
+1441 1 quad 1527 1590 1589 1588
+1442 1 quad 1527 1528 1591 1590
+1443 1 quad 1528 1529 1592 1591
+1444 1 quad 1529 1530 1593 1592
+1445 1 quad 1530 1531 1594 1593
+1446 1 quad 1531 1532 1595 1594
+1447 1 quad 1532 1533 1596 1595
+1448 1 quad 1533 1534 1566 1596
+1449 1 quad 1535 1536 1597 1582
+1450 1 quad 1536 1537 1598 1597
+1451 1 quad 1537 1538 1599 1598
+1452 1 quad 1538 1539 1600 1599
+1453 1 quad 1539 1540 1601 1600
+1454 1 quad 1540 1541 1602 1601
+1455 1 quad 1541 1542 1603 1602
+1456 1 quad 1542 1543 1604 1603
+1457 1 quad 1543 1544 1605 1604
+1458 1 quad 1544 1545 1606 1605
+1459 1 quad 1545 1546 1607 1606
+1460 1 quad 1546 1547 1608 1607
+1461 1 quad 1547 1548 1609 1608
+1462 1 quad 1548 1549 1610 1609
+1463 1 quad 1549 1550 1611 1610
+1464 1 quad 1550 1551 1515 1611
+1465 1 quad 1551 1552 1651 1515
+1466 1 quad 1552 1553 1612 1651
+1467 1 quad 1553 1554 1613 1612
+1468 1 quad 1554 1555 1614 1613
+1469 1 quad 1555 1556 1615 1614
+1470 1 quad 1556 1557 1616 1615
+1471 1 quad 1557 1558 1617 1616
+1472 1 quad 1558 1559 1618 1617
+1473 1 quad 1559 1560 1619 1618
+1474 1 quad 1560 1561 1620 1619
+1475 1 quad 1561 1562 1621 1620
+1476 1 quad 1562 1563 1622 1621
+1477 1 quad 1563 1564 1623 1622
+1478 1 quad 1564 1565 1624 1623
+1479 1 quad 1566 1567 1625 1596
+1480 1 quad 1567 1568 1626 1625
+1481 1 quad 1568 1569 1627 1626
+1482 1 quad 1569 1570 1571 1627
+1483 1 quad 1571 1572 1628 1627
+1484 1 quad 1572 1573 1629 1628
+1485 1 quad 1573 1574 1630 1629
+1486 1 quad 1574 1575 1631 1630
+1487 1 quad 1575 1576 1632 1631
+1488 1 quad 1576 1577 1633 1632
+1489 1 quad 1577 1578 1634 1633
+1490 1 quad 1578 1579 1635 1634
+1491 1 quad 1579 1580 1636 1635
+1492 1 quad 1580 1581 1637 1636
+1493 1 quad 1581 1582 1597 1637
+1494 1 quad 1624 1498 1638 1623
+1495 1 quad 1498 1499 1674 1638
+1496 1 quad 1499 1500 1673 1674
+1497 1 quad 1501 1502 1640 1639
+1498 1 quad 1502 1503 1641 1640
+1499 1 quad 1505 1506 1642 1641
+1500 1 quad 1506 1507 1643 1642
+1501 1 quad 1507 1508 1644 1643
+1502 1 quad 1508 1509 1645 1644
+1503 1 quad 1509 1510 1646 1645
+1504 1 quad 1510 1511 1647 1646
+1505 1 quad 1511 1512 1648 1647
+1506 1 quad 1512 1513 1649 1648
+1507 1 quad 1513 1514 1650 1649
+1508 1 quad 1514 1515 1651 1650
+1509 1 quad 1520 1583 1652 1519
+1510 1 quad 1583 1584 1653 1652
+1511 1 quad 1584 1585 1654 1653
+1512 1 quad 1585 1586 1655 1654
+1513 1 quad 1586 1587 1656 1655
+1514 1 quad 1587 1588 1657 1656
+1515 1 quad 1588 1589 1658 1657
+1516 1 quad 1589 1660 1659 1658
+1517 1 quad 1589 1590 1661 1660
+1518 1 quad 1590 1591 1662 1661
+1519 1 quad 1591 1592 1663 1662
+1520 1 quad 1592 1593 1664 1663
+1521 1 quad 1593 1594 1665 1664
+1522 1 quad 1594 1595 1666 1665
+1523 1 quad 1595 1596 1625 1666
+1524 1 quad 1651 1612 1667 1650
+1525 1 quad 1612 1613 1696 1667
+1526 1 quad 1613 1614 1695 1696
+1527 1 quad 1614 1615 1694 1695
+1528 1 quad 1615 1616 1668 1694
+1529 1 quad 1616 1617 1669 1668
+1530 1 quad 1617 1618 1670 1669
+1531 1 quad 1618 1619 1671 1670
+1532 1 quad 1619 1620 1672 1671
+1533 1 quad 1620 1621 1673 1672
+1534 1 quad 1621 1622 1674 1673
+1535 1 quad 1674 1622 1623 1638
+1536 1 quad 1597 1598 1675 1637
+1537 1 quad 1598 1599 1676 1675
+1538 1 quad 1599 1600 1677 1676
+1539 1 quad 1600 1601 1678 1677
+1540 1 quad 1601 1602 1679 1678
+1541 1 quad 1602 1603 1680 1679
+1542 1 quad 1603 1604 1681 1680
+1543 1 quad 1604 1605 1682 1681
+1544 1 quad 1605 1606 1683 1682
+1545 1 quad 1606 1607 1684 1683
+1546 1 quad 1607 1608 1685 1684
+1547 1 quad 1608 1609 1686 1685
+1548 1 quad 1609 1610 1687 1686
+1549 1 quad 1610 1611 1688 1687
+1550 1 quad 1611 1515 1516 1688
+1551 1 quad 1625 1626 1689 1666
+1552 1 quad 1626 1627 1628 1689
+1553 1 quad 1628 1629 1719 1689
+1554 1 quad 1629 1630 1718 1719
+1555 1 quad 1630 1631 1717 1718
+1556 1 quad 1631 1632 1716 1717
+1557 1 quad 1632 1633 1715 1716
+1558 1 quad 1633 1634 1714 1715
+1559 1 quad 1634 1635 1713 1714
+1560 1 quad 1635 1636 1690 1713
+1561 1 quad 1636 1637 1675 1690
+1562 1 quad 1516 1517 1707 1688
+1563 1 quad 1517 1518 1706 1707
+1564 1 quad 1518 1519 1652 1706
+1565 1 quad 1641 1642 1711 1640
+1566 1 quad 1642 1643 1710 1711
+1567 1 quad 1643 1644 1691 1710
+1568 1 quad 1644 1645 1692 1691
+1569 1 quad 1645 1646 1693 1692
+1570 1 quad 1646 1647 1694 1693
+1571 1 quad 1647 1648 1695 1694
+1572 1 quad 1648 1649 1696 1695
+1573 1 quad 1696 1649 1650 1667
+1574 1 quad 1675 1676 1697 1690
+1575 1 quad 1676 1677 1698 1697
+1576 1 quad 1677 1678 1699 1698
+1577 1 quad 1678 1679 1700 1699
+1578 1 quad 1679 1680 1701 1700
+1579 1 quad 1680 1681 1702 1701
+1580 1 quad 1681 1682 1703 1702
+1581 1 quad 1682 1683 1704 1703
+1582 1 quad 1683 1684 1705 1704
+1583 1 quad 1684 1685 1706 1705
+1584 1 quad 1685 1686 1707 1706
+1585 1 quad 1688 1707 1686 1687
+1586 1 quad 1694 1668 1708 1693
+1587 1 quad 1668 1669 1709 1708
+1588 1 quad 1669 1670 1727 1709
+1589 1 quad 1670 1671 1729 1727
+1590 1 quad 1671 1672 1712 1729
+1591 1 quad 1652 1653 1705 1706
+1592 1 quad 1653 1654 1704 1705
+1593 1 quad 1654 1655 1703 1704
+1594 1 quad 1655 1656 1702 1703
+1595 1 quad 1656 1657 1701 1702
+1596 1 quad 1657 1658 1700 1701
+1597 1 quad 1658 1659 1699 1700
+1598 1 quad 1699 1659 1713 1698
+1599 1 quad 1659 1660 1714 1713
+1600 1 quad 1660 1661 1715 1714
+1601 1 quad 1661 1662 1716 1715
+1602 1 quad 1662 1663 1717 1716
+1603 1 quad 1663 1664 1718 1717
+1604 1 quad 1664 1665 1719 1718
+1605 1 quad 1719 1665 1666 1689
+1606 1 quad 1713 1690 1697 1698
+1607 1 quad 1691 1692 1721 1720
+1608 1 quad 1693 1708 1721 1692
+1609 1 quad 1708 1709 1720 1721
+1610 1 quad 420 419 344 345
+1611 1 quad 345 346 347 420
+1612 1 quad 872 926 1722 846
+1613 1 quad 926 845 1723 1722
+1614 1 quad 1722 1723 1725 1724
+1615 1 quad 846 1722 1724 847
+1616 1 quad 847 1724 717 718
+1617 1 quad 1724 1725 716 717
+1618 1 quad 1504 1505 1641 1503
+1619 1 quad 1504 1408 1409 1410
+1620 1 quad 1521 1520 1424 1425
+1621 1 quad 1425 1427 1428 1521
+1622 1 quad 1473 1400 1565 1472
+1623 1 quad 1691 1720 1726 1710
+1624 1 quad 1720 1709 1727 1726
+1625 1 quad 1726 1727 1729 1728
+1626 1 quad 1710 1726 1728 1711
+1627 1 quad 1711 1728 1639 1640
+1628 1 quad 1728 1729 1712 1639
+1629 1 quad 1274 1730 1426 1731
+1630 1 quad 1730 1275 1427 1426
+1631 1 quad 1426 1427 1425 1731
+1632 1 quad 1712 1672 1673 1500
+1633 1 quad 1712 1500 1501 1639
+1633 1 line 67 180
+1634 1 line 180 179
+1635 1 line 179 178
+1636 1 line 178 177
+1637 1 line 177 176
+1638 1 line 176 175
+1639 1 line 175 174
+1640 1 line 174 173
+1641 1 line 173 172
+1642 1 line 172 171
+1643 1 line 171 170
+1644 1 line 170 169
+1645 1 line 169 168
+1646 1 line 168 167
+1647 1 line 167 166
+1648 1 line 166 165
+1649 1 line 165 164
+1650 1 line 164 163
+1651 1 line 163 162
+1652 1 line 162 161
+1653 1 line 161 160
+1654 1 line 160 159
+1655 1 line 159 158
+1656 1 line 158 157
+1657 1 line 157 156
+1658 1 line 156 155
+1659 1 line 155 154
+1660 1 line 154 153
+1661 1 line 153 152
+1662 1 line 152 151
+1663 1 line 151 150
+1664 1 line 150 149
+1665 1 line 149 148
+1666 1 line 148 147
+1667 1 line 147 146
+1668 1 line 146 145
+1669 1 line 145 144
+1670 1 line 144 143
+1671 1 line 143 142
+1672 1 line 142 141
+1673 1 line 141 140
+1674 1 line 140 139
+1675 1 line 139 138
+1676 1 line 138 137
+1677 1 line 137 136
+1678 1 line 136 135
+1679 1 line 135 134
+1680 2 line 66 34
+1681 2 line 9 35
+1682 2 line 35 36
+1683 2 line 36 37
+1684 2 line 37 38
+1685 2 line 38 39
+1686 2 line 39 40
+1687 2 line 40 41
+1688 2 line 41 42
+1689 2 line 42 43
+1690 2 line 43 44
+1691 2 line 44 45
+1692 2 line 45 46
+1693 2 line 46 47
+1694 2 line 47 48
+1695 2 line 48 49
+1696 2 line 49 50
+1697 2 line 50 51
+1698 2 line 51 52
+1699 2 line 52 53
+1700 2 line 53 54
+1701 2 line 54 55
+1702 2 line 55 56
+1703 2 line 56 57
+1704 2 line 57 58
+1705 2 line 58 59
+1706 2 line 59 60
+1707 2 line 60 61
+1708 2 line 61 62
+1709 2 line 62 63
+1710 2 line 63 64
+1711 2 line 64 65
+1712 2 line 65 66
+1713 3 line 133 67
+1714 3 line 132 133
+1715 3 line 131 132
+1716 3 line 34 68
+1717 3 line 68 69
+1718 3 line 69 70
+1719 3 line 70 71
+1720 3 line 71 72
+1721 3 line 72 73
+1722 3 line 73 74
+1723 3 line 74 75
+1724 3 line 75 76
+1725 3 line 76 77
+1726 3 line 77 78
+1727 3 line 78 79
+1728 3 line 79 80
+1729 3 line 80 81
+1730 3 line 81 82
+1731 3 line 82 83
+1732 3 line 83 84
+1733 3 line 84 85
+1734 3 line 85 86
+1735 3 line 86 87
+1736 3 line 87 88
+1737 3 line 88 89
+1738 3 line 89 90
+1739 3 line 90 91
+1740 3 line 91 92
+1741 3 line 92 93
+1742 3 line 93 94
+1743 3 line 94 95
+1744 3 line 95 96
+1745 3 line 96 97
+1746 3 line 97 98
+1747 3 line 98 99
+1748 3 line 99 100
+1749 3 line 100 101
+1750 3 line 101 102
+1751 3 line 102 103
+1752 3 line 103 104
+1753 3 line 104 105
+1754 3 line 105 106
+1755 3 line 106 107
+1756 3 line 107 108
+1757 3 line 108 109
+1758 3 line 109 110
+1759 3 line 110 111
+1760 3 line 111 112
+1761 3 line 112 113
+1762 3 line 113 114
+1763 3 line 114 115
+1764 3 line 115 116
+1765 3 line 116 117
+1766 3 line 117 118
+1767 3 line 118 119
+1768 3 line 119 120
+1769 3 line 120 121
+1770 3 line 121 122
+1771 3 line 122 123
+1772 3 line 123 124
+1773 3 line 124 125
+1774 3 line 125 126
+1775 3 line 126 127
+1776 3 line 127 128
+1777 3 line 128 129
+1778 3 line 129 130
+1779 3 line 130 131
+1780 4 line 185 181
+1781 4 line 134 182
+1782 4 line 182 183
+1783 4 line 183 184
+1784 4 line 184 185
+1785 4 line 181 187
+1786 4 line 187 188
+1787 4 line 188 189
+1788 4 line 189 190
+1789 4 line 190 191
+1790 4 line 191 192
+1791 4 line 192 186
+1792 4 line 186 193
+1793 4 line 193 194
+1794 4 line 194 1
+1795 4 line 1 3
+1796 4 line 3 4
+1797 4 line 4 5
+1798 4 line 5 6
+1799 4 line 6 7
+1800 4 line 7 8
+1801 4 line 8 2
+1802 4 line 2 10
+1803 4 line 10 11
+1804 4 line 11 12
+1805 4 line 12 13
+1806 4 line 13 14
+1807 4 line 14 15
+1808 4 line 15 16
+1809 4 line 16 17
+1810 4 line 17 18
+1811 4 line 18 19
+1812 4 line 19 20
+1813 4 line 20 21
+1814 4 line 21 22
+1815 4 line 22 23
+1816 4 line 23 24
+1817 4 line 24 25
+1818 4 line 25 26
+1819 4 line 26 27
+1820 4 line 27 28
+1821 4 line 28 29
+1822 4 line 29 30
+1823 4 line 30 31
+1824 4 line 31 32
+1825 4 line 32 33
+1826 4 line 33 9
+
--- /dev/null
+/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */\r
+/* Author: David Neckels, Boulder Colorado 2007 */\r
+/* step-12.cc,v 1.33 2005/08/08 16:41:40 wolf Exp */\r
+/* Version: Version-5-2-0 */\r
+/* */\r
+/* Copyright (C) 2001, 2002, 2003, 2004, 2005 by the deal.II authors */\r
+/* */\r
+/* This file is subject to QPL and may not be distributed */\r
+/* without copyright and license information. Please refer */\r
+/* to the file deal.II/doc/license.html for the text and */\r
+/* further information on this license. */\r
+\r
+ // This program solves the Euler equations\r
+ // of gas dynamics for a given configuration\r
+ // file. It uses a standard Galerkin approach\r
+ // with weakly applied boundary conditions.\r
+ \r
+ // @sect3{Include files}\r
+\r
+ // Aztecoo require mpi (even though we run on only\r
+ // one processor in this example).\r
+#include <mpi.h>\r
+\r
+ // Here we have the necessary TRILINOS includes.\r
+ //\r
+ // Epetra is the basic trilinos vector/matrix library.\r
+#include <Epetra_MpiComm.h>\r
+#include <Epetra_Map.h>\r
+#include <Epetra_CrsGraph.h>\r
+#include <Epetra_CrsMatrix.h>\r
+#include <Epetra_Vector.h>\r
+ // Teuchos is a Trilinos utility library that is used\r
+ // to set parameters within the Aztec solver library.\r
+#include "Teuchos_ParameterList.hpp"\r
+ // Aztec is the iterative solver library.\r
+#include <AztecOO.h>\r
+#include <AztecOO_Operator.h>\r
+#define HAVE_IFPACK_TEUCHOS\r
+#include <Ifpack.h>\r
+\r
+ // Amesos is a direct solver package within Trilinos.\r
+#include <Amesos.h>\r
+ // Sacado is the automatic differentiation package, which\r
+ // is used to find the jacobian for a fully implicit Newton\r
+ // iteration.\r
+#include <Sacado.hpp>\r
+\r
+ // A standard set of dealii includes. Nothing special to\r
+ // comment on here.\r
+#include <base/quadrature_lib.h>\r
+#include <base/function.h>\r
+#include <base/parameter_handler.h>\r
+#include <base/function_parser.h>\r
+\r
+#include <lac/vector.h>\r
+#include <lac/sparse_matrix.h>\r
+#include <lac/vector_memory.h>\r
+\r
+#include <grid/tria.h>\r
+#include <grid/grid_generator.h>\r
+#include <grid/grid_out.h>\r
+#include <grid/grid_refinement.h>\r
+#include <grid/tria_accessor.h>\r
+#include <grid/tria_iterator.h>\r
+#include <grid/grid_in.h>\r
+\r
+#include <fe/fe_values.h>\r
+#include <fe/fe_system.h>\r
+\r
+#include <dofs/dof_handler.h>\r
+#include <dofs/dof_accessor.h>\r
+#include <dofs/dof_tools.h>\r
+\r
+#include <numerics/data_out.h>\r
+#include <numerics/vectors.h>\r
+#include <numerics/solution_transfer.h>\r
+\r
+#include <fe/mapping_q1.h>\r
+#include <fe/fe_q.h>\r
+#include <numerics/derivative_approximation.h>\r
+ // And this again is C++:\r
+#include <iostream>\r
+#include <fstream>\r
+#include <vector>\r
+\r
+ // Introduce the dealii library into the current namespace.\r
+using namespace dealii;\r
+\r
+ // We define a shorter name for the automatic differentiation\r
+ // type.\r
+typedef Sacado::Fad::DFad<double> fad_double;\r
+typedef unsigned int UInt;\r
+ // The Epetra library requires a 'communicator', which describes\r
+ // the layout of a parallel (or serial) set of processors.\r
+Epetra_MpiComm *Comm;\r
+\r
+ //@sect3{Flux function definition}\r
+ // Here we define the flux function for this system of conservation\r
+ // laws. Note: it would be terribly difficult to use this example\r
+ // to solve some other system of conservation laws.\r
+ //\r
+ // We define the number of components in the system. Euler's has\r
+ // one entry for momenta in each spatial direction, plus the energy\r
+ // and density components.\r
+#define N_COMP (2 + DIMENSION)\r
+ // Define a handle to the density and energy indices. We have arrange\r
+ // the momenta to be first, then density, and, lastly, energy.\r
+#define DENS_IDX DIMENSION\r
+#define ENERGY_IDX (DIMENSION+1)\r
+\r
+ // The gas constant. This value is representative of air.\r
+const double GAMMA = 1.4;\r
+ // We define the flux functions as one large matrix. Each row of this\r
+ // matrix represents a scalar conservation law for the component in\r
+ // that row. We template the numerical type of the flux function\r
+ // so that we may use the automatic differentiation type here.\r
+ // The flux functions are defined in terms of the\r
+ // conserved variables $\rho w_0, \dots, \rho w_{d-1}, \rho, E$,\r
+ // so they do not look exactly like the Euler equations one is\r
+ // used to seeing. We evaluate the flux at a single quadrature\r
+ // point.\r
+template <typename number, int dim>\r
+void Flux(std::vector<std::vector<number> > &flux, \r
+ const Point<dim> &point, \r
+ const std::vector<number> &W)\r
+{\r
+\r
+ // Pressure is a dependent variable: $p = \r
+ // (\gamma - 1)(E-\frac{1}{2} \rho |v|^2)$.\r
+ number rho_normVsqr;\r
+ for (int d0 = 0; d0 < dim; d0++) rho_normVsqr += W[d0]*W[d0];\r
+ // Since W are $\rho v$, we get a $\rho^2$ in the\r
+ // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$.\r
+ rho_normVsqr /= W[DENS_IDX];\r
+\r
+ number pressure = (GAMMA-1.0)*(W[ENERGY_IDX] - number(0.5)*(rho_normVsqr));\r
+\r
+ // We compute the momentum terms. We divide by the\r
+ // density here to get $v_i \rho v_j$\r
+ for (int d = 0; d < dim; d++) {\r
+ for (int d1 = 0; d1 < dim; d1++) {\r
+ flux[d][d1] = W[d]*W[d1]/W[DENS_IDX];\r
+ }\r
+ // The pressure contribution, along the diagonal:\r
+ flux[d][d] += pressure;\r
+ // Advection/conservation of density:\r
+ flux[DENS_IDX][d] = W[d]; \r
+ // And, lastly, conservation of energy.\r
+ flux[ENERGY_IDX][d] = W[d]/W[DENS_IDX]*\r
+ (W[ENERGY_IDX] + pressure); // energy\r
+ }\r
+}\r
+\r
+ // On the boundaries of the domain and across `hanging nodes` we use\r
+ // a numerical flux function to enforce boundary conditions. This routine\r
+ // is the basic Lax-Friedrich's flux with a stabilization parameter\r
+ // $\alpha$.\r
+template <typename number, int dim>\r
+void LFNumFlux(\r
+ std::vector<std::vector<fad_double> > &nflux,\r
+ const std::vector<Point<dim> > &points, \r
+ const std::vector<Point<dim> > &normals,\r
+ const std::vector<std::vector<number> > &Wplus,\r
+ const std::vector<std::vector<number> > &Wminus,\r
+ double alpha)\r
+{\r
+ int n_q_points = points.size();\r
+\r
+ // We evaluate the flux at each of the quadrature points.\r
+ for (int q = 0; q < n_q_points; q++) {\r
+ std::vector<std::vector<fad_double> > iflux(N_COMP,\r
+ std::vector<fad_double>(dim, 0));\r
+ std::vector<std::vector<fad_double> > oflux(N_COMP,\r
+ std::vector<fad_double>(dim, 0));\r
+\r
+ Flux<number, dim>(iflux, points[q], Wplus[q]);\r
+ Flux<number, dim>(oflux, points[q], Wminus[q]);\r
+\r
+ for (int di = 0; di < N_COMP; di++) {\r
+ nflux[q][di] = 0;\r
+ for (int d = 0; d < dim; d++) {\r
+ nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d);\r
+ }\r
+ nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);\r
+ }\r
+ }\r
+\r
+}\r
+\r
+ // @sect3{Initial and side condition parsing}\r
+ // For the initial condition we use the expression parser function\r
+ // object.\r
+template <int dim>\r
+class InitialCondition : public FunctionParser<dim> \r
+{\r
+ public:\r
+ InitialCondition ();\r
+ \r
+ // This function should be called after parsing, but before using\r
+ // the object. It formalizes the expressions and initializes the\r
+ // function parser with the appropriate expressions.\r
+ void Init();\r
+\r
+ // During parsing we call this function as the initial condition\r
+ // for one of the $\mathbf{w}$ variables is encountered.\r
+\r
+ void set_ic(int _row, std::string &expr) {\r
+ expressions[_row] = expr;\r
+ }\r
+\r
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,\r
+ std::vector<Vector<double> > &value_list) const;\r
+ private:\r
+ std::vector<std::string> expressions;\r
+};\r
+\r
+template <int dim>\r
+InitialCondition<dim>::InitialCondition () :\r
+ FunctionParser<dim> (N_COMP),\r
+ expressions(N_COMP, "0.0")\r
+{}\r
+\r
+ // Here we set up x,y,z as the variables that one should use in the input\r
+ // deck to describe their initial condition.\r
+template<int dim>\r
+void InitialCondition<dim>::Init() {\r
+ std::map<std::string, double> constants;\r
+ constants["M_PI"] = M_PI;\r
+ std::string variables = (dim == 2 ? "x,y" : "x,y,z");\r
+\r
+ FunctionParser<dim>::initialize(variables, expressions, constants);\r
+\r
+}\r
+\r
+template <int dim>\r
+void InitialCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,\r
+ std::vector<Vector<double> > &value_list) const \r
+{\r
+ const unsigned int n_points = points.size();\r
+\r
+ Assert (value_list.size() == n_points, \r
+ ExcDimensionMismatch (value_list.size(), n_points));\r
+\r
+ for (unsigned int p=0; p<n_points; ++p)\r
+ InitialCondition<dim>::vector_value (points[p],\r
+ value_list[p]);\r
+}\r
+\r
+ // As above, we use the expression function parser for boundary conditions.\r
+template <int dim>\r
+class SideCondition : public FunctionParser<dim> \r
+{\r
+ public:\r
+ SideCondition (int ncomp);\r
+ ~SideCondition ();\r
+\r
+ // As above.\r
+ void Init();\r
+ // As above.\r
+ void set_coeff_row(int _row_n, std::string &expr);\r
+ \r
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,\r
+ std::vector<Vector<double> > &value_list) const;\r
+ private:\r
+ std::vector<std::string> expressions;\r
+};\r
+\r
+template <int dim>\r
+SideCondition<dim>::SideCondition (int ncomp) :\r
+ FunctionParser<dim> (ncomp),\r
+ expressions(ncomp, "0.0")\r
+{\r
+}\r
+template <int dim>\r
+void SideCondition<dim>::set_coeff_row (int _row_n, std::string &expr) \r
+{\r
+ expressions[_row_n] = expr;\r
+}\r
+\r
+template <int dim>\r
+void SideCondition<dim>::Init() {\r
+ std::map<std::string, double> constants;\r
+ constants["M_PI"] = M_PI;\r
+ std::string variables = (dim == 2 ? "x,y" : "x,y,z");\r
+\r
+ FunctionParser<dim>::initialize(variables, expressions, constants);\r
+\r
+}\r
+\r
+template <int dim>\r
+SideCondition<dim>::~SideCondition () \r
+{\r
+}\r
+\r
+template <int dim>\r
+void SideCondition<dim>::vector_value_list (const std::vector<Point<dim> > &points,\r
+ std::vector<Vector<double> > &value_list) const \r
+{\r
+ const unsigned int n_points = points.size();\r
+\r
+ Assert (value_list.size() == n_points, \r
+ ExcDimensionMismatch (value_list.size(), n_points));\r
+\r
+ for (unsigned int p=0; p<n_points; ++p)\r
+ SideCondition<dim>::vector_value (points[p],\r
+ value_list[p]);\r
+}\r
+ //@sect3{Conservation Law class}\r
+ // Here we define a Conservation Law class that helps group\r
+ // operations and data for our Euler equations into a manageable\r
+ // entity. Functions will be described as their definitions appear.\r
+template <int dim>\r
+class ConsLaw\r
+{\r
+ public:\r
+ ConsLaw ();\r
+ ~ConsLaw ();\r
+\r
+ void run ();\r
+ void declare_parameters();\r
+ void load_parameters(const char *);\r
+ \r
+ private:\r
+ void build_fe();\r
+ void setup_system ();\r
+ void initialize_system ();\r
+ void assemble_system (double &res_norm);\r
+ void solve (Vector<double> &solution, int &, double &);\r
+ void refine_grid ();\r
+ void output_results (const unsigned int cycle) const;\r
+ void initialize();\r
+ void zero_matrix();\r
+ void estimate();\r
+ void postprocess();\r
+ void compute_predictor();\r
+ \r
+ Triangulation<dim> triangulation;\r
+ const MappingQ1<dim> mapping;\r
+ \r
+ \r
+ FESystem<dim> *fe_ptr;\r
+\r
+ DoFHandler<dim> dof_handler;\r
+\r
+ SparsityPattern sparsity_pattern;\r
+ const QGauss<dim> quadrature;\r
+ const QGauss<dim-1> face_quadrature;\r
+ \r
+ // The actual solution to the Euler equation\r
+ Vector<double> solution;\r
+ // The current value of the solution during the Newton iteration\r
+ Vector<double> nlsolution;\r
+ // An estimate of the next time value; used for adaptivity and as a\r
+ // guess for the next Newton iteration.\r
+ Vector<double> predictor;\r
+ // Values after post-processing (used to output the physical variables).\r
+ Vector<double> ppsolution;\r
+ // The solution to the linear problem during the Newton iteration\r
+ Vector<double> dsolution;\r
+ Vector<double> right_hand_side;\r
+ \r
+ public:\r
+\r
+ void assemble_cell_term(const FEValues<dim>& fe_v,\r
+ std::vector<unsigned int> &dofs,\r
+ unsigned int cell_no\r
+ );\r
+ \r
+ void assemble_face_term(\r
+ int face_no,\r
+ const FEFaceValuesBase<dim>& fe_v,\r
+ const FEFaceValuesBase<dim>& fe_v_neighbor,\r
+ std::vector<unsigned int> &dofs,\r
+ std::vector<unsigned int> &dofs_neighbor,\r
+ int boundary = -1\r
+ );\r
+\r
+ unsigned int get_n_components() const { return N_COMP;}\r
+\r
+ private:\r
+ // T = current time, dT = time step, TF = final time.\r
+ double T, dT, TF;\r
+ double face_diameter;\r
+ double cell_diameter;\r
+ // An object to handle parsing the input deck.\r
+ ParameterHandler prm;\r
+ // Name of the mesh to read in.\r
+ string mesh;\r
+ InitialCondition<dim> ic;\r
+\r
+ // Enums for the various supported boundary conditions.\r
+ typedef enum {INFLOW_BC = 1, OUTFLOW_BC=2, NO_PENETRATION_BC=3, PRESSURE_BC=4} bc_type;\r
+\r
+ // For each boundary we store a map from boundary # to the type\r
+ // of boundary condition. If the boundary condition is prescribed,\r
+ // we store a pointer to a function object that will hold the expression\r
+ // for that boundary condition.\r
+ typedef typename std::map<unsigned int, std::pair<std::vector<bc_type>, Function<dim>*> > bdry_map_type;\r
+ bdry_map_type bdry_map;\r
+\r
+ void add_boundary(unsigned int bd, std::vector<bc_type>& flags, Function<dim> *bf);\r
+\r
+ // An object to store parameter information about the Aztec solver.\r
+ typedef struct {\r
+ int LIN_OUTPUT;\r
+ typedef enum { GMRES = 0, DIRECT = 1} solver_type;\r
+ solver_type SOLVER;\r
+ typedef enum { QUIET = 0, VERBOSE = 1 } output_type;\r
+ output_type OUTPUT;\r
+ // Linear residual tolerance.\r
+ double RES;\r
+ int MAX_ITERS;\r
+ // We use the ILUT preconditioner. This is similar\r
+ // to the ILU. FILL is the number of extra entries\r
+ // to add when forming the ILU decomposition.\r
+ double ILUT_FILL;\r
+ // When forming the preconditioner, for certain problems\r
+ // bad conditioning (or just bad luck) can cause the\r
+ // preconditioner to be very poorly conditioned. Hence\r
+ // it can help to add diagonal perturbations to the\r
+ // original matrix and form the preconditioner for this\r
+ // slightly better matrix. ATOL is an absolute perturbation\r
+ // that is added to the diagonal before forming the\r
+ // prec, and RTOL is a scaling factor $rtol >= 1$.\r
+ double ILUT_ATOL;\r
+ double ILUT_RTOL;\r
+ // The ILUT will drop any values that have magnitude less\r
+ // than this value. This is a way to\r
+ // manage the amount of memory used by this preconditioner.\r
+ double ILUT_DROP;\r
+ } solver_params_type;\r
+\r
+ solver_params_type solver_params;\r
+\r
+ // Some refinement parameters.\r
+ typedef struct {\r
+ typedef enum { NONE = 0, FIXED_NUMBER = 1, SHOCK = 2} refine_type;\r
+ double high_frac;\r
+ double low_frac;\r
+ refine_type refine;\r
+ double high_frac_sav;\r
+ double max_cells;\r
+ double shock_val;\r
+ double shock_levels;\r
+ } refinement_params_type;\r
+\r
+ refinement_params_type refinement_params;\r
+\r
+ // The user can set the stabilization parameter $\alpha$ \r
+ // in the Lax-Friedrich's flux. \r
+ typedef struct {\r
+ typedef enum {CONSTANT=1,MESH=2} LF_stab_type;\r
+ LF_stab_type LF_stab;\r
+ double LF_stab_value;\r
+ } flux_params_type;\r
+\r
+ flux_params_type flux_params;\r
+\r
+ bool is_stationary;\r
+\r
+ // Power for the mesh stabilization term.\r
+ double diffusion_power;\r
+ double gravity;\r
+ // If true, we output the squared gradient of the\r
+ // density instead of density. Using this one can\r
+ // create shock plots.\r
+ bool schlieren_plot;\r
+ // How often to create an output file.\r
+ double output_step;\r
+\r
+ Epetra_CrsMatrix *Matrix;\r
+ Epetra_Map *Map;\r
+ Vector<double> indicator;\r
+ \r
+ // Crank-Nicolson value\r
+ const double theta; \r
+\r
+};\r
+\r
+\r
+ // Asign a row of the conservation law a specified\r
+ // boundary type and (possibly) function.\r
+template <int dim>\r
+void ConsLaw<dim>::add_boundary(unsigned int bd,\r
+ std::vector<bc_type> &flags, Function<dim> *bf) {\r
+\r
+ std::pair<std::vector<bc_type>, Function<dim> *> entry(flags, bf);\r
+ bdry_map[bd] = entry;\r
+}\r
+\r
+\r
+ // Apply the initialial condition. Simultaneously\r
+ // initialize the non-linear solution.\r
+template <int dim>\r
+void ConsLaw<dim>::initialize() {\r
+ VectorTools::interpolate(dof_handler,\r
+ ic, solution);\r
+ VectorTools::interpolate(dof_handler,\r
+ ic, nlsolution);\r
+\r
+}\r
+\r
+ // @sect3{Assembly}\r
+ // @sect4{Function: assemble_cell_term}\r
+ //\r
+ // Assembles the cell term, adding minus the residual\r
+ // to the right hand side, and adding in the Jacobian\r
+ // contributions.\r
+template <int dim>\r
+void ConsLaw<dim>::assemble_cell_term(\r
+ const FEValues<dim> &fe_v,\r
+ std::vector<unsigned int> &dofs,\r
+ unsigned int cell_no\r
+ ) \r
+{\r
+ // The residual for each row (i) will be accumulating \r
+ // into this fad variable. At the end of the assembly\r
+ // for this row, we will query for the sensitivities\r
+ // to this variable and add them into the Jacobian.\r
+ fad_double F_i;\r
+ int dofs_per_cell = fe_v.dofs_per_cell;\r
+ int n_q_points = fe_v.n_quadrature_points;\r
+\r
+ // We will define the dofs on this cell in these fad variables.\r
+ std::vector<fad_double> DOF(dofs_per_cell);\r
+\r
+ // Values of the conservative variables at the quadrature points.\r
+ std::vector<std::vector<fad_double > > W (n_q_points,\r
+ std::vector<fad_double >(get_n_components()));\r
+\r
+ // Values at the last time step of the conservative variables.\r
+ // Note that these do not use fad variables, since they do\r
+ // not depend on the 'variables to be sought'=DOFS.\r
+ std::vector<std::vector<double > > Wl (n_q_points,\r
+ std::vector<double >(get_n_components()));\r
+\r
+ // Here we will hold the averaged values of the conservative\r
+ // variables that we will linearize around (cn=Crank Nicholson).\r
+ std::vector<std::vector<fad_double > > Wcn (n_q_points,\r
+ std::vector<fad_double >(get_n_components()));\r
+\r
+ // Gradients of the current variables. It is a\r
+ // bit of a shame that we have to compute these; we almost don't.\r
+ // The nice thing about a simple conservation law is that the\r
+ // the flux doesn't generally involve any gradients. We do\r
+ // need these, however, for the diffusion stabilization. \r
+ std::vector<std::vector<std::vector<fad_double> > > Wgrads (n_q_points,\r
+ std::vector<std::vector<fad_double> >(get_n_components(),\r
+ std::vector<fad_double>(dim)));\r
+\r
+\r
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
+\r
+ \r
+ // Here is the magical point where we declare a subset\r
+ // of the fad variables as degrees of freedom. All \r
+ // calculations that reference these variables (either\r
+ // directly or indirectly) will accumulate sensitivies\r
+ // with respect to these dofs.\r
+ for (int in = 0; in < dofs_per_cell; in++) {\r
+ DOF[in] = nlsolution(dofs[in]);\r
+ DOF[in].diff(in, dofs_per_cell);\r
+ }\r
+\r
+ // Here we compute the shape function values and gradients\r
+ // at the quadrature points. Ideally, we could call into \r
+ // something like get_function_values, get_function_grads,\r
+ // but since we don't want to make the entire solution vector\r
+ // fad types, only the local cell variables, we explicitly\r
+ // code this loop;\r
+ for (int q = 0; q < n_q_points; q++) {\r
+ for (int di = 0; di < get_n_components(); di++) {\r
+ W[q][di] = 0;\r
+ Wl[q][di] = 0;\r
+ Wcn[q][di] = 0;\r
+ for (int d = 0; d < dim; d++) {\r
+ Wgrads[q][di][d] = 0;\r
+ }\r
+ }\r
+ for (int sf = 0; sf < dofs_per_cell; sf++) {\r
+ int di = fe_v.get_fe().system_to_component_index(sf).first;\r
+ W[q][di] +=\r
+ DOF[sf]*fe_v.shape_value_component(sf, q, di);\r
+ Wl[q][di] +=\r
+ solution(dofs[sf])*fe_v.shape_value_component(sf, q, di);\r
+ Wcn[q][di] +=\r
+ (theta*DOF[sf]+(1-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);\r
+\r
+ for (int d = 0; d < dim; d++) {\r
+ Wgrads[q][di][d] += DOF[sf]*\r
+ fe_v.shape_grad_component(sf, q, di)[d];\r
+ } // for d\r
+\r
+ }\r
+\r
+ } // for q\r
+\r
+ // Gather the flux values for all components at\r
+ // all of the quadrature points. This also\r
+ // computes the matrix of sensitivities. Perhaps\r
+ // this could be done in a better way, since this\r
+ // could be a rather large object, but for now it \r
+ // seems to work just fine.\r
+ std::vector<std::vector<std::vector<fad_double> > > flux(n_q_points, \r
+ std::vector<std::vector<fad_double> >(get_n_components(),\r
+ std::vector<fad_double>(dim, 0)));\r
+\r
+ for (unsigned int q=0; q < n_q_points; ++q) {\r
+ Flux<fad_double, dim>(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]);\r
+ }\r
+\r
+ // We now have all of the function values/grads/fluxes,\r
+ // so perform the assembly. We have an outer loop\r
+ // through the components of the system, and an\r
+ // inner loop over the quadrature points, where we\r
+ // accumulate contributions to the ith residual.\r
+ //\r
+ // We initialy sum all contributions of the residual\r
+ // in the positive sense, so that we don't need to\r
+ // negative the Jacobian entries. Then, when we sum\r
+ // into the <code> right_hand_side </code> vector,\r
+ // we negate this residual.\r
+ for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) \r
+ {\r
+ // Find which component this dof contributes to.\r
+ const unsigned int\r
+ component_i = fe_v.get_fe().system_to_component_index(i).first;\r
+\r
+ // Initialize the fad residual to zero (removes\r
+ // any previous sensitivities.\r
+ F_i = 0;\r
+\r
+ // Loop quadrature points.\r
+ for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point) {\r
+\r
+ fad_double fdotgv = 0;\r
+\r
+ // Integrate the flux times gradient of the test function\r
+ for (int d = 0; d < dim; d++) \r
+ fdotgv += flux[point][component_i][d]*fe_v.shape_grad_component(i, point, component_i)[d];\r
+ \r
+ F_i -= fdotgv*JxW[point];\r
+\r
+ // The mass term (if the simulation is non-stationary).\r
+ fad_double delta_t= 1.0/dT*(W[point][component_i] - Wl[point][component_i]);\r
+ if (!is_stationary) F_i += delta_t*\r
+ fe_v.shape_value_component(i, point, component_i)*JxW[point];\r
+\r
+ // Stabilization (cell wise diffusion)\r
+ fad_double guv = 0;\r
+ for (int d = 0; d < dim; d++) {\r
+ guv += fe_v.shape_grad_component(i, point, component_i)[d]*\r
+ Wgrads[point][component_i][d];\r
+ }\r
+\r
+ F_i += 1.0*std::pow(cell_diameter, diffusion_power)*guv*JxW[point];\r
+ \r
+ // The gravity component only enters into the energy \r
+ // equation and into the vertical component of the \r
+ // velocity.\r
+ if (component_i == dim - 1) {\r
+ F_i += gravity*Wcn[point][DENS_IDX]*fe_v.shape_value_component(i,point, component_i)*JxW[point];\r
+ } else if (component_i == ENERGY_IDX) {\r
+ F_i += gravity*Wcn[point][DENS_IDX]*Wcn[point][dim-1]*\r
+ fe_v.shape_value_component(i,point, component_i)*JxW[point];\r
+ }\r
+ } // for q\r
+\r
+ // Here we gain access to the array of sensitivities\r
+ // of the residual. We then sum these into the\r
+ // Epetra matrix.\r
+ double *values = &(F_i.fastAccessDx(0));\r
+ Matrix->SumIntoGlobalValues(dofs[i],\r
+ dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));\r
+ \r
+ // Add minus the residual to the right hand side.\r
+ right_hand_side(dofs[i]) -= F_i.val();\r
+\r
+ } // for i\r
+}\r
+ // @sect4{Function: assemble_face_term}\r
+ // These are either\r
+ // boundary terms or terms across differing \r
+ // levels of refinement. In the first case,\r
+ // fe_v==fe_v_neighbor and dofs==dofs_neighbor.\r
+ // The int boundary < 0 if not at a boundary,\r
+ // otherwise it is the boundary indicator.\r
+template <int dim>\r
+void ConsLaw<dim>::assemble_face_term(\r
+ int face_no,\r
+ const FEFaceValuesBase<dim>& fe_v,\r
+ const FEFaceValuesBase<dim>& fe_v_neighbor, \r
+ std::vector<unsigned int> &dofs,\r
+ std::vector<unsigned int> &dofs_neighbor,\r
+ int boundary\r
+ ) \r
+{\r
+ fad_double F_i;\r
+ const unsigned int n_q_points = fe_v.n_quadrature_points;\r
+ const unsigned int dofs_per_cell = fe_v.get_fe().dofs_per_cell;\r
+ const unsigned int ndofs_per_cell = fe_v_neighbor.get_fe().dofs_per_cell;\r
+ Assert(dofs_per_cell == ndofs_per_cell,\r
+ ExcDimensionMismatch(dofs_per_cell, ndofs_per_cell));\r
+\r
+ // As above, the fad degrees of freedom\r
+ std::vector<fad_double> DOF(dofs_per_cell+ndofs_per_cell);\r
+\r
+ // The conservative variables for this cell,\r
+ // and for \r
+ std::vector<std::vector<fad_double > > Wplus (n_q_points,\r
+ std::vector<fad_double >(get_n_components()));\r
+ std::vector<std::vector<fad_double > > Wminus (n_q_points,\r
+ std::vector<fad_double >(get_n_components()));\r
+\r
+\r
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
+ const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();\r
+\r
+\r
+ // If we are at a boundary, then dofs_neighbor are\r
+ // the same as dofs, so we do not want to duplicate them.\r
+ // If there is a neighbor cell, then we want to include \r
+ // them.\r
+ int ndofs = (boundary < 0 ? dofs_per_cell + ndofs_per_cell : dofs_per_cell);\r
+ // Set the local DOFS.\r
+ for (int in = 0; in < dofs_per_cell; in++) {\r
+ DOF[in] = nlsolution(dofs[in]);\r
+ DOF[in].diff(in, ndofs);\r
+ }\r
+ // If present, set the neighbor dofs.\r
+ if (boundary < 0)\r
+ for (int in = 0; in < ndofs_per_cell; in++) {\r
+ DOF[in+dofs_per_cell] = nlsolution(dofs_neighbor[in]);\r
+ DOF[in+dofs_per_cell].diff(in+dofs_per_cell, ndofs);\r
+ }\r
+\r
+ // Set the values of the local conservative variables.\r
+ // Initialize all variables to zero.\r
+ for (int q = 0; q < n_q_points; q++) {\r
+ for (int di = 0; di < get_n_components(); di++) {\r
+ Wplus[q][di] = 0;\r
+ Wminus[q][di] = 0;\r
+ }\r
+ for (int sf = 0; sf < dofs_per_cell; sf++) {\r
+ int di = fe_v.get_fe().system_to_component_index(sf).first;\r
+ Wplus[q][di] +=\r
+ (theta*DOF[sf]+(1.0-theta)*solution(dofs[sf]))*fe_v.shape_value_component(sf, q, di);\r
+ }\r
+\r
+\r
+ // If there is a cell across, then initialize\r
+ // the exterior trace as a function of the other\r
+ // cell degrees of freedom.\r
+ if (boundary < 0) {\r
+ for (int sf = 0; sf < ndofs_per_cell; sf++) {\r
+ int di = fe_v_neighbor.get_fe().system_to_component_index(sf).first;\r
+ Wminus[q][di] +=\r
+ (theta*DOF[sf+dofs_per_cell]+(1.0-theta)*solution(dofs_neighbor[sf]))*\r
+ fe_v_neighbor.shape_value_component(sf, q, di);\r
+ }\r
+ } \r
+ } // for q\r
+\r
+ // If this is a boundary, then the values of $W^-$ will\r
+ // be either functions of $W^+$, or they will be prescribed.\r
+ // This switch sets them appropriately. Since we are\r
+ // using fad variables here, sensitivities will be updated \r
+ // appropriately. These sensitivities would be tremendously\r
+ // difficult to manage without fad!!!\r
+ if (boundary >= 0) {\r
+ // Get the boundary descriptor.\r
+ typename bdry_map_type::iterator bme = bdry_map.find(boundary);\r
+ assert(bme != bdry_map.end());\r
+\r
+ // Evaluate the function object. This is a bit\r
+ // tricky; a given boundary might have both prescribed\r
+ // and implicit values. If a particular component is not\r
+ // prescribed, the values evaluate to zero and are\r
+ // ignored, below.\r
+ std::vector<Vector<double> > bvals(n_q_points, Vector<double>(N_COMP));\r
+ bme->second.second->vector_value_list(fe_v.get_quadrature_points(), bvals);\r
+\r
+ // We loop the quadrature points, and we treat each\r
+ // component individualy.\r
+ for (int q = 0; q < n_q_points; q++) {\r
+ for (int di = 0; di < get_n_components(); di++) {\r
+\r
+ // An inflow/dirichlet type of boundary condition\r
+ if (bme->second.first[di] == INFLOW_BC) {\r
+ Wminus[q][di] = bvals[q](di);\r
+ } else if (bme->second.first[di] == PRESSURE_BC) {\r
+ // A prescribed pressure boundary condition. This boundary\r
+ // condition is complicated by the fact that even though\r
+ // the pressure is prescribed, we really are setting\r
+ // the energy index here, which will depend on velocity\r
+ // and pressure. So even though this seems like a dirichlet\r
+ // type boundary condition, we get sensitivities of\r
+ // energy to velocity and density (unless these\r
+ // are also prescribed.\r
+ fad_double rho_vel_sqr = 0;\r
+ fad_double dens;\r
+ \r
+ dens = bme->second.first[DENS_IDX] == INFLOW_BC ? bvals[q](DENS_IDX) :\r
+ Wplus[q][DENS_IDX];\r
+\r
+ for (int d=0; d < dim; d++) {\r
+ if (bme->second.first[d] == INFLOW_BC)\r
+ rho_vel_sqr += bvals[q](d)*bvals[q](d);\r
+ else\r
+ rho_vel_sqr += Wplus[q][d]*Wplus[q][d];\r
+ }\r
+ rho_vel_sqr /= dens;\r
+ // Finally set the energy value as determined by the\r
+ // prescribed pressure and the other variables.\r
+ Wminus[q][di] = bvals[q](di)/(GAMMA-1.0) +\r
+ 0.5*rho_vel_sqr;\r
+\r
+ } else if (bme->second.first[di] == OUTFLOW_BC) {\r
+ // A free/outflow boundary, very simple.\r
+ Wminus[q][di] = Wplus[q][di];\r
+\r
+ } else { \r
+ // We must be at a no-penetration boundary. We\r
+ // prescribe the velocity (we are dealing with a\r
+ // particular component here so that the average\r
+ // of the velocities is orthogonal to the surface\r
+ // normal. This creates sensitivies of across\r
+ // the velocity components.\r
+ fad_double vdotn = 0;\r
+ for (int d = 0; d < dim; d++) {\r
+ vdotn += Wplus[q][d]*normals[q](d);\r
+ }\r
+\r
+ Wminus[q][di] = Wplus[q][di] - 2.0*vdotn*normals[q](di);\r
+ }\r
+ }\r
+ } // for q\r
+ } // b>= 0\r
+ \r
+ // Determine the Lax-Friedrich's stability parameter,\r
+ // and evaluate the numerical flux function at the quadrature points\r
+ std::vector<std::vector<fad_double> > nflux(n_q_points, std::vector<fad_double>(get_n_components(), 0));\r
+ double alpha = 1;\r
+\r
+ switch(flux_params.LF_stab) {\r
+ case flux_params_type::CONSTANT:\r
+ alpha = flux_params.LF_stab_value;\r
+ break;\r
+ case flux_params_type::MESH:\r
+ alpha = face_diameter/(2.0*dT);\r
+ break;\r
+ }\r
+\r
+ LFNumFlux<fad_double, dim>(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,\r
+ alpha);\r
+\r
+ // Now assemble the face term\r
+ for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {\r
+ if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;\r
+ F_i = 0;\r
+ for (unsigned int point=0; point<n_q_points; ++point)\r
+ {\r
+ const unsigned int\r
+ component_i = fe_v.get_fe().system_to_component_index(i).first;\r
+\r
+ F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*JxW[point];\r
+\r
+ } \r
+\r
+ // Retrieve a pointer to the jacobian.\r
+ double *values = &(F_i.fastAccessDx(0));\r
+\r
+ // Honestly, I forget why this can happen, but \r
+ // for some reason it can!!\r
+ if (!values) continue;\r
+\r
+ // Update the matrix. Depending on whether there\r
+ // is/isn't a neighboring cell, we add more/less\r
+ // entries.\r
+ Matrix->SumIntoGlobalValues(dofs[i],\r
+ dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));\r
+ if (boundary < 0) {\r
+ Matrix->SumIntoGlobalValues(dofs[i],\r
+ dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&dofs_neighbor[0]));\r
+ }\r
+\r
+ // And add into the residual\r
+ right_hand_side(dofs[i]) -= F_i.val();\r
+ } \r
+\r
+}\r
+ // @sect4{Assembling the whole system}\r
+ // Now we put all of the assembly pieces together\r
+ // in a routine that dispatches the correct\r
+ // piece for each cell/face. We keep track of\r
+ // the norm of the resdual for the Newton iteration.\r
+template <int dim>\r
+void ConsLaw<dim>::assemble_system (double &res_norm) \r
+{\r
+ FESystem<dim> &fe = *fe_ptr;\r
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
+\r
+ // We track the dofs on this cell and (if necessary)\r
+ // the adjacent cell.\r
+ std::vector<unsigned int> dofs (dofs_per_cell);\r
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);\r
+\r
+ // First we create the\r
+ // ``UpdateFlags'' for the\r
+ // ``FEValues'' and the\r
+ // ``FEFaceValues'' objects.\r
+ UpdateFlags update_flags = update_values\r
+ | update_gradients\r
+ | update_q_points\r
+ | update_JxW_values;\r
+\r
+ // Note, that on faces we do not\r
+ // need gradients but we need\r
+ // normal vectors.\r
+ UpdateFlags face_update_flags = update_values\r
+ | update_q_points\r
+ | update_JxW_values\r
+ | update_normal_vectors;\r
+ \r
+ // On the neighboring cell we only\r
+ // need the shape values. Given a\r
+ // specific face, the quadrature\r
+ // points and `JxW values' are the\r
+ // same as for the current cells,\r
+ // the normal vectors are known to\r
+ // be the negative of the normal\r
+ // vectors of the current cell.\r
+ UpdateFlags neighbor_face_update_flags = update_values;\r
+ \r
+ // Then we create the ``FEValues''\r
+ // object. Note, that since version\r
+ // 3.2.0 of deal.II the constructor\r
+ // of this class takes a\r
+ // ``Mapping'' object as first\r
+ // argument. Although the\r
+ // constructor without ``Mapping''\r
+ // argument is still supported it\r
+ // is recommended to use the new\r
+ // constructor. This reduces the\r
+ // effect of `hidden magic' (the\r
+ // old constructor implicitely\r
+ // assumes a ``MappingQ1'' mapping)\r
+ // and makes it easier to change\r
+ // the mapping object later.\r
+ FEValues<dim> fe_v (\r
+ mapping, fe, quadrature, update_flags);\r
+ \r
+ // Similarly we create the\r
+ // ``FEFaceValues'' and\r
+ // ``FESubfaceValues'' objects for\r
+ // both, the current and the\r
+ // neighboring cell. Within the\r
+ // following nested loop over all\r
+ // cells and all faces of the cell\r
+ // they will be reinited to the\r
+ // current cell and the face (and\r
+ // subface) number.\r
+ FEFaceValues<dim> fe_v_face (\r
+ mapping, fe, face_quadrature, face_update_flags);\r
+ FESubfaceValues<dim> fe_v_subface (\r
+ mapping, fe, face_quadrature, face_update_flags);\r
+ FEFaceValues<dim> fe_v_face_neighbor (\r
+ mapping, fe, face_quadrature, neighbor_face_update_flags);\r
+ FESubfaceValues<dim> fe_v_subface_neighbor (\r
+ mapping, fe, face_quadrature, neighbor_face_update_flags);\r
+\r
+ // Furthermore we need some cell\r
+ // iterators.\r
+ typename DoFHandler<dim>::active_cell_iterator\r
+ cell = dof_handler.begin_active(),\r
+ endc = dof_handler.end();\r
+\r
+ // Now we start the loop over all\r
+ // active cells.\r
+ int fdofs_per_cell = fe_v.dofs_per_cell;\r
+ int fn_q_points = face_quadrature.n_quadrature_points;\r
+\r
+ unsigned int cell_no = 0;\r
+ for (;cell!=endc; ++cell, ++cell_no) \r
+ {\r
+ \r
+ // Now we reinit the ``FEValues''\r
+ // object for the current cell\r
+ fe_v.reinit (cell);\r
+\r
+ // Collect the local dofs and\r
+ // asssemble the cell term.\r
+ cell->get_dof_indices (dofs);\r
+\r
+ cell_diameter = cell->diameter();\r
+\r
+ assemble_cell_term(fe_v,\r
+ dofs,\r
+ cell_no);\r
+\r
+ // We use the DG style loop through faces\r
+ // to determine if we need to apply a\r
+ // 'hanging node' flux calculation or a boundary\r
+ // computation.\r
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)\r
+ {\r
+ // First we set the face\r
+ // iterator\r
+ typename DoFHandler<dim>::face_iterator face=cell->face(face_no);\r
+ face_diameter = face->diameter();\r
+ \r
+ if (face->at_boundary())\r
+ {\r
+ // We reinit the\r
+ // ``FEFaceValues''\r
+ // object to the\r
+ // current face\r
+ fe_v_face.reinit (cell, face_no);\r
+\r
+ // and assemble the\r
+ // corresponding face\r
+ // terms. We send the same\r
+ // fe_v and dofs as described\r
+ // in the assembly routine.\r
+ assemble_face_term(\r
+ face_no, fe_v_face,\r
+ fe_v_face,\r
+ dofs,\r
+ dofs,\r
+ face->boundary_indicator());\r
+ }\r
+ else\r
+ {\r
+ // Now we are not on\r
+ // the boundary of the\r
+ // domain, therefore\r
+ // there must exist a\r
+ // neighboring cell.\r
+ typename DoFHandler<dim>::cell_iterator neighbor=\r
+ cell->neighbor(face_no);;\r
+\r
+ if (face->has_children())\r
+ {\r
+ // case I: This cell refined compared to neighbor\r
+\r
+ const unsigned int neighbor2=\r
+ cell->neighbor_of_neighbor(face_no);\r
+ \r
+ \r
+ // We loop over\r
+ // subfaces\r
+ for (unsigned int subface_no=0;\r
+ subface_no<GeometryInfo<dim>::subfaces_per_face;\r
+ ++subface_no)\r
+ {\r
+ typename DoFHandler<dim>::active_cell_iterator\r
+ neighbor_child\r
+ = cell->neighbor_child_on_subface (face_no, subface_no);\r
+\r
+ face_diameter = neighbor_child->diameter(); // working on subface\r
+ \r
+ Assert (neighbor_child->face(neighbor2) == face->child(subface_no),\r
+ ExcInternalError());\r
+ Assert (!neighbor_child->has_children(), ExcInternalError());\r
+\r
+ fe_v_subface.reinit (cell, face_no, subface_no);\r
+ fe_v_face_neighbor.reinit (neighbor_child, neighbor2);\r
+ neighbor_child->get_dof_indices (dofs_neighbor);\r
+\r
+ // Assemble as if we are working with\r
+ // a DG element.\r
+ assemble_face_term(\r
+ face_no, fe_v_subface,\r
+ fe_v_face_neighbor,\r
+ dofs,\r
+ dofs_neighbor);\r
+ \r
+ }\r
+ // End of ``if\r
+ // (face->has_children())''\r
+ }\r
+ else\r
+ {\r
+ // We have no children, but \r
+ // the neighbor cell may be refine\r
+ // compared to use\r
+ neighbor->get_dof_indices (dofs_neighbor);\r
+ if (neighbor->level() != cell->level()) \r
+ {\r
+ // case II: This is refined compared to neighbor\r
+ Assert(neighbor->level() < cell->level(), ExcInternalError());\r
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=\r
+ cell->neighbor_of_coarser_neighbor(face_no);\r
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,\r
+ neighbor_subface_no=faceno_subfaceno.second;\r
+\r
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,\r
+ neighbor_subface_no)\r
+ == cell,\r
+ ExcInternalError());\r
+\r
+ // Reinit the\r
+ // appropriate\r
+ // ``FEFaceValues''\r
+ // and assemble\r
+ // the face\r
+ // terms.\r
+ fe_v_face.reinit (cell, face_no);\r
+ fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,\r
+ neighbor_subface_no);\r
+ \r
+ assemble_face_term(\r
+ face_no, fe_v_face,\r
+ fe_v_subface_neighbor,\r
+ dofs,\r
+ dofs_neighbor);\r
+\r
+ }\r
+\r
+ } \r
+ // End of ``face not at boundary'':\r
+ }\r
+ // End of loop over all faces:\r
+ } \r
+ \r
+ // End iteration through cells.\r
+ } \r
+\r
+ // Notify Epetra that the matrix is done.\r
+ Matrix->FillComplete();\r
+\r
+ // Compute the nonlinear residual.\r
+ res_norm = right_hand_side.l2_norm();\r
+ \r
+}\r
+\r
+ // Create a conservation law with some defaults.\r
+template <int dim>\r
+ConsLaw<dim>::ConsLaw ()\r
+ :\r
+ mapping (),\r
+ fe_ptr(NULL),\r
+ dof_handler (triangulation),\r
+ quadrature (2),\r
+ face_quadrature (2),\r
+ T(0),\r
+ dT(0.05),\r
+ TF(10),\r
+ Map(NULL),\r
+ Matrix(NULL),\r
+ is_stationary(false),\r
+ theta(0.5)\r
+{}\r
+\r
+ // At one time this example could work for both DG and\r
+ // continuous finite elements. The choice was made here.\r
+template <int dim>\r
+void ConsLaw<dim>::build_fe() {\r
+ fe_ptr = new FESystem<dim>(FE_Q<dim>(1), N_COMP);\r
+}\r
+\r
+ // Bye bye Conservation law.\r
+template <int dim>\r
+ConsLaw<dim>::~ConsLaw () \r
+{\r
+ dof_handler.clear ();\r
+ delete fe_ptr;\r
+}\r
+\r
+ // @sect3{Initialize System}\r
+ // Sizes all of the vectors and sets up the\r
+ // sparsity patter. This function is called at\r
+ // the very beginning of a simulation. The function\r
+ // <code> setup_system </code> repeats some of these\r
+ // chores and is called after adaptivity in leiu\r
+ // of this function.\r
+template <int dim>\r
+void ConsLaw<dim>::initialize_system ()\r
+{\r
+ // First we need to distribute the\r
+ // DoFs.\r
+ dof_handler.clear();\r
+ dof_handler.distribute_dofs (*fe_ptr);\r
+ \r
+ // Size all of the fields.\r
+ solution.reinit (dof_handler.n_dofs());\r
+ nlsolution.reinit (dof_handler.n_dofs());\r
+ predictor.reinit (dof_handler.n_dofs());\r
+ ppsolution.reinit (dof_handler.n_dofs());\r
+ dsolution.reinit (dof_handler.n_dofs());\r
+ right_hand_side.reinit (dof_handler.n_dofs());\r
+ indicator.reinit(triangulation.n_active_cells());\r
+}\r
+\r
+ // @sect3{Setup System}\r
+ // We call this function to build the sparsity\r
+ // and the matrix.\r
+template <int dim>\r
+void ConsLaw<dim>::setup_system ()\r
+{\r
+\r
+ // The DoFs of a cell are coupled\r
+ // with all DoFs of all neighboring\r
+ // cells. Therefore the maximum\r
+ // number of matrix entries per row\r
+ // is needed when all neighbors of\r
+ // a cell are once more refined\r
+ // than the cell under\r
+ // consideration.\r
+ sparsity_pattern.reinit (dof_handler.n_dofs(),\r
+ dof_handler.n_dofs(),\r
+ (GeometryInfo<dim>::faces_per_cell\r
+ *GeometryInfo<dim>::subfaces_per_face+1)*fe_ptr->dofs_per_cell);\r
+ \r
+ // Since the continuous sparsity pattern is\r
+ // a subset of the DG one, and since we need\r
+ // the DG terms for handling hanging nodes, we use\r
+ // the flux pattern.\r
+ DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);\r
+ \r
+ sparsity_pattern.compress();\r
+ \r
+ // Rebuild the map. In serial this doesn't do much,\r
+ // but is needed. In parallel, this would desribe\r
+ // the parallel dof layout.\r
+ if (Map) delete Map;\r
+ Map = new Epetra_Map(dof_handler.n_dofs(), 0, *Comm);\r
+\r
+ // Epetra can build a more efficient matrix if\r
+ // one knows ahead of time the maxiumum number of\r
+ // columns in any row entry. We traverse the sparsity\r
+ // to discover this.\r
+ int cur_row = 0;\r
+ int cur_col = 0;\r
+ int max_df = -1;\r
+ for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); \r
+ s_i != sparsity_pattern.end(); s_i++) {\r
+ if (s_i->row() != cur_row) {\r
+ cur_col = 0;\r
+ cur_row = s_i->row();\r
+ }\r
+ cur_col++;\r
+ if (cur_col >= max_df) max_df = cur_col;\r
+ }\r
+\r
+ if (cur_col >= max_df) max_df = cur_col;\r
+ std::cout << "max_df:" << max_df << std::endl;\r
+\r
+ // Now we build the matrix, using the constructor\r
+ // that optimizes with the <code> max_df </code> variable.\r
+ if (Matrix) delete Matrix;\r
+ Matrix = new Epetra_CrsMatrix(Copy, *Map, max_df+1, true);\r
+\r
+ // We add the sparsity pattern to the matrix by\r
+ // inserting zeros.\r
+ std::vector<double> vals(max_df, 0);\r
+ std::vector<int> row_indices(max_df);\r
+ \r
+ cur_row = 0;\r
+ cur_col = 0;\r
+ for (SparsityPattern::iterator s_i = sparsity_pattern.begin(); \r
+ s_i != sparsity_pattern.end(); s_i++) {\r
+ if (s_i->row() != cur_row) {\r
+ Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);\r
+ cur_col = 0;\r
+ cur_row = s_i->row();\r
+ }\r
+ row_indices[cur_col++] = s_i->column();\r
+ }\r
+ // The last row.\r
+ Matrix->InsertGlobalValues(cur_row, cur_col, &vals[0], &row_indices[0]);\r
+\r
+ // Epetra requires this function after building or\r
+ // filling a matrix. It typically does some parallel\r
+ // bookeeping; perhaps more.\r
+ Matrix->FillComplete();\r
+\r
+}\r
+\r
+ // @sect3{Solving the linear system}\r
+ // Actually solve the linear system, using either\r
+ // Aztec of Amesos.\r
+template <int dim>\r
+void ConsLaw<dim>::solve (Vector<double> &dsolution, int &niter, double &lin_residual) \r
+{\r
+\r
+ // We must hand the solvers Epetra vectors.\r
+ // Luckily, they support the concept of a \r
+ // 'view', so we just send in a pointer to our\r
+ // dealii vectors.\r
+ Epetra_Vector x(View, *Map, dsolution.begin());\r
+ Epetra_Vector b(View, *Map, right_hand_side.begin());\r
+\r
+ // The Direct option selects the Amesos solver.\r
+ if (solver_params.SOLVER == solver_params_type::DIRECT) {\r
+ \r
+ // Setup for solving with\r
+ // Amesos.\r
+ Epetra_LinearProblem prob;\r
+ prob.SetOperator(Matrix);\r
+ Amesos_BaseSolver *solver;\r
+ Amesos Factory;\r
+\r
+ // Other solvers are available\r
+ // and may be selected by changing this\r
+ // string.\r
+ char *stype = "Amesos_Klu";\r
+\r
+ solver = Factory.Create(stype, prob);\r
+\r
+ Assert (solver != NULL, ExcInternalError());\r
+\r
+ // There are two parts to the direct solve.\r
+ // As I understand, the symbolic part figures\r
+ // out the sparsity patterns, and then the\r
+ // numerical part actually performs Gaussian\r
+ // elimination or whatever the approach is.\r
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
+ std::cout << "Starting Symbolic fact\n" << std::flush;\r
+\r
+ solver->SymbolicFactorization();\r
+\r
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
+ std::cout << "Starting Numeric fact\n" << std::flush;\r
+\r
+ solver->NumericFactorization();\r
+\r
+ \r
+ // Define the linear problem by setting the\r
+ // right hand and left hand sides.\r
+ prob.SetRHS(&b);\r
+ prob.SetLHS(&x);\r
+ // And finally solve the problem.\r
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
+ std::cout << "Starting solve\n" << std::flush;\r
+ solver->Solve();\r
+ niter = 0;\r
+ lin_residual = 0;\r
+\r
+ // We must free the solver that was created\r
+ // for us.\r
+ delete solver;\r
+\r
+ } else if (solver_params.SOLVER == solver_params_type::GMRES) {\r
+\r
+ // For the iterative solvers, we use Aztec.\r
+ AztecOO Solver;\r
+\r
+ // Select the appropriate level of verbosity.\r
+ if (solver_params.OUTPUT == solver_params_type::QUIET)\r
+ Solver.SetAztecOption(AZ_output, AZ_none);\r
+\r
+ if (solver_params.OUTPUT == solver_params_type::VERBOSE)\r
+ Solver.SetAztecOption(AZ_output, AZ_all);\r
+\r
+ // Select gmres. Other solvers are available.\r
+ Solver.SetAztecOption(AZ_solver, AZ_gmres);\r
+ Solver.SetRHS(&b);\r
+ Solver.SetLHS(&x);\r
+\r
+ // Set up the ILUT preconditioner. I do not know\r
+ // why, but we must pretend like we are in parallel\r
+ // using domain decomposition or the preconditioner\r
+ // refuses to activate.\r
+ Solver.SetAztecOption(AZ_precond, AZ_dom_decomp);\r
+ Solver.SetAztecOption(AZ_subdomain_solve, AZ_ilut);\r
+ Solver.SetAztecOption(AZ_overlap, 0);\r
+ Solver.SetAztecOption(AZ_reorder, 0);\r
+\r
+ // ILUT parameters as described above.\r
+ Solver.SetAztecParam(AZ_drop, solver_params.ILUT_DROP);\r
+ Solver.SetAztecParam(AZ_ilut_fill, solver_params.ILUT_FILL);\r
+ Solver.SetAztecParam(AZ_athresh, solver_params.ILUT_ATOL);\r
+ Solver.SetAztecParam(AZ_rthresh, solver_params.ILUT_RTOL);\r
+ Solver.SetUserMatrix(Matrix);\r
+\r
+ // Run the solver iteration. Collect the number\r
+ // of iterations and the residual.\r
+ Solver.Iterate(solver_params.MAX_ITERS, solver_params.RES);\r
+ niter = Solver.NumIters();\r
+ lin_residual = Solver.TrueResidual();\r
+ }\r
+}\r
+\r
+ // @sect3{Postprocessing and Output}\r
+ // Recover the physical variables from the conservative\r
+ // variables so that output will be (perhaps) more\r
+ // meaningfull.\r
+template <int dim>\r
+void ConsLaw<dim>::postprocess() {\r
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
+ std::vector<unsigned int> dofs (dofs_per_cell);\r
+ UpdateFlags update_flags = update_values\r
+ | update_gradients\r
+ | update_q_points\r
+ | update_JxW_values;\r
+ UpdateFlags update_flags1 = update_values\r
+ | update_gradients\r
+ | update_q_points\r
+ | update_JxW_values;\r
+\r
+ QGauss<dim> quadrature_formula(4);\r
+\r
+ const std::vector<Point<dim> > &us = fe_ptr->base_element(0).get_unit_support_points();\r
+\r
+\r
+ Quadrature<dim> unit_support(us);\r
+\r
+ int n_q_points = quadrature_formula.n_quadrature_points;\r
+ int n_uq_points = unit_support.n_quadrature_points;\r
+\r
+ FEValues<dim> fe_v (\r
+ mapping, *fe_ptr, quadrature_formula, update_flags);\r
+\r
+ FEValues<dim> fe_v_unit (\r
+ mapping, *fe_ptr, unit_support, update_flags1);\r
+\r
+ std::vector<Vector<double> > U(n_uq_points,\r
+ Vector<double>(get_n_components()));\r
+ std::vector<Vector<double> > UU(n_q_points,\r
+ Vector<double>(get_n_components()));\r
+ std::vector<std::vector<Tensor<1,dim> > > dU(n_uq_points,\r
+ std::vector<Tensor<1,dim> >(get_n_components()));\r
+ \r
+ typename DoFHandler<dim>::active_cell_iterator\r
+ cell = dof_handler.begin_active(),\r
+ endc = dof_handler.end();\r
+\r
+ // Loop the cells\r
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
+ cell->get_dof_indices (dofs);\r
+ fe_v_unit.reinit(cell);\r
+ fe_v.reinit(cell);\r
+\r
+ fe_v_unit.get_function_values(solution, U);\r
+ fe_v_unit.get_function_grads(solution, dU);\r
+ fe_v.get_function_values(solution, UU);\r
+\r
+ const std::vector<double> &JxW = fe_v.get_JxW_values ();\r
+\r
+ for (int q = 0; q < fe_v.get_fe().base_element(0).n_dofs_per_cell(); q++) {\r
+ unsigned int didx = fe_v.get_fe().component_to_system_index(DENS_IDX, q);\r
+ unsigned int eidx = fe_v.get_fe().component_to_system_index(ENERGY_IDX, q);\r
+ double rho_normVsqr = 0;\r
+ for (int d = 0; d < dim; d++) {\r
+ unsigned int vidx = fe_v.get_fe().component_to_system_index(d, q);\r
+ ppsolution(dofs[vidx]) = solution(dofs[vidx])/solution(dofs[didx]);\r
+ rho_normVsqr += solution(dofs[vidx])*solution(dofs[vidx]);\r
+ }\r
+ rho_normVsqr /= solution(dofs[didx]);\r
+ // Pressure\r
+ ppsolution(dofs[eidx]) = (GAMMA-1.0)*(solution(dofs[eidx]) - 0.5*rho_normVsqr);\r
+\r
+ // Either output density or gradient squared of density,\r
+ // depending on what the user wants.\r
+ if (!schlieren_plot) {\r
+ ppsolution(dofs[didx]) = solution(dofs[didx]);\r
+ } else {\r
+ double ng = 0;\r
+ for (int i = 0; i < dim; i++) ng += dU[q][DENS_IDX][i]*dU[q][DENS_IDX][i];\r
+ ng = std::sqrt(ng);\r
+ ppsolution(dofs[didx]) = ng;\r
+ }\r
+ }\r
+\r
+ } // cell\r
+\r
+}\r
+\r
+ // Loop and assign a value for refinement. We\r
+ // simply use the density squared, which selects\r
+ // shocks with some success.\r
+template <int dim>\r
+void ConsLaw<dim>::estimate() {\r
+ \r
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;\r
+ std::vector<unsigned int> dofs (dofs_per_cell);\r
+ UpdateFlags update_flags = update_values\r
+ | update_gradients\r
+ | update_q_points\r
+ | update_JxW_values;\r
+\r
+ QGauss<dim> quadrature_formula(1);\r
+ int n_q_points = quadrature_formula.n_quadrature_points;\r
+\r
+\r
+ FEValues<dim> fe_v (\r
+ mapping, *fe_ptr, quadrature_formula, update_flags);\r
+\r
+ std::vector<Vector<double> > U(n_q_points,\r
+ Vector<double>(get_n_components()));\r
+ std::vector<std::vector<Tensor<1,dim> > > dU(n_q_points,\r
+ std::vector<Tensor<1,dim> >(get_n_components()));\r
+ \r
+ typename DoFHandler<dim>::active_cell_iterator\r
+ cell = dof_handler.begin_active(),\r
+ endc = dof_handler.end();\r
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
+ fe_v.reinit(cell);\r
+\r
+ fe_v.get_function_values(predictor, U);\r
+ fe_v.get_function_grads(predictor, dU);\r
+\r
+ indicator(cell_no) = 0;\r
+ for (int q = 0; q < n_q_points; q++) {\r
+ double ng = 0;\r
+ for (int d = 0; d < dim; d++) ng += dU[q][DENS_IDX][d]*dU[q][DENS_IDX][d];\r
+\r
+ indicator(cell_no) += std::log(1+std::sqrt(ng));\r
+ \r
+ } \r
+ indicator(cell_no) /= n_q_points;\r
+\r
+ } \r
+}\r
+\r
+template <int dim>\r
+void ConsLaw<dim>::refine_grid ()\r
+{\r
+\r
+ SolutionTransfer<dim, double> soltrans(dof_handler);\r
+\r
+ typename DoFHandler<dim>::active_cell_iterator\r
+ cell = dof_handler.begin_active(),\r
+ endc = dof_handler.end();\r
+\r
+ // Loop cells. If the indicator\r
+ // for the cell matches the refinement criterion,\r
+ // refine, else unrefine. The unrefinement has\r
+ // a slight hysterisis to avoid 'flashing' from refined\r
+ // to unrefined.\r
+ for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) {\r
+ cell->clear_coarsen_flag();\r
+ cell->clear_refine_flag();\r
+ if (cell->level() < refinement_params.shock_levels &&\r
+ std::fabs(indicator(cell_no)) > refinement_params.shock_val ) {\r
+ cell->set_refine_flag();\r
+ } else {\r
+ if (cell->level() > 0 &&\r
+ std::fabs(indicator(cell_no)) < 0.75*refinement_params.shock_val)\r
+ cell->set_coarsen_flag();\r
+ }\r
+ }\r
+\r
+ // The following code prolongs the solution\r
+ // to the new grid and carries out the refinement.\r
+ std::vector<Vector<double> > interp_in;\r
+ std::vector<Vector<double> > interp_out;\r
+\r
+ interp_in.push_back(solution);\r
+ interp_in.push_back(predictor);\r
+\r
+ triangulation.prepare_coarsening_and_refinement();\r
+ soltrans.prepare_for_coarsening_and_refinement(interp_in);\r
+\r
+ triangulation.execute_coarsening_and_refinement ();\r
+\r
+ dof_handler.clear();\r
+ dof_handler.distribute_dofs (*fe_ptr);\r
+\r
+ {\r
+ Vector<double> new_solution(1);\r
+ Vector<double> new_predictor(1);\r
+\r
+ interp_out.push_back(new_solution);\r
+ interp_out.push_back(new_predictor);\r
+ interp_out[0].reinit(dof_handler.n_dofs());\r
+ interp_out[1].reinit(dof_handler.n_dofs());\r
+ }\r
+\r
+ soltrans.interpolate(interp_in, interp_out);\r
+ \r
+ // Let the vector delete a very small vector\r
+ solution.reinit(1);\r
+ predictor.reinit(1);\r
+ solution.swap(interp_out[0]);\r
+ predictor.swap(interp_out[1]);\r
+\r
+ // resize these vectors for the new grid.\r
+ nlsolution.reinit(dof_handler.n_dofs());\r
+ ppsolution.reinit(dof_handler.n_dofs());\r
+ nlsolution = solution;\r
+ dsolution.reinit (dof_handler.n_dofs());\r
+ right_hand_side.reinit (dof_handler.n_dofs());\r
+\r
+ indicator.reinit(triangulation.n_active_cells());\r
+\r
+}\r
+\r
+template <int dim>\r
+void ConsLaw<dim>::output_results (const unsigned int cycle) const\r
+{\r
+ char filename[512];\r
+ std::sprintf(filename, "solution-%03d.vtk", cycle);\r
+ std::ofstream output (filename);\r
+\r
+ DataOut<dim> data_out;\r
+ data_out.attach_dof_handler (dof_handler);\r
+ std::vector<std::string> solution_names;\r
+\r
+ // Rename the output with the physical variable\r
+ // names. Send the post-processed values.\r
+ solution_names.clear();\r
+ for (int i = 0; i < dim; i++) {\r
+ char buf[512];\r
+ std::sprintf(buf, "v_%d", i);\r
+ solution_names.push_back (buf); \r
+ }\r
+ solution_names.push_back("density");\r
+ solution_names.push_back("pressure");\r
+ data_out.add_data_vector (ppsolution, solution_names);\r
+\r
+ data_out.add_data_vector (indicator, "error");\r
+ data_out.build_patches ();\r
+ data_out.write_vtk (output);\r
+\r
+ output.close();\r
+}\r
+\r
+ // @sect3{Parsing the Input Deck}\r
+ // Declare the parameters for the\r
+ // input deck. We assume a certain\r
+ // maximum number of boundaries and process\r
+ // any boundary the user supplies up to\r
+ // that maximum number. We\r
+ // leave a detailed explanation of these \r
+ // parameters to our description of the input\r
+ // sample file.\r
+const UInt MAX_BD = 10;\r
+template <int dim>\r
+void ConsLaw<dim>::declare_parameters() {\r
+\r
+ // Global scope parameters/\r
+ prm.declare_entry("mesh", "grid.inp",\r
+ Patterns::Anything(),\r
+ "intput file");\r
+\r
+ prm.declare_entry("diffusion power", "2.0",\r
+ Patterns::Double(),\r
+ "power of mesh size for diffusion");\r
+\r
+ prm.declare_entry("gravity", "0.0",\r
+ Patterns::Double(),\r
+ "gravity forcing");\r
+\r
+ // Time stepping block\r
+ prm.enter_subsection("time stepping");\r
+ prm.declare_entry("time step", "0.1",\r
+ Patterns::Double(),\r
+ "simulation time step");\r
+ prm.declare_entry("final time", "10.0",\r
+ Patterns::Double(),\r
+ "simulation end time");\r
+ prm.leave_subsection();\r
+\r
+\r
+ // Declare the boundary parameters\r
+ for (int b = 0; b < MAX_BD; b++) {\r
+ char bd[512];\r
+ std::sprintf(bd, "boundary_%d", b);\r
+ prm.enter_subsection(bd);\r
+ prm.declare_entry("no penetration", "false",\r
+ Patterns::Selection("true|false"),\r
+ "<true|false>");\r
+ // declare a slot for each of the conservative\r
+ // variables.\r
+ for (int di = 0; di < N_COMP; di++) {\r
+ char var[512];\r
+ std::sprintf(var, "w_%d", di);\r
+ prm.declare_entry(var, "outflow",\r
+ Patterns::Selection(\r
+ "inflow|outflow|pressure"),\r
+ "<inflow|outflow|pressure>");\r
+ \r
+ // for dirichlet, a function in x,y,z\r
+ std::sprintf(var, "w_%d value", di);\r
+ prm.declare_entry(var, "0.0",\r
+ Patterns::Anything(),\r
+ "expression in x,y,z");\r
+ }\r
+\r
+ prm.leave_subsection();\r
+ }\r
+\r
+ // Initial condition block.\r
+ prm.enter_subsection("initial condition");\r
+ for (int di = 0; di < N_COMP; di++) {\r
+ char var[512];\r
+ std::sprintf(var, "w_%d", di);\r
+ \r
+ // for dirichlet, a function in x,y,z\r
+ std::sprintf(var, "w_%d value", di);\r
+ prm.declare_entry(var, "0.0",\r
+ Patterns::Anything(),\r
+ "expression in x,y,z");\r
+ }\r
+ prm.leave_subsection();\r
+\r
+ // The linear solver block.\r
+ prm.enter_subsection("linear solver");\r
+ prm.declare_entry("output", "quiet",\r
+ Patterns::Selection(\r
+ "quiet|verbose"),\r
+ "<quiet|verbose>");\r
+ prm.declare_entry("method", "gmres",\r
+ Patterns::Selection(\r
+ "gmres|direct"),\r
+ "<gmres|direct>");\r
+ prm.declare_entry("residual", "1e-10",\r
+ Patterns::Double(),\r
+ "linear solver residual");\r
+ prm.declare_entry("max iters", "300",\r
+ Patterns::Double(),\r
+ "maximum solver iterations");\r
+ prm.declare_entry("ilut fill", "2",\r
+ Patterns::Double(),\r
+ "ilut preconditioner fill");\r
+ prm.declare_entry("ilut absolute tolerance", "1e-9",\r
+ Patterns::Double(),\r
+ "ilut preconditioner tolerance");\r
+ prm.declare_entry("ilut relative tolerance", "1.1",\r
+ Patterns::Double(),\r
+ "rel tol");\r
+ prm.declare_entry("ilut drop tolerance", "1e-10",\r
+ Patterns::Double(),\r
+ "ilut drop tol");\r
+ prm.leave_subsection();\r
+\r
+\r
+ // A refinement controller block.\r
+ prm.enter_subsection("refinement");\r
+ prm.declare_entry("refinement", "none",\r
+ Patterns::Selection(\r
+ "none|fixed number|shock"),\r
+ "<on|off>");\r
+ prm.declare_entry("refinement fraction", "0.1",\r
+ Patterns::Double(),\r
+ "Fraction of high refinement");\r
+ prm.declare_entry("unrefinement fraction", "0.1",\r
+ Patterns::Double(),\r
+ "Fraction of low unrefinement");\r
+ prm.declare_entry("max elements", "1000000",\r
+ Patterns::Double(),\r
+ "maximum number of elements");\r
+ prm.declare_entry("shock value", "4.0",\r
+ Patterns::Double(),\r
+ "value for shock indicator");\r
+ prm.declare_entry("shock levels", "3.0",\r
+ Patterns::Double(),\r
+ "number of shock refinement levels");\r
+ prm.leave_subsection();\r
+\r
+ // Output control.\r
+ prm.enter_subsection("output");\r
+ prm.declare_entry("density", "standard",\r
+ Patterns::Selection(\r
+ "standard|schlieren"),\r
+ "<standard|schlieren>");\r
+ prm.declare_entry("step", "-1",\r
+ Patterns::Double(),\r
+ "output once per this period");\r
+ prm.leave_subsection();\r
+\r
+ // Flux control\r
+ prm.enter_subsection("flux");\r
+ prm.declare_entry("stab", "alpha",\r
+ Patterns::Selection(\r
+ "alpha|constant|mesh"),\r
+ "<alpha|constant|mesh>");\r
+ prm.declare_entry("stab value", "1",\r
+ Patterns::Double(),\r
+ "alpha stabilization");\r
+ prm.leave_subsection();\r
+\r
+\r
+}\r
+\r
+ // Code to actually parse an input file. This function\r
+ // matches the declarations above.\r
+template <int dim>\r
+void ConsLaw<dim>::load_parameters(const char *infile){\r
+\r
+ prm.read_input(infile);\r
+\r
+ // The global parameters.\r
+ mesh = prm.get("mesh");\r
+\r
+ diffusion_power = prm.get_double("diffusion power");\r
+\r
+ gravity = prm.get_double("gravity");\r
+\r
+ // The time stepping.\r
+ prm.enter_subsection("time stepping");\r
+ dT = prm.get_double("time step");\r
+ std::cout << "dT=" << dT << std::endl;\r
+ if (dT == 0) {\r
+ is_stationary = true;\r
+ dT = 1.0;\r
+ TF = 1.0;\r
+ std::cout << "Stationary mode" << std::endl;\r
+ }\r
+ TF = prm.get_double("final time");\r
+ std::cout << "TF=" << TF << std::endl;\r
+ prm.leave_subsection();\r
+\r
+ // The boundary info\r
+ for (int b = 0; b < MAX_BD; b++) {\r
+ std::vector<bc_type> flags(N_COMP, OUTFLOW_BC);\r
+\r
+ // Define a parser for every boundary, though it may be\r
+ // unused.\r
+ SideCondition<dim> *sd = new SideCondition<dim>(N_COMP);\r
+ char bd[512];\r
+ std::sprintf(bd, "boundary_%d", b);\r
+ prm.enter_subsection(bd);\r
+\r
+ const std::string &nopen = prm.get("no penetration");\r
+\r
+ // Determine how each component is handled.\r
+ for (int di = 0; di < N_COMP; di++) {\r
+ char var[512];\r
+ std::sprintf(var, "w_%d", di);\r
+ std::string btype = prm.get(var);\r
+ std::sprintf(var, "w_%d value", di);\r
+ std::string var_value = prm.get(var);\r
+\r
+ if (di < dim && nopen == "true") {\r
+ flags[di] = NO_PENETRATION_BC;\r
+ } else if (btype == "inflow") {\r
+ flags[di] = INFLOW_BC;\r
+ sd->set_coeff_row(di, var_value); \r
+ } else if (btype == "pressure") {\r
+ flags[di] = PRESSURE_BC;\r
+ sd->set_coeff_row(di, var_value); \r
+ }\r
+ } \r
+ prm.leave_subsection();\r
+\r
+ // Add the boundary condition to the law.\r
+ sd->Init();\r
+ add_boundary(b, flags, sd);\r
+ }\r
+\r
+ // Initial conditions.\r
+ prm.enter_subsection("initial condition");\r
+ for (int di = 0; di < N_COMP; di++) {\r
+ char var[512];\r
+\r
+ std::sprintf(var, "w_%d value", di);\r
+ std::string var_value = prm.get(var);\r
+ ic.set_ic(di, var_value); \r
+ }\r
+ ic.Init();\r
+ prm.leave_subsection();\r
+\r
+ // The linear solver.\r
+ prm.enter_subsection("linear solver");\r
+ const std::string &op = prm.get("output");\r
+ if (op == "verbose") solver_params.OUTPUT = solver_params_type::VERBOSE;\r
+ if (op == "quiet") solver_params.OUTPUT = solver_params_type::QUIET;\r
+ const std::string &sv = prm.get("method");\r
+ if (sv == "direct") {\r
+ solver_params.SOLVER = solver_params_type::DIRECT;\r
+ } else if (sv == "gmres") {\r
+ solver_params.SOLVER = solver_params_type::GMRES;\r
+ } \r
+\r
+ solver_params.RES = prm.get_double("residual");\r
+ solver_params.MAX_ITERS = (int) prm.get_double("max iters");\r
+ solver_params.ILUT_FILL = prm.get_double("ilut fill");\r
+ solver_params.ILUT_ATOL = prm.get_double("ilut absolute tolerance");\r
+ solver_params.ILUT_RTOL = prm.get_double("ilut relative tolerance");\r
+ solver_params.ILUT_DROP = prm.get_double("ilut drop tolerance");\r
+ solver_params.RES = prm.get_double("residual");\r
+ prm.leave_subsection();\r
+\r
+\r
+ // And refiement.\r
+ prm.enter_subsection("refinement");\r
+ const std::string &ref = prm.get("refinement");\r
+ if (ref == "none") {\r
+ refinement_params.refine = refinement_params_type::NONE;\r
+ } else if (ref == "fixed number") {\r
+ refinement_params.refine = refinement_params_type::FIXED_NUMBER;\r
+ } else if (ref == "shock") {\r
+ refinement_params.refine = refinement_params_type::SHOCK;\r
+ } else\r
+ refinement_params.high_frac = prm.get_double("refinement fraction");\r
+ refinement_params.high_frac_sav = refinement_params.high_frac;\r
+ refinement_params.low_frac = prm.get_double("unrefinement fraction");\r
+ refinement_params.max_cells = prm.get_double("max elements");\r
+ refinement_params.shock_val = prm.get_double("shock value");\r
+ refinement_params.shock_levels = prm.get_double("shock levels");\r
+ prm.leave_subsection();\r
+ \r
+ // Output control.\r
+ prm.enter_subsection("output");\r
+ const std::string &dens = prm.get("density");\r
+ schlieren_plot = dens == "schlieren" ? true : false;\r
+ output_step = prm.get_double("step");\r
+ prm.leave_subsection();\r
+\r
+ // Flux control.\r
+ prm.enter_subsection("flux");\r
+ const std::string &stab = prm.get("stab");\r
+ if (stab == "constant") {\r
+ flux_params.LF_stab = flux_params_type::CONSTANT;\r
+ } else if (stab == "mesh ") {\r
+ flux_params.LF_stab = flux_params_type::MESH;\r
+ }\r
+ flux_params.LF_stab_value = prm.get_double("stab value");\r
+ prm.leave_subsection();\r
+\r
+\r
+}\r
+\r
+template<int dim>\r
+void ConsLaw<dim>::zero_matrix() {\r
+ Matrix->PutScalar(0); Matrix->FillComplete();\r
+}\r
+\r
+ // We use a predictor to try and make adaptivity\r
+ // work better. The idea is to try and refine ahead\r
+ // of a front, rather than stepping into a coarse\r
+ // set of elements and smearing the solution. This\r
+ // simple time extrapolator does the job.\r
+template<int dim>\r
+void ConsLaw<dim>::compute_predictor() {\r
+ predictor = nlsolution;\r
+ predictor.sadd(3/2.0, -1/2.0, solution);\r
+}\r
+\r
+ // @sect3{Run the simulation} Contains the initialization,\r
+ // the time loop, and the inner Newton iteration.\r
+template <int dim>\r
+void ConsLaw<dim>::run () \r
+{\r
+\r
+ // Open and load the mesh.\r
+ GridIn<dim> grid_in;\r
+ grid_in.attach_triangulation(triangulation);\r
+ std::cout << "Opening mesh <" << mesh << ">" << std::endl;\r
+ std::ifstream input_file(mesh.c_str(), std::ios::in);\r
+\r
+ Assert (infile,\r
+ ExcFileNotOpen());\r
+\r
+ grid_in.read_ucd(input_file); \r
+ input_file.close();\r
+ \r
+ build_fe();\r
+\r
+ unsigned int nstep = 0;\r
+ \r
+ // Initialize fields and matrices.\r
+ initialize_system (); \r
+ setup_system();\r
+ initialize(); \r
+ predictor = solution;\r
+\r
+ // Initial refinement. We apply the ic,\r
+ // estimate, refine, and repeat until\r
+ // happy.\r
+ if (refinement_params.refine != refinement_params_type::NONE)\r
+ for (int i = 0; i < refinement_params.shock_levels; i++) {\r
+ estimate();\r
+ refine_grid();\r
+ setup_system();\r
+ initialize(); \r
+ predictor = solution;\r
+ }\r
+ postprocess();\r
+ output_results (nstep);\r
+\r
+ // Determine when we will output next.\r
+ double next_output = T + output_step;\r
+\r
+ // @sect4{Main time stepping loop}\r
+ predictor = solution;\r
+ while(T < TF)\r
+ {\r
+ std::cout << "T=" << T << ", ";\r
+\r
+\r
+ std::cout << " Number of active cells: "\r
+ << triangulation.n_active_cells()\r
+ << std::endl;\r
+\r
+\r
+ std::cout << " Number of degrees of freedom: "\r
+ << dof_handler.n_dofs()\r
+ << std::endl;\r
+\r
+ bool nonlin_done = false;\r
+ double res_norm;\r
+ int lin_iter;\r
+\r
+ // Print some relevant information during the\r
+ // Newton iteration.\r
+ std::cout << "NonLin Res: Lin Iter Lin Res" << std::endl;\r
+ std::cout << "______________________________________" << std::endl;\r
+\r
+ int max_nonlin = 7;\r
+ int nonlin_iter = 0;\r
+ double lin_res;\r
+\r
+ // @sect5{Newton iteration}\r
+ nlsolution = predictor;\r
+ while (!nonlin_done) {\r
+ lin_iter = 0;\r
+ zero_matrix();\r
+ right_hand_side = 0;\r
+ assemble_system (res_norm);\r
+ // Flash a star to the screen so one can\r
+ // know when the assembly has stopped and the linear\r
+ // solution is starting.\r
+ std::cout << "* " << std::flush;\r
+\r
+ // Test against a (hardcoded) nonlinear tolderance.\r
+ // Do not solve the linear system at the last step \r
+ // (since it would be a waste).\r
+ \r
+ if (fabs(res_norm) < 1e-10) {\r
+ nonlin_done = true;\r
+ } else {\r
+ // Solve the linear system and update with the\r
+ // delta.\r
+ dsolution = 0;\r
+ solve (dsolution, lin_iter, lin_res);\r
+ nlsolution.add(1.0, dsolution);\r
+ }\r
+\r
+ // Print the residuals.\r
+ std::printf("%-16.3e %04d %-5.2e\n",\r
+ res_norm, lin_iter, lin_res);\r
+\r
+ nonlin_iter++;\r
+ } \r
+\r
+ // Various post convergence tasks.\r
+ compute_predictor();\r
+\r
+ solution = nlsolution;\r
+\r
+\r
+ estimate();\r
+\r
+ postprocess();\r
+\r
+ T += dT;\r
+\r
+ // Output if it is time.\r
+ if (output_step < 0) {\r
+ output_results (++nstep);\r
+ } else if (T >= next_output) {\r
+ output_results (++nstep);\r
+ next_output += output_step;\r
+ }\r
+\r
+ // Refine, if refinement is selected.\r
+ if (refinement_params.refine != refinement_params_type::NONE) {\r
+ refine_grid();\r
+ setup_system();\r
+ }\r
+ }\r
+}\r
+\r
+ // The following ``main'' function is\r
+ // similar to previous examples and\r
+ // need not to be commented on.\r
+int main (int argc, char *argv[]) \r
+{\r
+\r
+ MPI_Init(&argc, &argv);\r
+ Comm = new Epetra_MpiComm(MPI_COMM_WORLD);\r
+\r
+ if (argc != 2) {\r
+ std::cout << "Usage:" << argv[0] << " infile" << std::endl;\r
+ std::exit(1);\r
+ }\r
+ try\r
+ {\r
+ ConsLaw<DIMENSION> cons;\r
+ cons.declare_parameters();\r
+ cons.load_parameters(argv[1]);\r
+ cons.run ();\r
+ }\r
+ catch (std::exception &exc)\r
+ {\r
+ std::cerr << std::endl << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ std::cerr << "Exception on processing: " << std::endl\r
+ << exc.what() << std::endl\r
+ << "Aborting!" << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ return 1;\r
+ }\r
+ catch (...) \r
+ {\r
+ std::cerr << std::endl << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ std::cerr << "Unknown exception!" << std::endl\r
+ << "Aborting!" << std::endl\r
+ << "----------------------------------------------------"\r
+ << std::endl;\r
+ return 1;\r
+ };\r
+ \r
+ return 0;\r
+}\r
+\r