Changeset 4729
 Timestamp:
 Sep 12, 2007, 10:06:01 AM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4727 r4729 60 60 // dq1=q(vertex1)q(vertex0) 61 61 // dq2=q(vertex2)q(vertex0) 62 62 63 if (dq0>=0.0){ 63 64 if (dq1>=dq2){ … … 114 115 // calculate a multiplicative factor phi by which the provisional 115 116 // vertex jumps are to be limited 117 116 118 int i; 117 119 double r=1000.0, r0=1.0, phi=1.0; … … 121 123 // Any provisional jump with magnitude < TINY does not contribute to 122 124 // the limiting process. 123 125 for (i=0;i<3;i++){ 124 126 if (dqv[i]<TINY) 125 127 r0=qmin/dqv[i]; 128 126 129 if (dqv[i]>TINY) 127 130 r0=qmax/dqv[i]; 131 128 132 r=min(r0,r); 129 133 } … … 132 136 for (i=0;i<3;i++) 133 137 dqv[i]=dqv[i]*phi; 138 134 139 return 0; 135 140 } … … 830 835 831 836 832 // FIXME (Ole): Split this function up into gateway and computational function 833 // as done with the others (most recently flux 11/9/7). 834 // That will make the code faster and more readable. 835 // 837 // Computational routine 838 int _extrapolate_second_order_sw(int number_of_elements, 839 double epsilon, 840 double minimum_allowed_height, 841 double beta_w, 842 double beta_w_dry, 843 double beta_uh, 844 double beta_uh_dry, 845 double beta_vh, 846 double beta_vh_dry, 847 long* surrogate_neighbours, 848 long* number_of_boundaries, 849 double* centroid_coordinates, 850 double* stage_centroid_values, 851 double* xmom_centroid_values, 852 double* ymom_centroid_values, 853 double* elevation_centroid_values, 854 double* vertex_coordinates, 855 double* stage_vertex_values, 856 double* xmom_vertex_values, 857 double* ymom_vertex_values, 858 double* elevation_vertex_values, 859 int optimise_dry_cells) { 860 861 862 863 // Local variables 864 double a, b; // Gradient vector used to calculate vertex values from centroids 865 int k,k0,k1,k2,k3,k6,coord_index,i; 866 double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle 867 double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2; 868 double dqv[3], qmin, qmax, hmin, hmax; 869 double hc, h0, h1, h2, beta_tmp; 870 871 872 for (k=0; k<number_of_elements; k++) { 873 k3=k*3; 874 k6=k*6; 875 876 877 if (number_of_boundaries[k]==3){ 878 // No neighbours, set gradient on the triangle to zero 879 880 stage_vertex_values[k3] = stage_centroid_values[k]; 881 stage_vertex_values[k3+1] = stage_centroid_values[k]; 882 stage_vertex_values[k3+2] = stage_centroid_values[k]; 883 xmom_vertex_values[k3] = xmom_centroid_values[k]; 884 xmom_vertex_values[k3+1] = xmom_centroid_values[k]; 885 xmom_vertex_values[k3+2] = xmom_centroid_values[k]; 886 ymom_vertex_values[k3] = ymom_centroid_values[k]; 887 ymom_vertex_values[k3+1] = ymom_centroid_values[k]; 888 ymom_vertex_values[k3+2] = ymom_centroid_values[k]; 889 890 continue; 891 } 892 else { 893 // Triangle k has one or more neighbours. 894 // Get centroid and vertex coordinates of the triangle 895 896 // Get the vertex coordinates 897 xv0 = vertex_coordinates[k6]; yv0=vertex_coordinates[k6+1]; 898 xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3]; 899 xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5]; 900 901 // Get the centroid coordinates 902 coord_index=2*k; 903 x=centroid_coordinates[coord_index]; 904 y=centroid_coordinates[coord_index+1]; 905 906 // Store x and y differentials for the vertices of 907 // triangle k relative to the centroid 908 dxv0=xv0x; dxv1=xv1x; dxv2=xv2x; 909 dyv0=yv0y; dyv1=yv1y; dyv2=yv2y; 910 } 911 912 913 914 if (number_of_boundaries[k]<=1){ 915 916 //============================================== 917 // Number of boundaries <= 1 918 //============================================== 919 920 921 // If no boundaries, auxiliary triangle is formed 922 // from the centroids of the three neighbours 923 // If one boundary, auxiliary triangle is formed 924 // from this centroid and its two neighbours 925 926 k0=surrogate_neighbours[k3]; 927 k1=surrogate_neighbours[k3+1]; 928 k2=surrogate_neighbours[k3+2]; 929 930 // Get the auxiliary triangle's vertex coordinates 931 // (really the centroids of neighbouring triangles) 932 coord_index=2*k0; 933 x0=centroid_coordinates[coord_index]; 934 y0=centroid_coordinates[coord_index+1]; 935 936 coord_index=2*k1; 937 x1=centroid_coordinates[coord_index]; 938 y1=centroid_coordinates[coord_index+1]; 939 940 coord_index=2*k2; 941 x2=centroid_coordinates[coord_index]; 942 y2=centroid_coordinates[coord_index+1]; 943 944 // Store x and y differentials for the vertices 945 // of the auxiliary triangle 946 dx1=x1x0; dx2=x2x0; 947 dy1=y1y0; dy2=y2y0; 948 949 // Calculate 2*area of the auxiliary triangle 950 // The triangle is guaranteed to be counterclockwise 951 area2 = dy2*dx1  dy1*dx2; 952 953 // If the mesh is 'weird' near the boundary, 954 // the triangle might be flat or clockwise: 955 if (area2<=0) { 956 PyErr_SetString(PyExc_RuntimeError, 957 "shallow_water_ext.c: negative triangle area encountered"); 958 return 1; 959 } 960 961 // Calculate heights of neighbouring cells 962 hc = stage_centroid_values[k]  elevation_centroid_values[k]; 963 h0 = stage_centroid_values[k0]  elevation_centroid_values[k0]; 964 h1 = stage_centroid_values[k1]  elevation_centroid_values[k1]; 965 h2 = stage_centroid_values[k2]  elevation_centroid_values[k2]; 966 hmin = min(h0,min(h1,h2)); 967 968 969 if (optimise_dry_cells) { 970 // Check if linear reconstruction is necessary for triangle k 971 // This check will exclude dry cells. 972 973 hmax = max(h0,max(h1,h2)); 974 if (hmax < epsilon) { 975 continue; 976 } 977 } 978 979 980 // 981 // stage 982 // 983 984 // Calculate the difference between vertex 0 of the auxiliary 985 // triangle and the centroid of triangle k 986 dq0=stage_centroid_values[k0]stage_centroid_values[k]; 987 988 // Calculate differentials between the vertices 989 // of the auxiliary triangle (centroids of neighbouring triangles) 990 dq1=stage_centroid_values[k1]stage_centroid_values[k0]; 991 dq2=stage_centroid_values[k2]stage_centroid_values[k0]; 992 993 // Calculate the gradient of stage on the auxiliary triangle 994 a = dy2*dq1  dy1*dq2; 995 a /= area2; 996 b = dx1*dq2  dx2*dq1; 997 b /= area2; 998 999 // Calculate provisional jumps in stage from the centroid 1000 // of triangle k to its vertices, to be limited 1001 dqv[0]=a*dxv0+b*dyv0; 1002 dqv[1]=a*dxv1+b*dyv1; 1003 dqv[2]=a*dxv2+b*dyv2; 1004 1005 // Now we want to find min and max of the centroid and the 1006 // vertices of the auxiliary triangle and compute jumps 1007 // from the centroid to the min and max 1008 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 1009 1010 // Playing with dry wet interface 1011 hmin = qmin; 1012 beta_tmp = beta_w; 1013 if (hmin<minimum_allowed_height) 1014 beta_tmp = beta_w_dry; 1015 1016 // Limit the gradient 1017 limit_gradient(dqv,qmin,qmax,beta_tmp); 1018 1019 for (i=0;i<3;i++) 1020 stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i]; 1021 1022 1023 // 1024 // xmomentum 1025 // 1026 1027 // Calculate the difference between vertex 0 of the auxiliary 1028 // triangle and the centroid of triangle k 1029 dq0=xmom_centroid_values[k0]xmom_centroid_values[k]; 1030 1031 // Calculate differentials between the vertices 1032 // of the auxiliary triangle 1033 dq1=xmom_centroid_values[k1]xmom_centroid_values[k0]; 1034 dq2=xmom_centroid_values[k2]xmom_centroid_values[k0]; 1035 1036 // Calculate the gradient of xmom on the auxiliary triangle 1037 a = dy2*dq1  dy1*dq2; 1038 a /= area2; 1039 b = dx1*dq2  dx2*dq1; 1040 b /= area2; 1041 1042 // Calculate provisional jumps in stage from the centroid 1043 // of triangle k to its vertices, to be limited 1044 dqv[0]=a*dxv0+b*dyv0; 1045 dqv[1]=a*dxv1+b*dyv1; 1046 dqv[2]=a*dxv2+b*dyv2; 1047 1048 // Now we want to find min and max of the centroid and the 1049 // vertices of the auxiliary triangle and compute jumps 1050 // from the centroid to the min and max 1051 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 1052 beta_tmp = beta_uh; 1053 if (hmin<minimum_allowed_height) 1054 beta_tmp = beta_uh_dry; 1055 1056 // Limit the gradient 1057 limit_gradient(dqv,qmin,qmax,beta_tmp); 1058 1059 for (i=0;i<3;i++) 1060 xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i]; 1061 1062 1063 // 1064 // ymomentum 1065 // 1066 1067 // Calculate the difference between vertex 0 of the auxiliary 1068 // triangle and the centroid of triangle k 1069 dq0=ymom_centroid_values[k0]ymom_centroid_values[k]; 1070 1071 // Calculate differentials between the vertices 1072 // of the auxiliary triangle 1073 dq1=ymom_centroid_values[k1]ymom_centroid_values[k0]; 1074 dq2=ymom_centroid_values[k2]ymom_centroid_values[k0]; 1075 1076 // Calculate the gradient of xmom on the auxiliary triangle 1077 a = dy2*dq1  dy1*dq2; 1078 a /= area2; 1079 b = dx1*dq2  dx2*dq1; 1080 b /= area2; 1081 1082 // Calculate provisional jumps in stage from the centroid 1083 // of triangle k to its vertices, to be limited 1084 dqv[0]=a*dxv0+b*dyv0; 1085 dqv[1]=a*dxv1+b*dyv1; 1086 dqv[2]=a*dxv2+b*dyv2; 1087 1088 // Now we want to find min and max of the centroid and the 1089 // vertices of the auxiliary triangle and compute jumps 1090 // from the centroid to the min and max 1091 find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); 1092 beta_tmp = beta_vh; 1093 1094 if (hmin<minimum_allowed_height) 1095 beta_tmp = beta_vh_dry; 1096 1097 // Limit the gradient 1098 limit_gradient(dqv,qmin,qmax,beta_tmp); 1099 1100 for (i=0;i<3;i++) 1101 ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i]; 1102 1103 } // End number_of_boundaries <=1 1104 else{ 1105 1106 //============================================== 1107 // Number of boundaries == 2 1108 //============================================== 1109 1110 // One internal neighbour and gradient is in direction of the neighbour's centroid 1111 1112 // Find the only internal neighbour 1113 for (k2=k3;k2<k3+3;k2++){ 1114 // Find internal neighbour of triangle k 1115 // k2 indexes the edges of triangle k 1116 1117 if (surrogate_neighbours[k2]!=k) 1118 break; 1119 } 1120 if ((k2==k3+3)) { 1121 // If we didn't find an internal neighbour 1122 PyErr_SetString(PyExc_RuntimeError, 1123 "shallow_water_ext.c: Internal neighbour not found"); 1124 return 1; 1125 } 1126 1127 k1=surrogate_neighbours[k2]; 1128 1129 // The coordinates of the triangle are already (x,y). 1130 // Get centroid of the neighbour (x1,y1) 1131 coord_index=2*k1; 1132 x1=centroid_coordinates[coord_index]; 1133 y1=centroid_coordinates[coord_index+1]; 1134 1135 // Compute x and y distances between the centroid of 1136 // triangle k and that of its neighbour 1137 dx1=x1x; dy1=y1y; 1138 1139 // Set area2 as the square of the distance 1140 area2=dx1*dx1+dy1*dy1; 1141 1142 // Set dx2=(x1x0)/((x1x0)^2+(y1y0)^2) 1143 // and dy2=(y1y0)/((x1x0)^2+(y1y0)^2) which 1144 // respectively correspond to the x and y gradients 1145 // of the conserved quantities 1146 dx2=1.0/area2; 1147 dy2=dx2*dy1; 1148 dx2*=dx1; 1149 1150 1151 1152 // 1153 // stage 1154 // 1155 1156 // Compute differentials 1157 dq1=stage_centroid_values[k1]stage_centroid_values[k]; 1158 1159 // Calculate the gradient between the centroid of triangle k 1160 // and that of its neighbour 1161 a=dq1*dx2; 1162 b=dq1*dy2; 1163 1164 // Calculate provisional vertex jumps, to be limited 1165 dqv[0]=a*dxv0+b*dyv0; 1166 dqv[1]=a*dxv1+b*dyv1; 1167 dqv[2]=a*dxv2+b*dyv2; 1168 1169 // Now limit the jumps 1170 if (dq1>=0.0){ 1171 qmin=0.0; 1172 qmax=dq1; 1173 } 1174 else{ 1175 qmin=dq1; 1176 qmax=0.0; 1177 } 1178 1179 // Limit the gradient 1180 limit_gradient(dqv,qmin,qmax,beta_w); 1181 1182 for (i=0;i<3;i++) 1183 stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i]; 1184 1185 1186 // 1187 // xmomentum 1188 // 1189 1190 // Compute differentials 1191 dq1=xmom_centroid_values[k1]xmom_centroid_values[k]; 1192 1193 // Calculate the gradient between the centroid of triangle k 1194 // and that of its neighbour 1195 a=dq1*dx2; 1196 b=dq1*dy2; 1197 1198 // Calculate provisional vertex jumps, to be limited 1199 dqv[0]=a*dxv0+b*dyv0; 1200 dqv[1]=a*dxv1+b*dyv1; 1201 dqv[2]=a*dxv2+b*dyv2; 1202 1203 // Now limit the jumps 1204 if (dq1>=0.0){ 1205 qmin=0.0; 1206 qmax=dq1; 1207 } 1208 else{ 1209 qmin=dq1; 1210 qmax=0.0; 1211 } 1212 1213 // Limit the gradient 1214 limit_gradient(dqv,qmin,qmax,beta_w); 1215 1216 for (i=0;i<3;i++) 1217 xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i]; 1218 1219 1220 // 1221 // ymomentum 1222 // 1223 1224 // Compute differentials 1225 dq1=ymom_centroid_values[k1]ymom_centroid_values[k]; 1226 1227 // Calculate the gradient between the centroid of triangle k 1228 // and that of its neighbour 1229 a=dq1*dx2; 1230 b=dq1*dy2; 1231 1232 // Calculate provisional vertex jumps, to be limited 1233 dqv[0]=a*dxv0+b*dyv0; 1234 dqv[1]=a*dxv1+b*dyv1; 1235 dqv[2]=a*dxv2+b*dyv2; 1236 1237 // Now limit the jumps 1238 if (dq1>=0.0){ 1239 qmin=0.0; 1240 qmax=dq1; 1241 } 1242 else{ 1243 qmin=dq1; 1244 qmax=0.0; 1245 } 1246 1247 // Limit the gradient 1248 limit_gradient(dqv,qmin,qmax,beta_w); 1249 1250 for (i=0;i<3;i++) 1251 ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i]; 1252 1253 } // else [number_of_boundaries==2] 1254 } // for k=0 to number_of_elements1 1255 1256 return 0; 1257 } 1258 1259 836 1260 PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) { 1261 /*Compute the vertex values based on a linear reconstruction on each triangle 1262 These values are calculated as follows: 1263 1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle 1264 formed by the centroids of its three neighbours. 1265 2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product 1266 of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average 1267 of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the 1268 original triangle has gradient (a,b). 1269 3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions 1270 find_qmin_and_qmax and limit_gradient 1271 1272 Python call: 1273 extrapolate_second_order_sw(domain.surrogate_neighbours, 1274 domain.number_of_boundaries 1275 domain.centroid_coordinates, 1276 Stage.centroid_values 1277 Xmom.centroid_values 1278 Ymom.centroid_values 1279 domain.vertex_coordinates, 1280 Stage.vertex_values, 1281 Xmom.vertex_values, 1282 Ymom.vertex_values) 1283 1284 Post conditions: 1285 The vertices of each triangle have values from a limited linear reconstruction 1286 based on centroid values 1287 1288 */ 1289 PyArrayObject *surrogate_neighbours, 1290 *number_of_boundaries, 1291 *centroid_coordinates, 1292 *stage_centroid_values, 1293 *xmom_centroid_values, 1294 *ymom_centroid_values, 1295 *elevation_centroid_values, 1296 *vertex_coordinates, 1297 *stage_vertex_values, 1298 *xmom_vertex_values, 1299 *ymom_vertex_values, 1300 *elevation_vertex_values; 1301 PyObject *domain, *Tmp; 1302 1303 1304 double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry; 1305 double minimum_allowed_height, epsilon; 1306 int optimise_dry_cells, number_of_elements, e; 1307 1308 // Provisional jumps from centroids to v'tices and safety factor re limiting 1309 // by which these jumps are limited 1310 // Convert Python arguments to C 1311 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi", 1312 &domain, 1313 &surrogate_neighbours, 1314 &number_of_boundaries, 1315 ¢roid_coordinates, 1316 &stage_centroid_values, 1317 &xmom_centroid_values, 1318 &ymom_centroid_values, 1319 &elevation_centroid_values, 1320 &vertex_coordinates, 1321 &stage_vertex_values, 1322 &xmom_vertex_values, 1323 &ymom_vertex_values, 1324 &elevation_vertex_values, 1325 &optimise_dry_cells)) { 1326 1327 PyErr_SetString(PyExc_RuntimeError, "Input arguments to extrapolate_second_order_sw failed"); 1328 return NULL; 1329 } 1330 1331 // FIXME (Ole): Investigate if it is quicker to obtain all input arguments using GetAttrString rather than ParseTuple. 1332 // It certainly looked as if passing domain.epsilon is slowed things down 1333 1334 // Get the safety factor beta_w, set in the config.py file. This is used in the limiting process 1335 Tmp = PyObject_GetAttrString(domain, "beta_w"); 1336 if (!Tmp) { 1337 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain"); 1338 return NULL; 1339 } 1340 beta_w = PyFloat_AsDouble(Tmp); 1341 Py_DECREF(Tmp); 1342 1343 Tmp = PyObject_GetAttrString(domain, "beta_w_dry"); 1344 if (!Tmp) { 1345 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain"); 1346 return NULL; 1347 } 1348 beta_w_dry = PyFloat_AsDouble(Tmp); 1349 Py_DECREF(Tmp); 1350 1351 Tmp = PyObject_GetAttrString(domain, "beta_uh"); 1352 if (!Tmp) { 1353 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain"); 1354 return NULL; 1355 } 1356 beta_uh = PyFloat_AsDouble(Tmp); 1357 Py_DECREF(Tmp); 1358 1359 Tmp = PyObject_GetAttrString(domain, "beta_uh_dry"); 1360 if (!Tmp) { 1361 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain"); 1362 return NULL; 1363 } 1364 beta_uh_dry = PyFloat_AsDouble(Tmp); 1365 Py_DECREF(Tmp); 1366 1367 Tmp = PyObject_GetAttrString(domain, "beta_vh"); 1368 if (!Tmp) { 1369 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain"); 1370 return NULL; 1371 } 1372 beta_vh = PyFloat_AsDouble(Tmp); 1373 Py_DECREF(Tmp); 1374 1375 Tmp = PyObject_GetAttrString(domain, "beta_vh_dry"); 1376 if (!Tmp) { 1377 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain"); 1378 return NULL; 1379 } 1380 beta_vh_dry = PyFloat_AsDouble(Tmp); 1381 Py_DECREF(Tmp); 1382 1383 Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height"); 1384 if (!Tmp) { 1385 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height"); 1386 return NULL; 1387 } 1388 minimum_allowed_height = PyFloat_AsDouble(Tmp); 1389 Py_DECREF(Tmp); 1390 1391 Tmp = PyObject_GetAttrString(domain, "epsilon"); 1392 if (!Tmp) { 1393 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon"); 1394 return NULL; 1395 } 1396 epsilon = PyFloat_AsDouble(Tmp); 1397 Py_DECREF(Tmp); 1398 1399 // Call underlying computational routine 1400 number_of_elements = stage_centroid_values > dimensions[0]; 1401 e = _extrapolate_second_order_sw(number_of_elements, 1402 epsilon, 1403 minimum_allowed_height, 1404 beta_w, 1405 beta_w_dry, 1406 beta_uh, 1407 beta_uh_dry, 1408 beta_vh, 1409 beta_vh_dry, 1410 (long*) surrogate_neighbours > data, 1411 (long*) number_of_boundaries > data, 1412 (double*) centroid_coordinates > data, 1413 (double*) stage_centroid_values > data, 1414 (double*) xmom_centroid_values > data, 1415 (double*) ymom_centroid_values > data, 1416 (double*) elevation_centroid_values > data, 1417 (double*) vertex_coordinates > data, 1418 (double*) stage_vertex_values > data, 1419 (double*) xmom_vertex_values > data, 1420 (double*) ymom_vertex_values > data, 1421 (double*) elevation_vertex_values > data, 1422 optimise_dry_cells); 1423 if (e == 1) { 1424 // Use error string set inside computational routine 1425 return NULL; 1426 } 1427 1428 1429 return Py_BuildValue(""); 1430 1431 }// extrapolate_secondorder_sw 1432 1433 1434 1435 // FIXME (Ole): This function is obsolete as of 12 Sep 2007 1436 PyObject *extrapolate_second_order_sw_original(PyObject *self, PyObject *args) { 837 1437 /*Compute the vertex values based on a linear reconstruction on each triangle 838 1438 These values are calculated as follows: … … 887 1487 double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp; 888 1488 double minimum_allowed_height; 1489 889 1490 // Provisional jumps from centroids to v'tices and safety factor re limiting 890 1491 // by which these jumps are limited
Note: See TracChangeset
for help on using the changeset viewer.