Index: /inundation/ga/storm_surge/pyvolution/cornell_room.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/cornell_room.py	(revision 487)
+++ /inundation/ga/storm_surge/pyvolution/cornell_room.py	(revision 488)
@@ -12,16 +12,47 @@
      Transmissive_boundary, Time_boundary, Constant_height
 
-from mesh_factory import from_polyfile
+from mesh_factory import from_polyfile, rectangular
 from Numeric import array
+from math import sqrt
+from least_squares import Interpolation
     
 
 print 'Creating domain'
-#points, triangles, values = from_polyfile('cornell_room_medres')
-points, triangles, values = from_polyfile('hires2')
+#data_points, _, data_values = from_polyfile('cornell_room_medres')
+#points, triangles, values = from_polyfile('hires2')
+data_points, _, data_values = from_polyfile('hires2')
 
 
+#Regrid onto numerically stable mesh
+#
+#Compute regular mesh based on resolution and extent of data
+data_points = array(data_points)
+pmax = max(data_points)
+pmin = min(data_points)
+
+M = len(data_points)
+
+N = int(0.8*sqrt(M))
+
+#print N
+
+mesh_points, vertices, boundary = rectangular(N, N, 
+					      len1=pmax[0]-pmin[0],
+                                              len2=pmax[1]-pmin[1],
+					      origin = pmin)
+					      
+
+#Compute smooth surface on new mesh based on values from old (regrid)
+print 'Interp'
+interp = Interpolation(mesh_points, vertices, data_points, alpha=0.1)
+mesh_values = interp.fit(data_values)
+print 'Len mesh values', len(mesh_values)
+print 'Len mesh points', len(mesh_points)
+
+    
 #Create shallow water domain
-domain = Domain(points, triangles)
-    
+print 'Creating domain'
+domain = Domain(mesh_points, vertices) #, boundary)
+
 domain.check_integrity()
 domain.default_order = 2
@@ -47,5 +78,5 @@
 
 print 'Field values'
-domain.set_quantity('elevation', values)
+domain.set_quantity('elevation', mesh_values)
 domain.set_quantity('friction', manning)
 
Index: /inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 487)
+++ /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 488)
@@ -3,5 +3,5 @@
    Implements a penalised least-squares fit and associated interpolations.
 
-   The panalty term (or smoothing term) is controlled by the smoothing
+   The penalty term (or smoothing term) is controlled by the smoothing
    parameter alpha.
    With a value of alpha=0, the fit function will attempt
@@ -20,5 +20,5 @@
 
 #FIXME (Ole): Currently datapoints outside the triangular mesh are ignored.
-#             Is there a clean way of inlcuding them?
+#             Is there a clean way of including them?
 
 
@@ -152,12 +152,14 @@
         Inputs:
         
-          vertex_coordinates: List of coordinate pairs [xi, eta] of points
-          constituting mesh (or a an m x 2 Numeric array)
+          vertex_coordinates: List of coordinate pairs [xi, eta] of 
+	  points constituting mesh (or a an m x 2 Numeric array)
         
           triangles: List of 3-tuples (or a Numeric array) of
           integers representing indices of all vertices in the mesh.
 
-          point_coordinates: List of coordinate pairs [x, y] of data points
-          (or an nx2 Numeric array)
+          point_coordinates: List of coordinate pairs [x, y] of 
+	  data points (or an nx2 Numeric array)
+	  If point_coordinates is absent, only smoothing matrix will 
+	  be built 
 
           alpha: Smoothing parameter
@@ -504,6 +506,6 @@
     def fit_points(self, z):
         """Like fit, but more robust when each point has two or more attributes
-	FIXME(Ole): The name fit_points doesn't carry any meaning for me.
-	How about something like fit_multiple or fit_columns?
+	FIXME (Ole): The name fit_points doesn't carry any meaning 
+	for me. How about something like fit_multiple or fit_columns?
         """
 	
@@ -529,5 +531,5 @@
         
     def interpolate(self, f):
-        """Compute predicted values at data points implied in self.A. 
+        """Evaluate smooth surface f at data points implied in self.A. 
 
         The mesh values representing a smooth surface are
@@ -546,20 +548,4 @@
 	
 	    
-     
-    #FIXME: We will need a method 'evaluate(self):' that will interpolate
-    #a computed surface living on the mesh onto a collection of
-    #arbitrary data points
-    #
-    #Precondition: self.fit(z) has stored its result in self.f.
-    #
-    #Input: data_points
-    #Algorithm:
-    #  1 Build a new temporary A matrix based on mesh and new data points
-    #  2 Apply it to self.f (return A*self.f)
-    #
-    # ON
-
-
-
 #-------------------------------------------------------------
 if __name__ == "__main__":
Index: /inundation/ga/storm_surge/pyvolution/mesh_factory.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 487)
+++ /inundation/ga/storm_surge/pyvolution/mesh_factory.py	(revision 488)
@@ -86,5 +86,5 @@
     points = []    #x, y
     values = []    #z
-    vertex_values = []    #Repeated z    
+    ##vertex_values = []    #Repeated z    
     triangles = [] #v0, v1, v2
     
@@ -172,8 +172,8 @@
                 
             triangles.append([j0, j1, j2])
-            vertex_values.append([values[j0], values[j1], values[j2]])
+            ##vertex_values.append([values[j0], values[j1], values[j2]])
         else:
             offending +=1 
 
     print 'Removed %d offending triangles out of %d' %(offending, len(lines))
-    return points, triangles, vertex_values
+    return points, triangles, values
