Changeset 488

Timestamp:

Nov 5, 2004, 4:42:19 PM (20 years ago)

Author:

ole

Message:

Played with regridding of Cornell data

Location:

inundation/ga/storm_surge/pyvolution

Files:

• cornell_room.py (modified) (2 diffs)
• least_squares.py (modified) (6 diffs)
• mesh_factory.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/cornell_room.py

@@ -12 +12 @@
      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 +78 @@
 
 print 'Field values'
-domain.set_quantity('elevation', values)
+domain.set_quantity('elevation', mesh_values)
 domain.set_quantity('friction', manning)
 

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -3 +3 @@
    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 +20 @@
 
 #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 +152 @@
         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 +506 @@
     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 +531 @@
         
     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 +548 @@
         
             
-     
-    #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__":

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -86 +86 @@
     points = []    #x, y
     values = []    #z
-    vertex_values = []    #Repeated z    
+    ##vertex_values = []    #Repeated z    
     triangles = [] #v0, v1, v2
     
@@ -172 +172 @@
                 
             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