Changeset 6840


Ignore:
Timestamp:
Apr 20, 2009, 2:31:32 PM (15 years ago)
Author:
nariman
Message:

Optimisation of _extrapolate_second_order_sw and limit_gradient.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r6778 r6840  
    133133 
    134134  phi=min(r*beta_w,1.0);
    135   for (i=0;i<3;i++)
    136     dqv[i]=dqv[i]*phi;
     135  //for (i=0;i<3;i++)
     136  dqv[0]=dqv[0]*phi;
     137  dqv[1]=dqv[1]*phi;
     138  dqv[2]=dqv[2]*phi;
    137139   
    138140  return 0;
     
    500502        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    501503        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    502         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     504        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
    503505        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    504506
     
    10621064// Computational routine
    10631065int _extrapolate_second_order_sw(int number_of_elements,
    1064                   double epsilon,
    1065                   double minimum_allowed_height,
    1066                   double beta_w,
    1067                   double beta_w_dry,
    1068                   double beta_uh,
    1069                   double beta_uh_dry,
    1070                   double beta_vh,
    1071                   double beta_vh_dry,
    1072                   long* surrogate_neighbours,
    1073                   long* number_of_boundaries,
    1074                   double* centroid_coordinates,
    1075                   double* stage_centroid_values,
    1076                   double* xmom_centroid_values,
    1077                   double* ymom_centroid_values,
    1078                   double* elevation_centroid_values,
    1079                   double* vertex_coordinates,
    1080                   double* stage_vertex_values,
    1081                   double* xmom_vertex_values,
    1082                   double* ymom_vertex_values,
    1083                   double* elevation_vertex_values,
    1084                   int optimise_dry_cells) {
     1066                                 double epsilon,
     1067                                 double minimum_allowed_height,
     1068                                 double beta_w,
     1069                                 double beta_w_dry,
     1070                                 double beta_uh,
     1071                                 double beta_uh_dry,
     1072                                 double beta_vh,
     1073                                 double beta_vh_dry,
     1074                                 long* surrogate_neighbours,
     1075                                 long* number_of_boundaries,
     1076                                 double* centroid_coordinates,
     1077                                 double* stage_centroid_values,
     1078                                 double* xmom_centroid_values,
     1079                                 double* ymom_centroid_values,
     1080                                 double* elevation_centroid_values,
     1081                                 double* vertex_coordinates,
     1082                                 double* stage_vertex_values,
     1083                                 double* xmom_vertex_values,
     1084                                 double* ymom_vertex_values,
     1085                                 double* elevation_vertex_values,
     1086                                 int optimise_dry_cells) {
    10851087                 
    10861088                 
     
    10881090  // Local variables
    10891091  double a, b; // Gradient vector used to calculate vertex values from centroids
    1090   int k,k0,k1,k2,k3,k6,coord_index,i;
    1091   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
    1092   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
     1092  int k, k0, k1, k2, k3, k6, coord_index, i;
     1093  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
     1094  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
    10931095  double dqv[3], qmin, qmax, hmin, hmax;
    10941096  double hc, h0, h1, h2, beta_tmp, hfactor;
    10951097
    10961098   
    1097   for (k=0; k<number_of_elements; k++) {
     1099  for (k = 0; k < number_of_elements; k++)
     1100  {
    10981101    k3=k*3;
    10991102    k6=k*6;
    11001103
    1101    
    1102     if (number_of_boundaries[k]==3){
     1104    if (number_of_boundaries[k]==3)
     1105    {
    11031106      // No neighbours, set gradient on the triangle to zero
    11041107     
     
    11151118      continue;
    11161119    }
    1117     else {
     1120    else
     1121    {
    11181122      // Triangle k has one or more neighbours.
    11191123      // Get centroid and vertex coordinates of the triangle
    11201124     
    11211125      // Get the vertex coordinates
    1122       xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
    1123       xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
    1124       xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
     1126      xv0 = vertex_coordinates[k6];   
     1127      yv0 = vertex_coordinates[k6+1];
     1128      xv1 = vertex_coordinates[k6+2];
     1129      yv1 = vertex_coordinates[k6+3];
     1130      xv2 = vertex_coordinates[k6+4];
     1131      yv2 = vertex_coordinates[k6+5];
    11251132     
    11261133      // Get the centroid coordinates
    1127       coord_index=2*k;
    1128       x=centroid_coordinates[coord_index];
    1129       y=centroid_coordinates[coord_index+1];
     1134      coord_index = 2*k;
     1135      x = centroid_coordinates[coord_index];
     1136      y = centroid_coordinates[coord_index+1];
    11301137     
    11311138      // Store x- and y- differentials for the vertices of
    11321139      // triangle k relative to the centroid
    1133       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    1134       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
     1140      dxv0 = xv0 - x;
     1141      dxv1 = xv1 - x;
     1142      dxv2 = xv2 - x;
     1143      dyv0 = yv0 - y;
     1144      dyv1 = yv1 - y;
     1145      dyv2 = yv2 - y;
    11351146    }
    11361147
    11371148
    11381149           
    1139     if (number_of_boundaries[k]<=1){
    1140    
     1150    if (number_of_boundaries[k]<=1)
     1151    {
    11411152      //==============================================
    11421153      // Number of boundaries <= 1
     
    11491160      // from this centroid and its two neighbours
    11501161     
    1151       k0=surrogate_neighbours[k3];
    1152       k1=surrogate_neighbours[k3+1];
    1153       k2=surrogate_neighbours[k3+2];
     1162      k0 = surrogate_neighbours[k3];
     1163      k1 = surrogate_neighbours[k3 + 1];
     1164      k2 = surrogate_neighbours[k3 + 2];
    11541165     
    11551166      // Get the auxiliary triangle's vertex coordinates
    11561167      // (really the centroids of neighbouring triangles)
    1157       coord_index=2*k0;
    1158       x0=centroid_coordinates[coord_index];
    1159       y0=centroid_coordinates[coord_index+1];
    1160      
    1161       coord_index=2*k1;
    1162       x1=centroid_coordinates[coord_index];
    1163       y1=centroid_coordinates[coord_index+1];
    1164      
    1165       coord_index=2*k2;
    1166       x2=centroid_coordinates[coord_index];
    1167       y2=centroid_coordinates[coord_index+1];
     1168      coord_index = 2*k0;
     1169      x0 = centroid_coordinates[coord_index];
     1170      y0 = centroid_coordinates[coord_index+1];
     1171     
     1172      coord_index = 2*k1;
     1173      x1 = centroid_coordinates[coord_index];
     1174      y1 = centroid_coordinates[coord_index+1];
     1175     
     1176      coord_index = 2*k2;
     1177      x2 = centroid_coordinates[coord_index];
     1178      y2 = centroid_coordinates[coord_index+1];
    11681179     
    11691180      // Store x- and y- differentials for the vertices
    11701181      // of the auxiliary triangle
    1171       dx1=x1-x0; dx2=x2-x0;
    1172       dy1=y1-y0; dy2=y2-y0;
     1182      dx1 = x1 - x0;
     1183      dx2 = x2 - x0;
     1184      dy1 = y1 - y0;
     1185      dy2 = y2 - y0;
    11731186     
    11741187      // Calculate 2*area of the auxiliary triangle
     
    11781191      // If the mesh is 'weird' near the boundary,
    11791192      // the triangle might be flat or clockwise:
    1180       if (area2<=0) {
    1181     PyErr_SetString(PyExc_RuntimeError,
    1182             "shallow_water_ext.c: negative triangle area encountered");
    1183     return -1;
     1193      if (area2 <= 0)
     1194      {
     1195        PyErr_SetString(PyExc_RuntimeError,
     1196                        "shallow_water_ext.c: negative triangle area encountered");
     1197   
     1198        return -1;
    11841199      } 
    11851200     
     
    11891204      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    11901205      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
    1191       hmin = min(min(h0,min(h1,h2)),hc);
     1206      hmin = min(min(h0, min(h1, h2)), hc);
    11921207      //hfactor = hc/(hc + 1.0);
    11931208
    11941209      hfactor = 0.0;
    1195       if (hmin > 0.001 ) {
    1196       hfactor = (hmin-0.001)/(hmin+0.004);
    1197       }
    1198      
    1199       if (optimise_dry_cells) {     
    1200     // Check if linear reconstruction is necessary for triangle k
    1201     // This check will exclude dry cells.
    1202 
    1203     hmax = max(h0,max(h1,h2));     
    1204     if (hmax < epsilon) {
    1205       continue;
    1206     }
    1207       }
    1208 
    1209            
     1210      if (hmin > 0.001)
     1211      {
     1212        hfactor = (hmin - 0.001)/(hmin + 0.004);
     1213      }
     1214     
     1215      if (optimise_dry_cells)
     1216      {     
     1217        // Check if linear reconstruction is necessary for triangle k
     1218        // This check will exclude dry cells.
     1219
     1220        hmax = max(h0, max(h1, h2));     
     1221        if (hmax < epsilon)
     1222        {
     1223            continue;
     1224        }
     1225      }
     1226   
    12101227      //-----------------------------------
    12111228      // stage
     
    12141231      // Calculate the difference between vertex 0 of the auxiliary
    12151232      // triangle and the centroid of triangle k
    1216       dq0=stage_centroid_values[k0]-stage_centroid_values[k];
     1233      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
    12171234     
    12181235      // Calculate differentials between the vertices
    12191236      // of the auxiliary triangle (centroids of neighbouring triangles)
    1220       dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
    1221       dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
    1222      
     1237      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
     1238      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
     1239     
     1240          inv_area2 = 1.0/area2;
    12231241      // Calculate the gradient of stage on the auxiliary triangle
    12241242      a = dy2*dq1 - dy1*dq2;
    1225       a /= area2;
     1243      a *= inv_area2;
    12261244      b = dx1*dq2 - dx2*dq1;
    1227       b /= area2;
     1245      b *= inv_area2;
    12281246     
    12291247      // Calculate provisional jumps in stage from the centroid
    12301248      // of triangle k to its vertices, to be limited
    1231       dqv[0]=a*dxv0+b*dyv0;
    1232       dqv[1]=a*dxv1+b*dyv1;
    1233       dqv[2]=a*dxv2+b*dyv2;
     1249      dqv[0] = a*dxv0 + b*dyv0;
     1250      dqv[1] = a*dxv1 + b*dyv1;
     1251      dqv[2] = a*dxv2 + b*dyv2;
    12341252     
    12351253      // Now we want to find min and max of the centroid and the
    12361254      // vertices of the auxiliary triangle and compute jumps
    12371255      // from the centroid to the min and max
    1238       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1256      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12391257     
    12401258      // Playing with dry wet interface
     
    12491267      //printf("beta_tmp = %f\n",beta_tmp);
    12501268      // Limit the gradient
    1251       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1252      
    1253       for (i=0;i<3;i++)
    1254     stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1269      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1270     
     1271      //for (i=0;i<3;i++)
     1272      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
     1273          stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
     1274          stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
    12551275     
    12561276     
     
    12611281      // Calculate the difference between vertex 0 of the auxiliary
    12621282      // triangle and the centroid of triangle k     
    1263       dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
     1283      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
    12641284     
    12651285      // Calculate differentials between the vertices
    12661286      // of the auxiliary triangle
    1267       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
    1268       dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
     1287      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
     1288      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
    12691289     
    12701290      // Calculate the gradient of xmom on the auxiliary triangle
    12711291      a = dy2*dq1 - dy1*dq2;
    1272       a /= area2;
     1292      a *= inv_area2;
    12731293      b = dx1*dq2 - dx2*dq1;
    1274       b /= area2;
     1294      b *= inv_area2;
    12751295     
    12761296      // Calculate provisional jumps in stage from the centroid
    12771297      // of triangle k to its vertices, to be limited     
    1278       dqv[0]=a*dxv0+b*dyv0;
    1279       dqv[1]=a*dxv1+b*dyv1;
    1280       dqv[2]=a*dxv2+b*dyv2;
     1298      dqv[0] = a*dxv0+b*dyv0;
     1299      dqv[1] = a*dxv1+b*dyv1;
     1300      dqv[2] = a*dxv2+b*dyv2;
    12811301     
    12821302      // Now we want to find min and max of the centroid and the
    12831303      // vertices of the auxiliary triangle and compute jumps
    12841304      // from the centroid to the min and max
    1285       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1305      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12861306      //beta_tmp = beta_uh;
    12871307      //if (hmin<minimum_allowed_height)
     
    12901310
    12911311      // Limit the gradient
    1292       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1293 
    1294       for (i=0;i<3;i++)
    1295     xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1296      
     1312      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1313
     1314      for (i=0; i < 3; i++)
     1315      {
     1316        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
     1317      }
    12971318     
    12981319      //-----------------------------------
     
    13021323      // Calculate the difference between vertex 0 of the auxiliary
    13031324      // triangle and the centroid of triangle k
    1304       dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
     1325      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
    13051326     
    13061327      // Calculate differentials between the vertices
    13071328      // of the auxiliary triangle
    1308       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
    1309       dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
     1329      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
     1330      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
    13101331     
    13111332      // Calculate the gradient of xmom on the auxiliary triangle
    13121333      a = dy2*dq1 - dy1*dq2;
    1313       a /= area2;
     1334      a *= inv_area2;
    13141335      b = dx1*dq2 - dx2*dq1;
    1315       b /= area2;
     1336      b *= inv_area2;
    13161337     
    13171338      // Calculate provisional jumps in stage from the centroid
    13181339      // of triangle k to its vertices, to be limited
    1319       dqv[0]=a*dxv0+b*dyv0;
    1320       dqv[1]=a*dxv1+b*dyv1;
    1321       dqv[2]=a*dxv2+b*dyv2;
     1340      dqv[0] = a*dxv0 + b*dyv0;
     1341      dqv[1] = a*dxv1 + b*dyv1;
     1342      dqv[2] = a*dxv2 + b*dyv2;
    13221343     
    13231344      // Now we want to find min and max of the centroid and the
    13241345      // vertices of the auxiliary triangle and compute jumps
    13251346      // from the centroid to the min and max
    1326       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1347      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    13271348     
    13281349      //beta_tmp = beta_vh;
     
    13331354
    13341355      // Limit the gradient
    1335       limit_gradient(dqv,qmin,qmax,beta_tmp);
     1356      limit_gradient(dqv, qmin, qmax, beta_tmp);
    13361357     
    13371358      for (i=0;i<3;i++)
    1338     ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1339    
     1359      {
     1360        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1361      }
    13401362    } // End number_of_boundaries <=1
    1341     else{
     1363    else
     1364    {
    13421365
    13431366      //==============================================
     
    13481371     
    13491372      // Find the only internal neighbour (k1?)
    1350       for (k2=k3;k2<k3+3;k2++){
    1351     // Find internal neighbour of triangle k     
    1352     // k2 indexes the edges of triangle k   
    1353      
    1354     if (surrogate_neighbours[k2]!=k)
    1355       break;
    1356       }
    1357      
    1358       if ((k2==k3+3)) {
    1359     // If we didn't find an internal neighbour
    1360     PyErr_SetString(PyExc_RuntimeError,
    1361             "shallow_water_ext.c: Internal neighbour not found");     
    1362     return -1;
    1363       }
    1364      
    1365       k1=surrogate_neighbours[k2];
     1373      for (k2 = k3; k2 < k3 + 3; k2++)
     1374      {
     1375      // Find internal neighbour of triangle k     
     1376      // k2 indexes the edges of triangle k   
     1377     
     1378          if (surrogate_neighbours[k2] != k)
     1379          {
     1380             break;
     1381          }
     1382      }
     1383     
     1384      if ((k2 == k3 + 3))
     1385      {
     1386        // If we didn't find an internal neighbour
     1387        PyErr_SetString(PyExc_RuntimeError,
     1388                        "shallow_water_ext.c: Internal neighbour not found");     
     1389        return -1;
     1390      }
     1391     
     1392      k1 = surrogate_neighbours[k2];
    13661393     
    13671394      // The coordinates of the triangle are already (x,y).
    13681395      // Get centroid of the neighbour (x1,y1)
    1369       coord_index=2*k1;
    1370       x1=centroid_coordinates[coord_index];
    1371       y1=centroid_coordinates[coord_index+1];
     1396      coord_index = 2*k1;
     1397      x1 = centroid_coordinates[coord_index];
     1398      y1 = centroid_coordinates[coord_index + 1];
    13721399     
    13731400      // Compute x- and y- distances between the centroid of
    13741401      // triangle k and that of its neighbour
    1375       dx1=x1-x; dy1=y1-y;
     1402      dx1 = x1 - x;
     1403      dy1 = y1 - y;
    13761404     
    13771405      // Set area2 as the square of the distance
    1378       area2=dx1*dx1+dy1*dy1;
     1406      area2 = dx1*dx1 + dy1*dy1;
    13791407     
    13801408      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
     
    13821410      // respectively correspond to the x- and y- gradients
    13831411      // of the conserved quantities
    1384       dx2=1.0/area2;
    1385       dy2=dx2*dy1;
    1386       dx2*=dx1;
     1412      dx2 = 1.0/area2;
     1413      dy2 = dx2*dy1;
     1414      dx2 *= dx1;
    13871415     
    13881416     
     
    13921420
    13931421      // Compute differentials
    1394       dq1=stage_centroid_values[k1]-stage_centroid_values[k];
     1422      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
    13951423     
    13961424      // Calculate the gradient between the centroid of triangle k
    13971425      // and that of its neighbour
    1398       a=dq1*dx2;
    1399       b=dq1*dy2;
     1426      a = dq1*dx2;
     1427      b = dq1*dy2;
    14001428     
    14011429      // Calculate provisional vertex jumps, to be limited
    1402       dqv[0]=a*dxv0+b*dyv0;
    1403       dqv[1]=a*dxv1+b*dyv1;
    1404       dqv[2]=a*dxv2+b*dyv2;
     1430      dqv[0] = a*dxv0 + b*dyv0;
     1431      dqv[1] = a*dxv1 + b*dyv1;
     1432      dqv[2] = a*dxv2 + b*dyv2;
    14051433     
    14061434      // Now limit the jumps
    1407       if (dq1>=0.0){
    1408     qmin=0.0;
    1409     qmax=dq1;
    1410       }
    1411       else{
    1412     qmin=dq1;
    1413     qmax=0.0;
     1435      if (dq1>=0.0)
     1436      {
     1437        qmin=0.0;
     1438        qmax=dq1;
     1439      }
     1440      else
     1441      {
     1442        qmin = dq1;
     1443        qmax = 0.0;
    14141444      }
    14151445     
    14161446      // Limit the gradient
    1417       limit_gradient(dqv,qmin,qmax,beta_w);
    1418      
    1419       for (i=0;i<3;i++)
    1420     stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
    1421      
     1447      limit_gradient(dqv, qmin, qmax, beta_w);
     1448     
     1449      //for (i=0; i < 3; i++)
     1450      //{
     1451      stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
     1452      stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
     1453      stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
     1454      //}
    14221455     
    14231456      //-----------------------------------
     
    14261459     
    14271460      // Compute differentials
    1428       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
     1461      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
    14291462     
    14301463      // Calculate the gradient between the centroid of triangle k
    14311464      // and that of its neighbour
    1432       a=dq1*dx2;
    1433       b=dq1*dy2;
     1465      a = dq1*dx2;
     1466      b = dq1*dy2;
    14341467     
    14351468      // Calculate provisional vertex jumps, to be limited
    1436       dqv[0]=a*dxv0+b*dyv0;
    1437       dqv[1]=a*dxv1+b*dyv1;
    1438       dqv[2]=a*dxv2+b*dyv2;
     1469      dqv[0] = a*dxv0+b*dyv0;
     1470      dqv[1] = a*dxv1+b*dyv1;
     1471      dqv[2] = a*dxv2+b*dyv2;
    14391472     
    14401473      // Now limit the jumps
    1441       if (dq1>=0.0){
    1442     qmin=0.0;
    1443     qmax=dq1;
    1444       }
    1445       else{
    1446     qmin=dq1;
    1447     qmax=0.0;
     1474      if (dq1 >= 0.0)
     1475      {
     1476        qmin = 0.0;
     1477        qmax = dq1;
     1478      }
     1479      else
     1480      {
     1481        qmin = dq1;
     1482        qmax = 0.0;
    14481483      }
    14491484     
    14501485      // Limit the gradient     
    1451       limit_gradient(dqv,qmin,qmax,beta_w);
    1452      
    1453       for (i=0;i<3;i++)
    1454     xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1455      
     1486      limit_gradient(dqv, qmin, qmax, beta_w);
     1487     
     1488      //for (i=0;i<3;i++)
     1489      //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
     1490      //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
     1491      //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
     1492     
     1493      for (i = 0; i < 3;i++)
     1494      {
     1495          xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
     1496      }
    14561497     
    14571498      //-----------------------------------
     
    14601501
    14611502      // Compute differentials
    1462       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
     1503      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
    14631504     
    14641505      // Calculate the gradient between the centroid of triangle k
    14651506      // and that of its neighbour
    1466       a=dq1*dx2;
    1467       b=dq1*dy2;
     1507      a = dq1*dx2;
     1508      b = dq1*dy2;
    14681509     
    14691510      // Calculate provisional vertex jumps, to be limited
    1470       dqv[0]=a*dxv0+b*dyv0;
    1471       dqv[1]=a*dxv1+b*dyv1;
    1472       dqv[2]=a*dxv2+b*dyv2;
     1511      dqv[0] = a*dxv0 + b*dyv0;
     1512      dqv[1] = a*dxv1 + b*dyv1;
     1513      dqv[2] = a*dxv2 + b*dyv2;
    14731514     
    14741515      // Now limit the jumps
    1475       if (dq1>=0.0){
    1476     qmin=0.0;
    1477     qmax=dq1;
    1478       }
    1479       else{
    1480     qmin=dq1;
    1481     qmax=0.0;
     1516      if (dq1>=0.0)
     1517      {
     1518        qmin = 0.0;
     1519        qmax = dq1;
     1520      }
     1521      else
     1522      {
     1523        qmin = dq1;
     1524        qmax = 0.0;
    14821525      }
    14831526     
    14841527      // Limit the gradient
    1485       limit_gradient(dqv,qmin,qmax,beta_w);
    1486      
    1487       for (i=0;i<3;i++)
    1488     ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1489    
     1528      limit_gradient(dqv, qmin, qmax, beta_w);
     1529     
     1530      //for (i=0;i<3;i++)
     1531      //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1532      //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1533      //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
     1534
     1535for (i=0;i<3;i++)
     1536        {
     1537ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1538        }
     1539//ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1540//ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1541//ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
    14901542    } // else [number_of_boundaries==2]
    14911543  } // for k=0 to number_of_elements-1
Note: See TracChangeset for help on using the changeset viewer.