Index: inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/mesh.py	(revision 485)
+++ inundation/ga/storm_surge/pyvolution/mesh.py	(revision 486)
@@ -432,12 +432,3 @@
         #    assert self.neighbours[id,edge] < 0
 
-#FIXME: May get rid of
-def rectangular_mesh(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
-    from mesh_factory import rectangular
-    points, triangles, boundary = rectangular(m, n, len1, len2, origin)
-
-    return Mesh(points, triangles, boundary)
-
-
-#FIXME: Get oblique and contracting and circular meshes from Chris
-
+
Index: inundation/ga/storm_surge/pyvolution/mesh_factory.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 485)
+++ inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 486)
@@ -62,4 +62,6 @@
     return points, elements, boundary       
 
+
+#FIXME: Get oblique and contracting and circular meshes from Chris
 
 
Index: inundation/ga/storm_surge/pyvolution/test_domain.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 485)
+++ inundation/ga/storm_surge/pyvolution/test_domain.py	(revision 486)
@@ -6,5 +6,5 @@
 from domain import *
 from config import epsilon
-from Numeric import allclose, array
+from Numeric import allclose, array, ones, Float
 
 
@@ -266,15 +266,16 @@
         domain.update_conserved_quantities()
 
-
-        #FIXME: Update this test once I understand the semi_implicit scheme    
-        #x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
-        #x /= array( [.9,.9,.9,.9] )
-        #
-        #for name in domain.conserved_quantities:
-        #    assert allclose(domain.quantities[name].centroid_values, x)   
+        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
+        denom = ones(4, Float)-domain.timestep*sem
+
+        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
+        x /= denom
+
+        for name in domain.conserved_quantities:
+            assert allclose(domain.quantities[name].centroid_values, x)   
 
 
     def testz_set_region(self):
-        """Domain implements a default first order gradient limiter
+        """Set quantities for sub region
         """
         
@@ -326,10 +327,10 @@
         domain.set_region([add_to_verts])
         
-        self.failUnless( domain.test == "Mound",
+        self.failUnless(domain.test == "Mound",
                         'set region failed')
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'testz')
+    suite = unittest.makeSuite(TestCase,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: inundation/ga/storm_surge/pyvolution/test_mesh.py
===================================================================
--- inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 485)
+++ inundation/ga/storm_surge/pyvolution/test_mesh.py	(revision 486)
@@ -7,4 +7,5 @@
 
 from mesh import *
+from mesh_factory import rectangular
 from config import epsilon
 from Numeric import allclose, array
@@ -343,5 +344,7 @@
         M=1
         N=1
-        mesh = rectangular_mesh(M, N)
+
+        points, vertices, boundary = rectangular(M, N)
+	mesh = Mesh(points, vertices, boundary)
 
         #Test that points are arranged in a counter clock wise order        
@@ -350,6 +353,7 @@
         M=2
         N=2
-        mesh = rectangular_mesh(M, N)
-        
+        points, vertices, boundary = rectangular(M, N)
+	mesh = Mesh(points, vertices, boundary)
+	
         #Test that points are arranged in a counter clock wise order        
         mesh.check_integrity()
@@ -363,5 +367,7 @@
         
 
-        mesh = rectangular_mesh(4, 4)
+        points, vertices, boundary = rectangular(4, 4)
+	mesh = Mesh(points, vertices, boundary)
+	
 
         #Test that points are arranged in a counter clock wise order        
@@ -389,5 +395,7 @@
         len1 = 100.0
         len2 = 17.0
-        mesh = rectangular_mesh(M, N, len1, len2)
+	
+        points, vertices, boundary = rectangular(M, N, len1, len2)
+	mesh = Mesh(points, vertices, boundary)
 
         assert len(mesh) == 2*M*N        
@@ -408,5 +416,8 @@
         #Check that integers don't cause trouble
         N = 16
-        mesh = rectangular_mesh(2*N, N, len1=10, len2=10)
+	
+        points, vertices, boundary = rectangular(2*N, N, len1=10, len2=10)
+	mesh = Mesh(points, vertices, boundary)
+
 
 
