Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 3718)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 3719)
@@ -209,5 +209,5 @@
 
 
-	//MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
+	//MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
 	for (k=0;k<N;k++){
 		semi_implicit_update[k]=0.0;
Index: /anuga_core/source/anuga/examples/island.py
===================================================================
--- /anuga_core/source/anuga/examples/island.py	(revision 3718)
+++ /anuga_core/source/anuga/examples/island.py	(revision 3719)
@@ -33,7 +33,8 @@
                                            'top': [2],
                                            'left': [3]},
-                          maximum_triangle_area = 5,
+                          maximum_triangle_area = 1.0,
                           filename = 'island.msh' ,
-                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
+                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 1.0)]
+                          #interior_holes=[[50,25], [70,25], [70,75], [50,75]],
                           )
 
@@ -54,5 +55,6 @@
 # beta_h == 0.0 means that constant (1st order) gradients are introduced
 # on h. This is equivalent to the constant depth used previously.
-domain.beta_h = 0.9
+domain.beta_h     = 0.2
+domain.beta_w_dry = 0.2
 
 
@@ -78,5 +80,5 @@
 
 #domain.set_quantity('friction', 0.1)  #Honky dory
-domain.set_quantity('friction', 100)     #Creep
+domain.set_quantity('friction', 100.0)     #Creep
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
@@ -94,5 +96,5 @@
 domain.check_integrity()
 
-
+domain.initialise_visualiser()
 #------------------------------------------------------------------------------
 # Evolve system through time
Index: /anuga_core/source/anuga/mesh_engine/test_triangmodule.py
===================================================================
--- /anuga_core/source/anuga/mesh_engine/test_triangmodule.py	(revision 3718)
+++ /anuga_core/source/anuga/mesh_engine/test_triangmodule.py	(revision 3719)
@@ -211,11 +211,11 @@
 
     def testbad_hole(self):
-
+        
         holelist = [(9.0)]
         regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
+        
+        points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
+        pointattlist = [[],[],[],[]]
+        seglist = [(0,1),(1,3),(3,2),(2,0)]
         segattlist = [0,0,0,0]
         trilist = []
Index: /anuga_core/source/anuga/pmesh/mesh_interface.py
===================================================================
--- /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 3718)
+++ /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 3719)
@@ -25,4 +25,5 @@
                              filename=None,
                              interior_regions=None,
+                             interior_holes=None,
                              poly_geo_reference=None,
                              mesh_geo_reference=None,
@@ -43,4 +44,6 @@
     Interior_regions is a list of tuples consisting of (polygon, resolution)
     for each region to be separately refined.
+    
+    Interior_holes is a list of ploygons for each hole.
 
     This function does not allow segments to share points - use underlying
@@ -71,4 +74,5 @@
               'filename': filename,
               'interior_regions': interior_regions,
+              'interior_holes': interior_holes,
               'poly_geo_reference': poly_geo_reference,
               'mesh_geo_reference': mesh_geo_reference,
@@ -106,4 +110,5 @@
                               filename=None,
                               interior_regions=None,
+                              interior_holes=None,
                               poly_geo_reference=None,
                               mesh_geo_reference=None,
@@ -146,5 +151,14 @@
                 raise PolygonError, msg
 
-
+    if interior_holes is not None:        
+        # Test that all the interior polygons are inside the bounding_poly
+        for interior_polygon, res in interior_holes:
+            indices = inside_polygon(interior_polygon, bounding_polygon,
+                                     closed = True, verbose = False)
+    
+            if len(indices) <> len(interior_polygon): 
+                msg = 'Interior polygon %s is outside bounding polygon: %s'\
+                      %(str(interior_polygon), str(bounding_polygon))
+                raise PolygonError, msg
 
     # Resolve geo referencing        
@@ -192,17 +206,32 @@
 
     # Do interior regions
+#    if interior_regions is not None:    
+#        for polygon, res in interior_regions:
+#            m.add_region_from_polygon(polygon,
+#                                      geo_reference=poly_geo_reference)
+#            # convert bounding poly to absolute values
+#            if poly_geo_reference is not None:
+#                polygon_absolute = \
+#                    poly_geo_reference.get_absolute(polygon)
+#            else:
+#                polygon_absolute = polygon
+#            inner_point = point_in_polygon(polygon_absolute)
+#            region = m.add_region(inner_point[0], inner_point[1])
+#            region.setMaxArea(res)
+            
+            
     if interior_regions is not None:    
         for polygon, res in interior_regions:
             m.add_region_from_polygon(polygon,
+                                      max_triangle_area = res,
                                       geo_reference=poly_geo_reference)
-            # convert bounding poly to absolute values
-            if poly_geo_reference is not None:
-                polygon_absolute = \
-                    poly_geo_reference.get_absolute(polygon)
-            else:
-                polygon_absolute = polygon
-            inner_point = point_in_polygon(polygon_absolute)
-            region = m.add_region(inner_point[0], inner_point[1])
-            region.setMaxArea(res)
+    
+            
+   # Do interior holes
+    if interior_holes is not None:    
+        for polygon, res in interior_holes:
+            m.add_hole_from_polygon(polygon,
+                                      geo_reference=poly_geo_reference)
+       
             
 
Index: /anuga_core/source/anuga/pmesh/test_mesh_interface.py
===================================================================
--- /anuga_core/source/anuga/pmesh/test_mesh_interface.py	(revision 3718)
+++ /anuga_core/source/anuga/pmesh/test_mesh_interface.py	(revision 3719)
@@ -137,5 +137,5 @@
         interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
 
-
+        interior_holes = None
 
         # Clear cache first
@@ -147,4 +147,5 @@
                'filename': None,
                'interior_regions': interior_regions,
+               'interior_holes': interior_holes,
                'poly_geo_reference': geo_ref_poly,
                'mesh_geo_reference': mesh_geo,
Index: /anuga_core/source/anuga/shallow_water/shallow_water_ext.c
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3718)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_ext.c	(revision 3719)
@@ -524,5 +524,5 @@
 
     if (hc < minimum_allowed_height) {
-
+    	
       //Old code: Set momentum to zero and ensure h is non negative
       //xmomc[k] = 0.0;
@@ -536,24 +536,25 @@
       //             has no Python equivalent
       if (hc <= 0.0) {
-	wc[k] = zc[k];
-	xmomc[k] = 0.0;
-	ymomc[k] = 0.0;
+      	wc[k] = zc[k];
+		xmomc[k] = 0.0;
+		ymomc[k] = 0.0;
       } else {
         //Reduce excessive speeds derived from division by small hc
+        
         u = xmomc[k]/hc;
-	if (fabs(u) > maximum_allowed_speed) {
-	  reduced_speed = maximum_allowed_speed * u/fabs(u);
-	  //printf("Speed (u) has been reduced from %.3f to %.3f\n",
-	  //	 u, reduced_speed);
-	  xmomc[k] = reduced_speed * hc;
-	}
+		if (fabs(u) > maximum_allowed_speed) {
+	  		reduced_speed = maximum_allowed_speed * u/fabs(u);
+	  		//printf("Speed (u) has been reduced from %.3f to %.3f\n",
+	  		//	 u, reduced_speed);
+	  		xmomc[k] = reduced_speed * hc;
+		}
 
         v = ymomc[k]/hc;
-	if (fabs(v) > maximum_allowed_speed) {
-	  reduced_speed = maximum_allowed_speed * v/fabs(v);
-	  //printf("Speed (v) has been reduced from %.3f to %.3f\n",
-	  //	 v, reduced_speed);
-	  ymomc[k] = reduced_speed * hc;
-	}
+		if (fabs(v) > maximum_allowed_speed) {
+	  		reduced_speed = maximum_allowed_speed * v/fabs(v);
+	  		//printf("Speed (v) has been reduced from %.3f to %.3f\n",
+	  		//	 v, reduced_speed);
+	  		ymomc[k] = reduced_speed * hc;
+		}
       }
     }
