Changeset 9015 for trunk

Timestamp:

Nov 3, 2013, 8:49:42 PM (11 years ago)

Author:

davies

Message:

Various changes including more aggressive velocity limiting
Seems to remove wb artefacts

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/bal_and/swb2_domain_ext.c (modified) (3 diffs)
• tests/dam_break/plotme.py (modified) (1 diff)

trunk/anuga_work/development/gareth/experimental/bal_and/swb2_domain_ext.c

@@ -958 +958 @@
 
   
-  
-  //
-  //for(k=0; k<number_of_elements;k++){ 
-  //    // Count the wet neighbours
-  //    count_wet_neighbours[k]=0;
-  //    for (i=0; i<3; i++){
-  //      ktmp = surrogate_neighbours[3*k+i];              
-  //      if( (ktmp >= 0) && cell_wetness_scale[ktmp]>0.){
-  //          count_wet_neighbours[k]+=1;
-  //      }else if(ktmp<=0){
-  //          // Boundary condition -- FIXME: assume wet
-  //          count_wet_neighbours[k]+=1;
-  //      }
-  //       
-  //      
-  //    }
-  //}
-
   if(extrapolate_velocity_second_order==1){
 
@@ -997 +979 @@
       }
 
-  // First treat all 'fully wet' cells, then treat 'partially dry' cells
-  //for(k_wetdry=0; k_wetdry<2; k_wetdry++){
-
-      // Begin extrapolation routine
-      for (k = 0; k < number_of_elements; k++) 
+  // Begin extrapolation routine
+  for (k = 0; k < number_of_elements; k++) 
+  {
+
+    // Useful indices
+    k3=k*3;
+    k6=k*6;
+
+    if (number_of_boundaries[k]==3)
+    {
+      // No neighbours, set gradient on the triangle to zero
+      
+      stage_edge_values[k3]   = stage_centroid_values[k];
+      stage_edge_values[k3+1] = stage_centroid_values[k];
+      stage_edge_values[k3+2] = stage_centroid_values[k];
+
+      //xmom_centroid_values[k] = 0.;
+      //ymom_centroid_values[k] = 0.;
+      
+      xmom_edge_values[k3]    = xmom_centroid_values[k];
+      xmom_edge_values[k3+1]  = xmom_centroid_values[k];
+      xmom_edge_values[k3+2]  = xmom_centroid_values[k];
+      ymom_edge_values[k3]    = ymom_centroid_values[k];
+      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
+      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
+      dk = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
+      height_edge_values[k3] = dk;
+      height_edge_values[k3+1] = dk;
+      height_edge_values[k3+2] = dk;
+      
+      continue;
+    }
+    else 
+    {
+      // Triangle k has one or more neighbours. 
+      // Get centroid and edge coordinates of the triangle
+      
+      // Get the edge coordinates
+      xv0 = edge_coordinates[k6];   
+      yv0 = edge_coordinates[k6+1];
+      xv1 = edge_coordinates[k6+2]; 
+      yv1 = edge_coordinates[k6+3];
+      xv2 = edge_coordinates[k6+4]; 
+      yv2 = edge_coordinates[k6+5];
+      
+      // Get the centroid coordinates
+      coord_index = 2*k;
+      x = centroid_coordinates[coord_index];
+      y = centroid_coordinates[coord_index+1];
+      
+      // Store x- and y- differentials for the edges of 
+      // triangle k relative to the centroid
+      dxv0 = xv0 - x; 
+      dxv1 = xv1 - x; 
+      dxv2 = xv2 - x;
+      dyv0 = yv0 - y; 
+      dyv1 = yv1 - y; 
+      dyv2 = yv2 - y;
+
+    }
+
+
+            
+    if (number_of_boundaries[k]<=1)
+    {
+      //==============================================
+      // Number of boundaries <= 1
+      // 'Typical case'
+      //==============================================    
+    
+    
+      // If no boundaries, auxiliary triangle is formed 
+      // from the centroids of the three neighbours
+      // If one boundary, auxiliary triangle is formed 
+      // from this centroid and its two neighbours
+      
+      k0 = surrogate_neighbours[k3];
+      k1 = surrogate_neighbours[k3 + 1];
+      k2 = surrogate_neighbours[k3 + 2];
+     
+      // Test to see whether we accept the surrogate neighbours 
+      // Note that if ki is replaced with k in more than 1 neighbour, then the
+      // triangle area will be zero, and a first order extrapolation will be
+      // used
+      // FIXME: Remove cell_wetness_scale if you don't need it
+      if(k2<0 || (cell_wetness_scale[k2]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          k2 = k ;
+      }
+      if(k0 < 0 || (cell_wetness_scale[k0]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          k0 = k ;
+      }
+      if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
+          k1 = k ;
+      }
+
+      // Take note if the max neighbour bed elevation is greater than the min
+      // neighbour stage -- suggests a 'steep' bed relative to the flow
+      //bedmax = max(elevation_centroid_values[k], 
+      //             max(elevation_centroid_values[k0],
+      //                 max(elevation_centroid_values[k1], 
+      //                     elevation_centroid_values[k2])));
+      //bedmin = min(elevation_centroid_values[k], 
+      //             min(elevation_centroid_values[k0],
+      //                 min(elevation_centroid_values[k1], 
+      //                     elevation_centroid_values[k2])));
+      //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
+      //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
+      //                   min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
+      //                       max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
+
+      // Get the auxiliary triangle's vertex coordinates 
+      // (really the centroids of neighbouring triangles)
+      coord_index = 2*k0;
+      x0 = centroid_coordinates[coord_index];
+      y0 = centroid_coordinates[coord_index+1];
+      
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index+1];
+      
+      coord_index = 2*k2;
+      x2 = centroid_coordinates[coord_index];
+      y2 = centroid_coordinates[coord_index+1];
+
+      // Store x- and y- differentials for the vertices 
+      // of the auxiliary triangle
+      dx1 = x1 - x0; 
+      dx2 = x2 - x0;
+      dy1 = y1 - y0; 
+      dy2 = y2 - y0;
+      
+      // Calculate 2*area of the auxiliary triangle
+      // The triangle is guaranteed to be counter-clockwise      
+      area2 = dy2*dx1 - dy1*dx2; 
+       
+      //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
+      if ((area2 <= 0.) )//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
       {
 
-        //// First treat all 'fully wet' cells, then treat 'partially dry' cells
-        //// For partially wet cells, we now know that the edges of neighbouring
-        //// fully wet cells have been defined
-        //if( cell_wetness_scale[k]==1.0*(1-k_wetdry) ){
-        //  continue;
-        //}
-
-        // Useful indices
-        k3=k*3;
-        k6=k*6;
-
-        if (number_of_boundaries[k]==3)
-        {
-          // No neighbours, set gradient on the triangle to zero
+          // Isolated wet cell -- constant stage/depth extrapolation 
+          //stage_edge_values[k3]   = stage_centroid_values[k];
+          //stage_edge_values[k3+1] = stage_centroid_values[k];
+          //stage_edge_values[k3+2] = stage_centroid_values[k];
+
+          dk=height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
+          height_edge_values[k3] = dk;
+          height_edge_values[k3+1] = dk;
+          height_edge_values[k3+2] = dk;
           
-          stage_edge_values[k3]   = stage_centroid_values[k];
-          stage_edge_values[k3+1] = stage_centroid_values[k];
-          stage_edge_values[k3+2] = stage_centroid_values[k];
-    
-          //xmom_centroid_values[k] = 0.;
-          //ymom_centroid_values[k] = 0.;
-          
+          stage_edge_values[k3]   = elevation_centroid_values[k]+dk;
+          stage_edge_values[k3+1] = elevation_centroid_values[k]+dk;
+          stage_edge_values[k3+2] = elevation_centroid_values[k]+dk;
+          //stage_edge_values[k3] = elevation_edge_values[k3]+dk;
+          //stage_edge_values[k3+1] = elevation_edge_values[k3+1]+dk;
+          //stage_edge_values[k3+2] = elevation_edge_values[k3+2]+dk;
+
           xmom_edge_values[k3]    = xmom_centroid_values[k];
           xmom_edge_values[k3+1]  = xmom_centroid_values[k];
@@ -1032 +1141 @@
           ymom_edge_values[k3+1]  = ymom_centroid_values[k];
           ymom_edge_values[k3+2]  = ymom_centroid_values[k];
-          dk = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
-          height_edge_values[k3] = dk;
-          height_edge_values[k3+1] = dk;
-          height_edge_values[k3+2] = dk;
-          
+
           continue;
-        }
-        else 
-        {
-          // Triangle k has one or more neighbours. 
-          // Get centroid and edge coordinates of the triangle
-          
-          // Get the edge coordinates
-          xv0 = edge_coordinates[k6];   
-          yv0 = edge_coordinates[k6+1];
-          xv1 = edge_coordinates[k6+2]; 
-          yv1 = edge_coordinates[k6+3];
-          xv2 = edge_coordinates[k6+4]; 
-          yv2 = edge_coordinates[k6+5];
-          
-          // Get the centroid coordinates
-          coord_index = 2*k;
-          x = centroid_coordinates[coord_index];
-          y = centroid_coordinates[coord_index+1];
-          
-          // Store x- and y- differentials for the edges of 
-          // triangle k relative to the centroid
-          dxv0 = xv0 - x; 
-          dxv1 = xv1 - x; 
-          dxv2 = xv2 - x;
-          dyv0 = yv0 - y; 
-          dyv1 = yv1 - y; 
-          dyv2 = yv2 - y;
-
-
-          // Compute bed slope as a reference
-          //area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0);
-          //inv_area_e=1.0/area_e;
-          //a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) -
-          //       (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
-          //a_bs *= inv_area_e;
-
-          //b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) +
-          //       (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);
-          //b_bs *= inv_area_e;
-
-          //printf("slopes: %f, %f \n", a_bs, b_bs);
-        }
-
-
-                
-        if (number_of_boundaries[k]<=1)
-        {
-          //==============================================
-          // Number of boundaries <= 1
-          // 'Typical case'
-          //==============================================    
+      }  
+      
+      // Calculate heights of neighbouring cells
+      hc = stage_centroid_values[k]  - elevation_centroid_values[k];
+      h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+      
+      hmin = min(min(h0, min(h1, h2)), hc);
+      hmax = max(max(h0, max(h1, h2)), hc);
+
+      // Look for strong changes in cell depth, or shallow cell depths, as an indicator of near-wet-dry
+      hfactor= max(0., min(2.0*max(hmin,0.0)/max(hc,1.0e-06)-0.1, 
+                           min(2.0*max(hc,0.)/max(hmax,1.0e-06)-0.1, 1.0))
+                  );
+      //hfactor=1.0;
+      hfactor=min( max(hmin,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
+
+      //-----------------------------------
+      // stage
+      //-----------------------------------      
+      
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k
+      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle (centroids of neighbouring triangles)
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
+      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
+     
+      inv_area2 = 1.0/area2;
+      // Calculate the gradient of stage on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+    
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+    
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+      //beta_tmp = beta_w_dry*0. + (beta_w - beta_w_dry*0.) * hfactor;
+      //beta_tmp=1.0;
+    
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+      stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
+
+
+      //-----------------------------------
+      // height
+      //-----------------------------------
+       
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k
+      dq0 = height_centroid_values[k0] - height_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle (centroids of neighbouring triangles)
+      dq1 = height_centroid_values[k1] - height_centroid_values[k0];
+      dq2 = height_centroid_values[k2] - height_centroid_values[k0];
+     
+      inv_area2 = 1.0/area2;
+      // Calculate the gradient of height on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      // Calculate provisional jumps in height from the centroid 
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+    
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+    
+      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
+      //beta_tmp=beta_w; 
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+      height_edge_values[k3+0] = height_centroid_values[k] + dqv[0];
+      height_edge_values[k3+1] = height_centroid_values[k] + dqv[1];
+      height_edge_values[k3+2] = height_centroid_values[k] + dqv[2];
+
+
+      // REDEFINE hfactor for momentum terms -- make MORE first order
+      hfactor= max(0., min(1.6*max(hmin,0.0)/max(hc,1.0e-06)-0.5, 
+                           min(1.6*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))
+                  );
+      hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), hfactor);
+      
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------            
+
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k      
+      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
+      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
+      
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited      
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+      
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+
+      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+
+      for (i=0; i < 3; i++)
+      {
+        xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
+      }
+      
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------                  
+
+      // Calculate the difference between vertex 0 of the auxiliary 
+      // triangle and the centroid of triangle k
+      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
+      
+      // Calculate differentials between the vertices 
+      // of the auxiliary triangle
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
+      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
+      
+      // Calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a *= inv_area2;
+      b = dx1*dq2 - dx2*dq1;
+      b *= inv_area2;
+      
+      // Calculate provisional jumps in stage from the centroid 
+      // of triangle k to its vertices, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now we want to find min and max of the centroid and the 
+      // vertices of the auxiliary triangle and compute jumps 
+      // from the centroid to the min and max
+      //
+      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
+      
+      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
+
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_tmp);
+      
+      for (i=0;i<3;i++)
+      {
+        ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+      }
+      
+    } // End number_of_boundaries <=1 
+    else
+    {
+
+      //==============================================
+      // Number of boundaries == 2
+      //==============================================        
         
-        
-          // If no boundaries, auxiliary triangle is formed 
-          // from the centroids of the three neighbours
-          // If one boundary, auxiliary triangle is formed 
-          // from this centroid and its two neighbours
-          
-          k0 = surrogate_neighbours[k3];
-          k1 = surrogate_neighbours[k3 + 1];
-          k2 = surrogate_neighbours[k3 + 2];
-         
-          // Test to see whether we accept the surrogate neighbours 
-          // Note that if ki is replaced with k in more than 1 neighbour, then the
-          // triangle area will be zero, and a first order extrapolation will be
-          // used
-          //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
-          //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
-          // FIXME: Remove cell_wetness_scale if you don't need it
-          if(k2<0 || (cell_wetness_scale[k2]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
-              k2 = k ;
+      // One internal neighbour and gradient is in direction of the neighbour's centroid
+      
+      // Find the only internal neighbour (k1?)
+      for (k2 = k3; k2 < k3 + 3; k2++)
+      {
+      // Find internal neighbour of triangle k      
+      // k2 indexes the edges of triangle k   
+      
+          if (surrogate_neighbours[k2] != k)
+          {
+             break;
           }
-          //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
-          //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
-          if(k0 < 0 || (cell_wetness_scale[k0]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
-              k0 = k ;
-          }
-          //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
-          //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
-          //if(k1 < 0 || (cell_wetness_scale[k1]==0.0 && stage_centroid_values[k]<elevation_edge_values[3*k+1])){
-          if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
-              k1 = k ;
-          }
-
-          // Take note if the max neighbour bed elevation is greater than the min
-          // neighbour stage -- suggests a 'steep' bed relative to the flow
-          bedmax = max(elevation_centroid_values[k], 
-                       max(elevation_centroid_values[k0],
-                           max(elevation_centroid_values[k1], 
-                               elevation_centroid_values[k2])));
-          bedmin = min(elevation_centroid_values[k], 
-                       min(elevation_centroid_values[k0],
-                           min(elevation_centroid_values[k1], 
-                               elevation_centroid_values[k2])));
-          //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 
-          //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
-          //                   min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
-          //                       max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
-
-          // Get the auxiliary triangle's vertex coordinates 
-          // (really the centroids of neighbouring triangles)
-          coord_index = 2*k0;
-          x0 = centroid_coordinates[coord_index];
-          y0 = centroid_coordinates[coord_index+1];
-          
-          coord_index = 2*k1;
-          x1 = centroid_coordinates[coord_index];
-          y1 = centroid_coordinates[coord_index+1];
-          
-          coord_index = 2*k2;
-          x2 = centroid_coordinates[coord_index];
-          y2 = centroid_coordinates[coord_index+1];
-
-          // Store x- and y- differentials for the vertices 
-          // of the auxiliary triangle
-          dx1 = x1 - x0; 
-          dx2 = x2 - x0;
-          dy1 = y1 - y0; 
-          dy2 = y2 - y0;
-          
-          // Calculate 2*area of the auxiliary triangle
-          // The triangle is guaranteed to be counter-clockwise      
-          area2 = dy2*dx1 - dy1*dx2; 
-           
-          //// Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
-          if ((area2 <= 0.) )//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
-          {
-
-              // Isolated wet cell -- constant stage/depth extrapolation 
-              //stage_edge_values[k3]   = stage_centroid_values[k];
-              //stage_edge_values[k3+1] = stage_centroid_values[k];
-              //stage_edge_values[k3+2] = stage_centroid_values[k];
-
-              dk=height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
-              height_edge_values[k3] = dk;
-              height_edge_values[k3+1] = dk;
-              height_edge_values[k3+2] = dk;
-              
-              stage_edge_values[k3]   = elevation_centroid_values[k]+dk;
-              stage_edge_values[k3+1] = elevation_centroid_values[k]+dk;
-              stage_edge_values[k3+2] = elevation_centroid_values[k]+dk;
-              //stage_edge_values[k3] = elevation_edge_values[k3]+dk;
-              //stage_edge_values[k3+1] = elevation_edge_values[k3+1]+dk;
-              //stage_edge_values[k3+2] = elevation_edge_values[k3+2]+dk;
-
-              xmom_edge_values[k3]    = xmom_centroid_values[k];
-              xmom_edge_values[k3+1]  = xmom_centroid_values[k];
-              xmom_edge_values[k3+2]  = xmom_centroid_values[k];
-              ymom_edge_values[k3]    = ymom_centroid_values[k];
-              ymom_edge_values[k3+1]  = ymom_centroid_values[k];
-              ymom_edge_values[k3+2]  = ymom_centroid_values[k];
-
-              continue;
-          }  
-          
-          // Calculate heights of neighbouring cells
-          hc = stage_centroid_values[k]  - elevation_centroid_values[k];
-          h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
-          h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
-          h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-          //h0=0.;
-          //h1=0.;
-          //h2=0.;
-
-          //if(k!=k0) h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
-          //if(k!=k1) h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
-          //if(k!=k2) h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-          
-          hmin = min(min(h0, min(h1, h2)), hc);
-          hmax = max(max(h0, max(h1, h2)), hc);
-          //// GD FIXME: No longer needed? 
-          //hfactor = 0.0;
-          ////if (hmin > 0.001) 
-          //if (hmin > 0.) 
-          ////if (hc>0.0)
-          //{
-          //  hfactor = 1.0 ;//hmin/(hmin + 0.004);
-            //hfactor=hmin/(hmin + 0.004);
-          //}
-          //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(1.0*(max_elevation_edgevalue[k]-min_elevation_edgevalue[k]),minimum_allowed_height*1.)), 1.0);
-
-          // Look for strong changes in cell wetness as an indicator of near-wet-dry
-          hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height), 
-                       min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height), 1.0));
-          //hfactor= min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height), 1.0);
-          //hfactor=0.;
-          //if(hmin==hc) hfactor=0.;
-          //hfactor= min(2.0*max(hmin,0.0)/max(hc,minimum_allowed_height+bedmax-bedmin), 
-          //             min(2.0*max(hc,0.)/max(hmax,minimum_allowed_height+bedmax-bedmin), 1.0));
-
-          //hfactor=1.0;
-          //hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), 1.0);
-          //if(hc<minimum_allowed_height*10.) hfactor=0.;
-          //if(cell_wetness_scale[k0]==0 | cell_wetness_scale[k1]==0 | cell_wetness_scale[k2]==0) hfactor=0;
-          //hfactor=1.0;
-          
-          //-----------------------------------
-          // stage
-          //-----------------------------------      
-          
-          // Calculate the difference between vertex 0 of the auxiliary 
-          // triangle and the centroid of triangle k
-          dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
-          
-          // Calculate differentials between the vertices 
-          // of the auxiliary triangle (centroids of neighbouring triangles)
-          dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
-          dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
-         
-          inv_area2 = 1.0/area2;
-          // Calculate the gradient of stage on the auxiliary triangle
-          a = dy2*dq1 - dy1*dq2;
-          a *= inv_area2;
-          b = dx1*dq2 - dx2*dq1;
-          b *= inv_area2;
-          // Calculate provisional jumps in stage from the centroid 
-          // of triangle k to its vertices, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-        
-          // Now we want to find min and max of the centroid and the 
-          // vertices of the auxiliary triangle and compute jumps 
-          // from the centroid to the min and max
-          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
-        
-          beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
-          //beta_tmp = beta_w_dry*0. + (beta_w - beta_w_dry*0.) * hfactor;
-          //beta_tmp=1.0;
-        
-          
-          // Limit the gradient
-          limit_gradient(dqv, qmin, qmax, beta_tmp);
-
-          stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
-          stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
-          stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
-
-
-          //-----------------------------------
-          // height
-          //-----------------------------------
-           
-          // Calculate the difference between vertex 0 of the auxiliary 
-          // triangle and the centroid of triangle k
-          dq0 = height_centroid_values[k0] - height_centroid_values[k];
-          
-          // Calculate differentials between the vertices 
-          // of the auxiliary triangle (centroids of neighbouring triangles)
-          dq1 = height_centroid_values[k1] - height_centroid_values[k0];
-          dq2 = height_centroid_values[k2] - height_centroid_values[k0];
-         
-          inv_area2 = 1.0/area2;
-          // Calculate the gradient of height on the auxiliary triangle
-          a = dy2*dq1 - dy1*dq2;
-          a *= inv_area2;
-          b = dx1*dq2 - dx2*dq1;
-          b *= inv_area2;
-          // Calculate provisional jumps in height from the centroid 
-          // of triangle k to its vertices, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-        
-          // Now we want to find min and max of the centroid and the 
-          // vertices of the auxiliary triangle and compute jumps 
-          // from the centroid to the min and max
-          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
-        
-          beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
-          //beta_tmp=beta_w; 
-          
-          // Limit the gradient
-          limit_gradient(dqv, qmin, qmax, beta_tmp);
-
-          height_edge_values[k3+0] = height_centroid_values[k] + dqv[0];
-          height_edge_values[k3+1] = height_centroid_values[k] + dqv[1];
-          height_edge_values[k3+2] = height_centroid_values[k] + dqv[2];
-
-          //-----------------------------------
-          // xmomentum
-          //-----------------------------------            
-
-          // Calculate the difference between vertex 0 of the auxiliary 
-          // triangle and the centroid of triangle k      
-          dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
-          
-          // Calculate differentials between the vertices 
-          // of the auxiliary triangle
-          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
-          dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
-          
-          // Calculate the gradient of xmom on the auxiliary triangle
-          a = dy2*dq1 - dy1*dq2;
-          a *= inv_area2;
-          b = dx1*dq2 - dx2*dq1;
-          b *= inv_area2;
-          
-          // Calculate provisional jumps in stage from the centroid 
-          // of triangle k to its vertices, to be limited      
-          dqv[0] = a*dxv0+b*dyv0;
-          dqv[1] = a*dxv1+b*dyv1;
-          dqv[2] = a*dxv2+b*dyv2;
-          
-          // Now we want to find min and max of the centroid and the 
-          // vertices of the auxiliary triangle and compute jumps 
-          // from the centroid to the min and max
-          //
-          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
-
-          beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
-
-          // Limit the gradient
-          //if((k!=k0)&(k!=k1)&(k!=k2))
-          //if(stagemin>=bedmax){
-          //if((count_wet_neighbours[k]>0)){ // && 
-            // (cell_wetness_scale[k]==1.0)){
-          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
-          
-          limit_gradient(dqv, qmin, qmax, beta_tmp);
-          //}else{
-          //  dqv[0]=0.;
-          //  dqv[1]=0.;
-          //  dqv[2]=0.;
-          ////  limit_gradient(dqv, qmin, qmax, 0.);
-          //} 
-          //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
-          
-
-          for (i=0; i < 3; i++)
-          {
-            xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
-          }
-          
-          //-----------------------------------
-          // ymomentum
-          //-----------------------------------                  
-
-          // Calculate the difference between vertex 0 of the auxiliary 
-          // triangle and the centroid of triangle k
-          dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
-          
-          // Calculate differentials between the vertices 
-          // of the auxiliary triangle
-          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
-          dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
-          
-          // Calculate the gradient of xmom on the auxiliary triangle
-          a = dy2*dq1 - dy1*dq2;
-          a *= inv_area2;
-          b = dx1*dq2 - dx2*dq1;
-          b *= inv_area2;
-          
-          // Calculate provisional jumps in stage from the centroid 
-          // of triangle k to its vertices, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-          
-          // Now we want to find min and max of the centroid and the 
-          // vertices of the auxiliary triangle and compute jumps 
-          // from the centroid to the min and max
-          //
-          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
-          
-          beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
-
-          // Limit the gradient
-          //if(stagemin>=bedmax){
-          //if((count_wet_neighbours[k]>0)){//&&
-             //(cell_wetness_scale[k]==1.0)){
-            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
-          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
-            limit_gradient(dqv, qmin, qmax, beta_tmp);
-          //}else{
-          //  dqv[0]=0.;
-          //  dqv[1]=0.;
-          //  dqv[2]=0.;
-          //} 
-          
-          for (i=0;i<3;i++)
-          {
-            ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
-          }
-          
-        } // End number_of_boundaries <=1 
-        else
-        {
-
-          //==============================================
-          // Number of boundaries == 2
-          //==============================================        
-            
-          // One internal neighbour and gradient is in direction of the neighbour's centroid
-          
-          // Find the only internal neighbour (k1?)
-          for (k2 = k3; k2 < k3 + 3; k2++)
-          {
-          // Find internal neighbour of triangle k      
-          // k2 indexes the edges of triangle k   
-          
-              if (surrogate_neighbours[k2] != k)
+      }
+      
+      if ((k2 == k3 + 3)) 
+      {
+        // If we didn't find an internal neighbour
+        report_python_error(AT, "Internal neighbour not found");
+        return -1;
+      }
+      
+      k1 = surrogate_neighbours[k2];
+      
+      // The coordinates of the triangle are already (x,y). 
+      // Get centroid of the neighbour (x1,y1)
+      coord_index = 2*k1;
+      x1 = centroid_coordinates[coord_index];
+      y1 = centroid_coordinates[coord_index + 1];
+      
+      // Compute x- and y- distances between the centroid of 
+      // triangle k and that of its neighbour
+      dx1 = x1 - x; 
+      dy1 = y1 - y;
+      
+      // Set area2 as the square of the distance
+      area2 = dx1*dx1 + dy1*dy1;
+      
+      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) 
+      // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+      // respectively correspond to the x- and y- gradients 
+      // of the conserved quantities
+      dx2 = 1.0/area2;
+      dy2 = dx2*dy1;
+      dx2 *= dx1;
+      
+      
+      //-----------------------------------
+      // stage
+      //-----------------------------------            
+
+      // Compute differentials
+      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
+      stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
+      stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
+
+      //-----------------------------------
+      // height
+      //-----------------------------------
+      
+      // Compute differentials
+      dq1 = height_centroid_values[k1] - height_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      height_edge_values[k3] = height_centroid_values[k] + dqv[0];
+      height_edge_values[k3 + 1] = height_centroid_values[k] + dqv[1];
+      height_edge_values[k3 + 2] = height_centroid_values[k] + dqv[2];
+
+      //-----------------------------------
+      // xmomentum
+      //-----------------------------------                        
+      
+      // Compute differentials
+      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0+b*dyv0;
+      dqv[1] = a*dxv1+b*dyv1;
+      dqv[2] = a*dxv2+b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1 >= 0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      }
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient      
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      //for (i=0;i<3;i++)
+      //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
+      //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
+      //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
+      
+      for (i = 0; i < 3;i++)
+      {
+          xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
+      }
+      
+      //-----------------------------------
+      // ymomentum
+      //-----------------------------------                        
+
+      // Compute differentials
+      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
+      
+      // Calculate the gradient between the centroid of triangle k 
+      // and that of its neighbour
+      a = dq1*dx2;
+      b = dq1*dy2;
+      
+      // Calculate provisional edge jumps, to be limited
+      dqv[0] = a*dxv0 + b*dyv0;
+      dqv[1] = a*dxv1 + b*dyv1;
+      dqv[2] = a*dxv2 + b*dyv2;
+      
+      // Now limit the jumps
+      if (dq1>=0.0)
+      {
+        qmin = 0.0;
+        qmax = dq1;
+      } 
+      else
+      {
+        qmin = dq1;
+        qmax = 0.0;
+      }
+      
+      // Limit the gradient
+      limit_gradient(dqv, qmin, qmax, beta_w);
+      
+      for (i=0;i<3;i++)
               {
-                 break;
+              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
               }
-          }
-          
-          if ((k2 == k3 + 3)) 
-          {
-            // If we didn't find an internal neighbour
-            report_python_error(AT, "Internal neighbour not found");
-            return -1;
-          }
-          
-          k1 = surrogate_neighbours[k2];
-          
-          // The coordinates of the triangle are already (x,y). 
-          // Get centroid of the neighbour (x1,y1)
-          coord_index = 2*k1;
-          x1 = centroid_coordinates[coord_index];
-          y1 = centroid_coordinates[coord_index + 1];
-          
-          // Compute x- and y- distances between the centroid of 
-          // triangle k and that of its neighbour
-          dx1 = x1 - x; 
-          dy1 = y1 - y;
-          
-          // Set area2 as the square of the distance
-          area2 = dx1*dx1 + dy1*dy1;
-          
-          // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) 
-          // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
-          // respectively correspond to the x- and y- gradients 
-          // of the conserved quantities
-          dx2 = 1.0/area2;
-          dy2 = dx2*dy1;
-          dx2 *= dx1;
-          
-          
-          //-----------------------------------
-          // stage
-          //-----------------------------------            
-
-          // Compute differentials
-          dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
-          
-          // Calculate the gradient between the centroid of triangle k 
-          // and that of its neighbour
-          a = dq1*dx2;
-          b = dq1*dy2;
-          
-          // Calculate provisional edge jumps, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-          
-          // Now limit the jumps
-          if (dq1>=0.0)
-          {
-            qmin=0.0;
-            qmax=dq1;
-          }
-          else
-          {
-            qmin = dq1;
-            qmax = 0.0;
-          }
-          
-          // Limit the gradient
-          limit_gradient(dqv, qmin, qmax, beta_w);
-          
-          stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
-          stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
-          stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
-
-          //-----------------------------------
-          // height
-          //-----------------------------------
-          
-          // Compute differentials
-          dq1 = height_centroid_values[k1] - height_centroid_values[k];
-          
-          // Calculate the gradient between the centroid of triangle k 
-          // and that of its neighbour
-          a = dq1*dx2;
-          b = dq1*dy2;
-          
-          // Calculate provisional edge jumps, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-          
-          // Now limit the jumps
-          if (dq1>=0.0)
-          {
-            qmin=0.0;
-            qmax=dq1;
-          }
-          else
-          {
-            qmin = dq1;
-            qmax = 0.0;
-          }
-          
-          // Limit the gradient
-          limit_gradient(dqv, qmin, qmax, beta_w);
-          
-          height_edge_values[k3] = height_centroid_values[k] + dqv[0];
-          height_edge_values[k3 + 1] = height_centroid_values[k] + dqv[1];
-          height_edge_values[k3 + 2] = height_centroid_values[k] + dqv[2];
-
-          //-----------------------------------
-          // xmomentum
-          //-----------------------------------                        
-          
-          // Compute differentials
-          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
-          
-          // Calculate the gradient between the centroid of triangle k 
-          // and that of its neighbour
-          a = dq1*dx2;
-          b = dq1*dy2;
-          
-          // Calculate provisional edge jumps, to be limited
-          dqv[0] = a*dxv0+b*dyv0;
-          dqv[1] = a*dxv1+b*dyv1;
-          dqv[2] = a*dxv2+b*dyv2;
-          
-          // Now limit the jumps
-          if (dq1 >= 0.0)
-          {
-            qmin = 0.0;
-            qmax = dq1;
-          }
-          else
-          {
-            qmin = dq1;
-            qmax = 0.0;
-          }
-          
-          // Limit the gradient      
-          limit_gradient(dqv, qmin, qmax, beta_w);
-          
-          //for (i=0;i<3;i++)
-          //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
-          //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
-          //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
-          
-          for (i = 0; i < 3;i++)
-          {
-              xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
-          }
-          
-          //-----------------------------------
-          // ymomentum
-          //-----------------------------------                        
-
-          // Compute differentials
-          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
-          
-          // Calculate the gradient between the centroid of triangle k 
-          // and that of its neighbour
-          a = dq1*dx2;
-          b = dq1*dy2;
-          
-          // Calculate provisional edge jumps, to be limited
-          dqv[0] = a*dxv0 + b*dyv0;
-          dqv[1] = a*dxv1 + b*dyv1;
-          dqv[2] = a*dxv2 + b*dyv2;
-          
-          // Now limit the jumps
-          if (dq1>=0.0)
-          {
-            qmin = 0.0;
-            qmax = dq1;
-          } 
-          else
-          {
-            qmin = dq1;
-            qmax = 0.0;
-          }
-          
-          // Limit the gradient
-          limit_gradient(dqv, qmin, qmax, beta_w);
-          
-          for (i=0;i<3;i++)
-                  {
-                  ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
-                  }
-        } // else [number_of_boundaries==2]
-      } // for k=0 to number_of_elements-1
-  //} 
-
-  //for(k=0;k<number_of_elements; k++){
-  //    // Hack for 'dry' cells
-  //    // Make sure edge value is not greater than neighbour edge value
-  //    //if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
-  //    if(cell_wetness_scale[k]==0.0){
-  //      for(i=0; i<3;i++){
-  //          if(surrogate_neighbours[3*k+i] >= 0){
-  //              ktmp=3*surrogate_neighbours[3*k+i] + neighbour_edges[3*k+i];
-  //          //if(cell_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
-  //              stage_edge_values[3*k+i] = min(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
-  //          }
-  //      }
-  //    }
-  //}
+    } // else [number_of_boundaries==2]
+  } // for k=0 to number_of_elements-1
+
 
   // Compute vertex values of quantities

trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

@@ -11 +11 @@
 
 
-p_dev = util2.get_output('dam_break_20131030_174145/dam_break.sww', 0.001)
+p_dev = util2.get_output('dam_break_20131103_204343/dam_break.sww', 0.001)
 p2_dev=util2.get_centroids(p_dev, velocity_extrapolation=True)