Index: /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 8499)
+++ /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 8500)
@@ -64,4 +64,5 @@
                  triangles,
                  geo_reference=None,
+                 use_inscribed_circle=False,
                  verbose=False):
         """Build triangular 2d mesh from nodes and triangle information
@@ -82,6 +83,7 @@
         """
 
-        if verbose: log.critical('General_mesh: Building basic mesh structure '
-                                 'in ANUGA domain')
+        if verbose: log.critical('General_mesh: Building basic mesh structure')
+
+        self.use_inscribed_circle = use_inscribed_circle
          
         self.triangles = num.array(triangles, num.int)
@@ -130,4 +132,10 @@
 
         # Get x,y coordinates for all triangle vertices and store
+        self.centroid_coordinates = num.zeros((N, 2), num.float)
+
+        #Allocate space for geometric quantities
+        self.radii = num.zeros(N, num.float)
+
+        # Get x,y coordinates for all triangle vertices and store
         self.vertex_coordinates = V = self.compute_vertex_coordinates()
 
@@ -137,65 +145,220 @@
         # Initialise each triangle
         if verbose:
-            log.critical('General_mesh: Computing areas, normals '
-                         'and edgelengths')
-
-        for i in range(N):
-            if verbose and i % ((N+10)/10) == 0: log.critical('(%d/%d)' % (i, N))
-
-            x0, y0 = V[3*i, :]
-            x1, y1 = V[3*i+1, :]
-            x2, y2 = V[3*i+2, :]
-
-
-            i0 = self.triangles[i][0]
-            i1 = self.triangles[i][1]
-            i2 = self.triangles[i][2]
-
-            assert x0 == self.nodes[i0][0]
-            assert y0 == self.nodes[i0][1]
-
-            assert x1 == self.nodes[i1][0]
-            assert y1 == self.nodes[i1][1]
-
-            assert x2 == self.nodes[i2][0]
-            assert y2 == self.nodes[i2][1]
-            
-            # Area
-            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
-
-            msg = 'Triangle %g (%f,%f), (%f,%f), (%f, %f)' % (i,x0,y0,x1,y1,x2,y2)
-            msg += ' is degenerate:  area == %f' % self.areas[i]
-            assert self.areas[i] > 0.0, msg
-
-            # Normals
-            # The normal vectors
-            #   - point outward from each edge
-            #   - are orthogonal to the edge
-            #   - have unit length
-            #   - Are enumerated according to the opposite corner:
-            #     (First normal is associated with the edge opposite
-            #     the first vertex, etc)
-            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
-            n0 = num.array([x2-x1, y2-y1], num.float)
-            l0 = num.sqrt(num.sum(n0**2))
-
-            n1 = num.array([x0-x2, y0-y2], num.float)
-            l1 = num.sqrt(num.sum(n1**2))
-
-            n2 = num.array([x1-x0, y1-y0], num.float)
-            l2 = num.sqrt(num.sum(n2**2))
-
-            # Normalise
-            n0 /= l0
-            n1 /= l1
-            n2 /= l2
-
-            # Compute and store
-            self.normals[i, :] = [n0[1], -n0[0],
-                                  n1[1], -n1[0],
-                                  n2[1], -n2[0]]
-
-            # Edgelengths
-            self.edgelengths[i, :] = [l0, l1, l2]
+            log.critical('General_mesh: Computing areas, normals, '
+                         'edgelengths, centroids and radii')
+
+
+        # Calculate Areas
+        V0 = V[0:3*N:3, :]
+        V1 = V[1:3*N:3, :]
+        V2 = V[2:3*N:3, :]
+
+
+        # Area
+        x0 = V0[:,0]
+        y0 = V0[:,1]
+        x1 = V1[:,0]
+        y1 = V1[:,1]
+        x2 = V2[:,0]
+        y2 = V2[:,1]
+
+        self.areas[:] = num.abs( (x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
+
+        msg = 'Degenerate Triangle'
+        assert num.all(self.areas > 0.0), msg
+        #print V.shape, V0.shape, V1.shape, V2.shape
+
+#        #print E.shape, E[0:3*M:3, :].shape, E[1:3*M:3, :].shape, E[2:3*M:3, :].shape
+#        E[0:3*M:3, :] = 0.5*(V1+V2)
+#        E[1:3*M:3, :] = 0.5*(V2+V0)
+#        E[2:3*M:3, :] = 0.5*(V0+V1)
+
+        i0 = self.triangles[:,0]
+        i1 = self.triangles[:,1]
+        i2 = self.triangles[:,2]
+
+        assert num.allclose( x0, self.nodes[i0,0] )
+        assert num.allclose( y0, self.nodes[i0,1] )
+
+        assert num.allclose( x1, self.nodes[i1,0] )
+        assert num.allclose( y1, self.nodes[i1,1] )
+
+        assert num.allclose( x2, self.nodes[i2,0] )
+        assert num.allclose( y2, self.nodes[i2,1] )
+
+
+        xn0 = x2-x1
+        yn0 = y2-y1
+        l0 = num.sqrt(xn0**2 + yn0**2)
+
+        xn0 /= l0
+        yn0 /= l0
+
+        xn1 = x0-x2
+        yn1 = y0-y2
+        l1 = num.sqrt(xn1**2 + yn1**2)
+
+        xn1 /= l1
+        yn1 /= l1
+
+        xn2 = x1-x0
+        yn2 = y1-y0
+        l2 = num.sqrt(xn2**2 + yn2**2)
+
+        xn2 /= l2
+        yn2 /= l2
+
+        # Compute and store
+
+        self.normals[:,0] =  yn0
+        self.normals[:,1] = -xn0
+
+        self.normals[:,2] =  yn1
+        self.normals[:,3] = -xn1
+
+        self.normals[:,4] =  yn2
+        self.normals[:,5] = -xn2
+        
+        self.edgelengths[:,0] = l0
+        self.edgelengths[:,1] = l1
+        self.edgelengths[:,2] = l2
+
+        self.centroid_coordinates[:,0] = (x0 + x1 + x2)/3
+        self.centroid_coordinates[:,1] = (y0 + y1 + y2)/3
+
+
+
+        if self.use_inscribed_circle == False:
+            #OLD code. Computed radii may exceed that of an
+            #inscribed circle
+
+            #Midpoints
+            xm0 = (x1 + x2)/2
+            ym0 = (y1 + y2)/2
+
+            xm1 = (x2 + x0)/2
+            ym1 = (y2 + y0)/2
+
+            xm2 = (x0 + x1)/2
+            ym2 = (y0 + y1)/2
+
+
+            #The radius is the distance from the centroid of
+            #a triangle to the midpoint of the side of the triangle
+            #closest to the centroid
+
+            d0 = num.sqrt((self.centroid_coordinates[:,0] - xm0)**2 + (self.centroid_coordinates[:,1] - ym0)**2)
+            d1 = num.sqrt((self.centroid_coordinates[:,0] - xm1)**2 + (self.centroid_coordinates[:,1] - ym1)**2)
+            d2 = num.sqrt((self.centroid_coordinates[:,0] - xm2)**2 + (self.centroid_coordinates[:,1] - ym2)**2)
+
+
+            self.radii[:] = num.minimum(num.minimum(d0, d1), d2)
+
+        else:
+            #NEW code added by Peter Row. True radius
+            #of inscribed circle is computed
+
+            a = num.sqrt((x0-x1)**2+(y0-y1)**2)
+            b = num.sqrt((x1-x2)**2+(y1-y2)**2)
+            c = num.sqrt((x2-x0)**2+(y2-y0)**2)
+
+            self.radii[:]=2.0*self.areas/(a+b+c)
+
+
+
+#        for i in range(N):
+#            if verbose and i % ((N+10)/10) == 0: log.critical('(%d/%d)' % (i, N))
+#
+#            x0, y0 = V[3*i, :]
+#            x1, y1 = V[3*i+1, :]
+#            x2, y2 = V[3*i+2, :]
+#
+#
+#            i0 = self.triangles[i][0]
+#            i1 = self.triangles[i][1]
+#            i2 = self.triangles[i][2]
+#
+##            assert x0 == self.nodes[i0][0]
+##            assert y0 == self.nodes[i0][1]
+##
+##            assert x1 == self.nodes[i1][0]
+##            assert y1 == self.nodes[i1][1]
+##
+##            assert x2 == self.nodes[i2][0]
+##            assert y2 == self.nodes[i2][1]
+#
+##            # Area
+##            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
+##
+##            msg = 'Triangle %g (%f,%f), (%f,%f), (%f, %f)' % (i,x0,y0,x1,y1,x2,y2)
+##            msg += ' is degenerate:  area == %f' % self.areas[i]
+##            assert self.areas[i] > 0.0, msg
+#
+#            # Normals
+#            # The normal vectors
+#            #   - point outward from each edge
+#            #   - are orthogonal to the edge
+#            #   - have unit length
+#            #   - Are enumerated according to the opposite corner:
+#            #     (First normal is associated with the edge opposite
+#            #     the first vertex, etc)
+#            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
+#            n0 = num.array([x2-x1, y2-y1], num.float)
+#            l0 = num.sqrt(num.sum(n0**2))
+#
+#            n1 = num.array([x0-x2, y0-y2], num.float)
+#            l1 = num.sqrt(num.sum(n1**2))
+#
+#            n2 = num.array([x1-x0, y1-y0], num.float)
+#            l2 = num.sqrt(num.sum(n2**2))
+#
+#            # Normalise
+#            n0 /= l0
+#            n1 /= l1
+#            n2 /= l2
+#
+##            # Compute and store
+##            self.normals[i, :] = [n0[1], -n0[0],
+##                                  n1[1], -n1[0],
+##                                  n2[1], -n2[0]]
+#
+#            # Edgelengths
+#            #self.edgelengths[i, :] = [l0, l1, l2]
+#
+#
+#
+#            #Compute centroid
+##            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float)
+###            self.centroid_coordinates[i] = centroid
+##
+##
+##            if self.use_inscribed_circle == False:
+##                #OLD code. Computed radii may exceed that of an
+##                #inscribed circle
+##
+##                #Midpoints
+##                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
+##                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
+##                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
+##
+##                #The radius is the distance from the centroid of
+##                #a triangle to the midpoint of the side of the triangle
+##                #closest to the centroid
+##                d0 = num.sqrt(num.sum( (centroid-m0)**2 ))
+##                d1 = num.sqrt(num.sum( (centroid-m1)**2 ))
+##                d2 = num.sqrt(num.sum( (centroid-m2)**2 ))
+##
+##                #self.radii[i] = min(d0, d1, d2)
+##
+##            else:
+##                #NEW code added by Peter Row. True radius
+##                #of inscribed circle is computed
+##
+##                a = num.sqrt((x0-x1)**2+(y0-y1)**2)
+##                b = num.sqrt((x1-x2)**2+(y1-y2)**2)
+##                c = num.sqrt((x2-x0)**2+(y2-y0)**2)
+##
+##                self.radii[i]=2.0*self.areas[i]/(a+b+c)
+
 
         # Build structure listing which triangles belong to which node.
@@ -348,4 +511,12 @@
         M = self.number_of_triangles
         vertex_coordinates = num.zeros((3*M, 2), num.float)
+
+#        k0 = self.triangles[:,0]
+#        k1 = self.triangles[:,1]
+#        k2 = self.triangles[:,2]
+#
+#        vertex_coordinates[3*k0,:]   = self.nodes[k0,:]
+#        vertex_coordinates[3*k1+1,:] = self.nodes[k1,:]
+#        vertex_coordinates[3*k2+2,:] = self.nodes[k2,:]
 
         for i in range(M):
Index: /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py	(revision 8499)
+++ /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py	(revision 8500)
@@ -82,4 +82,6 @@
             coordinates = source
 
+        if verbose: log.critical('Domain: Initialising')
+
         # In case a filename has been specified, extract content
         if mesh_filename is not None:
@@ -131,5 +133,5 @@
         self.geo_reference = self.mesh.geo_reference
 
-        if verbose: log.critical('Initialising Domain')
+
 
         # List of quantity names entering the conservation equations
@@ -1347,5 +1349,5 @@
         # Update centroid values of conserved quantites to satisfy
         # special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts to ensure all centroid values are available
@@ -1401,4 +1403,7 @@
             # Update time
             self.set_time(initial_time + self.timestep)
+
+            # Update ghosts
+            self.update_ghosts()
 
             # Update vertex and edge values
@@ -1474,10 +1479,9 @@
         self.update_conserved_quantities()
 
-
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts
-        self.update_ghosts()
+        #self.update_ghosts()
 
 
@@ -1509,5 +1513,5 @@
 
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts
@@ -1548,8 +1552,8 @@
 
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts
-        self.update_ghosts()
+        #self.update_ghosts()
 
 
@@ -1582,5 +1586,5 @@
 
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts
@@ -1622,5 +1626,5 @@
 
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
 
         # Update ghosts
@@ -1659,11 +1663,12 @@
 
         # Update special conditions
-        self.update_special_conditions()
+        #self.update_special_conditions()
         
-        # Update ghosts
-        self.update_ghosts()
 
         # Set new time
         self.set_time(initial_time + self.timestep)
+
+        # Update ghosts
+        #self.update_ghosts()
 
 
@@ -1786,32 +1791,4 @@
 
 
-#        for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects):
-#            if B is None:
-#                log.critical('WARNING: Ignored boundary segment (None)')
-#            else:
-#                q_bdry = B.evaluate(vol_id, edge_id)
-#
-#                if len(q_bdry) == len(self.evolved_quantities):
-#                    # conserved and evolved quantities are the same
-#                    q_evol = q_bdry
-#                elif len(q_bdry) == len(self.conserved_quantities):
-#                    # boundary just returns conserved quantities
-#                    # Need to calculate all the evolved quantities
-#                    # Use default conversion
-#
-#                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
-#
-#                    q_evol = self.conserved_values_to_evolved_values \
-#                                                            (q_bdry, q_evol)
-#                else:
-#                    msg = 'Boundary must return array of either conserved'
-#                    msg += ' or evolved quantities'
-#                    raise Exception(msg)
-#
-#                for j, name in enumerate(self.evolved_quantities):
-#                    Q = self.quantities[name]
-#                    Q.boundary_values[i] = q_evol[j]
-
-
     def update_boundary(self):
         """Go through list of boundary objects and update boundary values
@@ -1823,5 +1800,4 @@
         """
 
-
         for tag in self.tag_boundary_cells:
 
@@ -1838,5 +1814,4 @@
         
 
-
     def compute_fluxes(self):
         msg = 'Method compute_fluxes must be overridden by Domain subclass'
@@ -1848,5 +1823,4 @@
         for operator in self.fractional_step_operators:
             operator()
-
 
 
@@ -1972,11 +1946,11 @@
                 num.put(Q_cv, Idg, num.take(Q_cv, Idf, axis=0))
 
-    def update_special_conditions(self):
-        """There may be a need to change the values of the conserved 
-        quantities to satisfy special conditions at the very lowest level
-        the fluid flow calculation
-        """
-        
-        pass
+#    def update_special_conditions(self):
+#        """There may be a need to change the values of the conserved
+#        quantities to satisfy special conditions at the very lowest level
+#        the fluid flow calculation
+#        """
+#
+#        pass
 
 
Index: /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 8499)
+++ /trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 8500)
@@ -85,19 +85,15 @@
         General_mesh.__init__(self, coordinates, triangles,
                               geo_reference=geo_reference,
+                              use_inscribed_circle=use_inscribed_circle,
                               verbose=verbose)
 
-        if verbose: log.critical('Initialising mesh')
+        if verbose: log.critical('Mesh: Initialising')
 
         N = len(self) #Number_of_triangles
 
-        self.use_inscribed_circle = use_inscribed_circle
-
-        #Allocate space for geometric quantities
-        self.centroid_coordinates = num.zeros((N, 2), num.float)
-
-        self.radii = num.zeros(N, num.float)
-
-        self.neighbours = num.zeros((N, 3), num.int)
-        self.neighbour_edges = num.zeros((N, 3), num.int)
+        # Allocate arrays for neighbour data
+
+        self.neighbours = -1*num.ones((N, 3), num.int)
+        self.neighbour_edges = -1*num.ones((N, 3), num.int)
         self.number_of_boundaries = num.zeros(N, num.int)
         self.surrogate_neighbours = num.zeros((N, 3), num.int)
@@ -106,57 +102,57 @@
         V = self.vertex_coordinates # Relative coordinates
 
-        #Initialise each triangle
-        if verbose: log.critical('Mesh: Computing centroids and radii')
-        for i in range(N):
-            if verbose and i % ((N+10)/10) == 0: log.critical('(%d/%d)' % (i, N))
-
-            x0, y0 = V[3*i, :]
-            x1, y1 = V[3*i+1, :]
-            x2, y2 = V[3*i+2, :]
-
-            #x0 = V[i, 0]; y0 = V[i, 1]
-            #x1 = V[i, 2]; y1 = V[i, 3]
-            #x2 = V[i, 4]; y2 = V[i, 5]
-
-            #Compute centroid
-            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float)
-            self.centroid_coordinates[i] = centroid
-
-
-            if self.use_inscribed_circle == False:
-                #OLD code. Computed radii may exceed that of an
-                #inscribed circle
-
-                #Midpoints
-                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
-                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
-                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
-
-                #The radius is the distance from the centroid of
-                #a triangle to the midpoint of the side of the triangle
-                #closest to the centroid
-                d0 = num.sqrt(num.sum( (centroid-m0)**2 ))
-                d1 = num.sqrt(num.sum( (centroid-m1)**2 ))
-                d2 = num.sqrt(num.sum( (centroid-m2)**2 ))
-
-                self.radii[i] = min(d0, d1, d2)
-
-            else:
-                #NEW code added by Peter Row. True radius
-                #of inscribed circle is computed
-
-                a = num.sqrt((x0-x1)**2+(y0-y1)**2)
-                b = num.sqrt((x1-x2)**2+(y1-y2)**2)
-                c = num.sqrt((x2-x0)**2+(y2-y0)**2)
-
-                self.radii[i]=2.0*self.areas[i]/(a+b+c)
-
-
-            #Initialise Neighbours (-1 means that it is a boundary neighbour)
-            self.neighbours[i, :] = [-1, -1, -1]
-
-            #Initialise edge ids of neighbours
-            #In case of boundaries this slot is not used
-            self.neighbour_edges[i, :] = [-1, -1, -1]
+#        #Initialise each triangle
+#        if verbose: log.critical('Mesh: Computing centroids and radii')
+#        for i in range(N):
+#            if verbose and i % ((N+10)/10) == 0: log.critical('(%d/%d)' % (i, N))
+#
+#            x0, y0 = V[3*i, :]
+#            x1, y1 = V[3*i+1, :]
+#            x2, y2 = V[3*i+2, :]
+#
+#            #x0 = V[i, 0]; y0 = V[i, 1]
+#            #x1 = V[i, 2]; y1 = V[i, 3]
+#            #x2 = V[i, 4]; y2 = V[i, 5]
+#
+#            #Compute centroid
+#            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float)
+#            self.centroid_coordinates[i] = centroid
+#
+#
+#            if self.use_inscribed_circle == False:
+#                #OLD code. Computed radii may exceed that of an
+#                #inscribed circle
+#
+#                #Midpoints
+#                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
+#                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
+#                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
+#
+#                #The radius is the distance from the centroid of
+#                #a triangle to the midpoint of the side of the triangle
+#                #closest to the centroid
+#                d0 = num.sqrt(num.sum( (centroid-m0)**2 ))
+#                d1 = num.sqrt(num.sum( (centroid-m1)**2 ))
+#                d2 = num.sqrt(num.sum( (centroid-m2)**2 ))
+#
+#                self.radii[i] = min(d0, d1, d2)
+#
+#            else:
+#                #NEW code added by Peter Row. True radius
+#                #of inscribed circle is computed
+#
+#                a = num.sqrt((x0-x1)**2+(y0-y1)**2)
+#                b = num.sqrt((x1-x2)**2+(y1-y2)**2)
+#                c = num.sqrt((x2-x0)**2+(y2-y0)**2)
+#
+#                self.radii[i]=2.0*self.areas[i]/(a+b+c)
+#
+#
+#            #Initialise Neighbours (-1 means that it is a boundary neighbour)
+#            self.neighbours[i, :] = [-1, -1, -1]
+#
+#            #Initialise edge ids of neighbours
+#            #In case of boundaries this slot is not used
+#            self.neighbour_edges[i, :] = [-1, -1, -1]
 
 
Index: /trunk/anuga_core/source/anuga/operators/erosion_operators.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/erosion_operators.py	(revision 8499)
+++ /trunk/anuga_core/source/anuga/operators/erosion_operators.py	(revision 8500)
@@ -68,5 +68,5 @@
 
         #-----------------------------------------
-        # Some stras for reporting
+        # Some extras for reporting
         #-----------------------------------------
         self.max_change = 0
Index: /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 8499)
+++ /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 8500)
@@ -240,5 +240,6 @@
         self.alpha_balance = alpha_balance
         self.tight_slope_limiters = tight_slope_limiters
-        self.optimise_dry_cells = int(optimise_dry_cells)
+        
+        self.set_use_optimise_dry_cells(optimise_dry_cells)
         self.set_extrapolate_velocity(extrapolate_velocity_second_order)
 
@@ -503,4 +504,15 @@
         elif flag is False:
             self.use_edge_limiter = False
+
+    def set_use_optimise_dry_cells(self, flag=True):
+        """ Try to optimize calculations where region is dry
+        """
+
+        if flag is True:
+            self.optimise_dry_cells = int(True)
+        elif flag is False:
+            self.optimise_dry_cells = int(False)
+
+
 
 
Index: /trunk/anuga_core/source/pypar_dist/source/pypar.py
===================================================================
--- /trunk/anuga_core/source/pypar_dist/source/pypar.py	(revision 8499)
+++ /trunk/anuga_core/source/pypar_dist/source/pypar.py	(revision 8500)
@@ -898,8 +898,8 @@
          reduce_array,\
          allreduce_array,\
-	 bsend_string, bsend_array, \
-	 mpi_alloc_and_attach, mpi_detach_and_dealloc, \
-	 mpi_alloc, mpi_dealloc, mpi_attach, mpi_detach, \
-	 string_push_for_alloc_and_attach, array_push_for_alloc_and_attach, \
+         bsend_string, bsend_array, \
+	     mpi_alloc_and_attach, mpi_detach_and_dealloc, \
+	     mpi_alloc, mpi_dealloc, mpi_attach, mpi_detach, \
+	     string_push_for_alloc_and_attach, array_push_for_alloc_and_attach, \
          MPI_ANY_TAG as any_tag, MPI_TAG_UB as max_tag,\
          MPI_ANY_SOURCE as any_source,\
