Changeset 488


Ignore:
Timestamp:
Nov 5, 2004, 4:42:19 PM (19 years ago)
Author:
ole
Message:

Played with regridding of Cornell data

Location:
inundation/ga/storm_surge/pyvolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/cornell_room.py

    r320 r488  
    1212     Transmissive_boundary, Time_boundary, Constant_height
    1313
    14 from mesh_factory import from_polyfile
     14from mesh_factory import from_polyfile, rectangular
    1515from Numeric import array
     16from math import sqrt
     17from least_squares import Interpolation
    1618   
    1719
    1820print 'Creating domain'
    19 #points, triangles, values = from_polyfile('cornell_room_medres')
    20 points, triangles, values = from_polyfile('hires2')
     21#data_points, _, data_values = from_polyfile('cornell_room_medres')
     22#points, triangles, values = from_polyfile('hires2')
     23data_points, _, data_values = from_polyfile('hires2')
    2124
    2225
     26#Regrid onto numerically stable mesh
     27#
     28#Compute regular mesh based on resolution and extent of data
     29data_points = array(data_points)
     30pmax = max(data_points)
     31pmin = min(data_points)
     32
     33M = len(data_points)
     34
     35N = int(0.8*sqrt(M))
     36
     37#print N
     38
     39mesh_points, vertices, boundary = rectangular(N, N,
     40                                              len1=pmax[0]-pmin[0],
     41                                              len2=pmax[1]-pmin[1],
     42                                              origin = pmin)
     43                                             
     44
     45#Compute smooth surface on new mesh based on values from old (regrid)
     46print 'Interp'
     47interp = Interpolation(mesh_points, vertices, data_points, alpha=0.1)
     48mesh_values = interp.fit(data_values)
     49print 'Len mesh values', len(mesh_values)
     50print 'Len mesh points', len(mesh_points)
     51
     52   
    2353#Create shallow water domain
    24 domain = Domain(points, triangles)
    25    
     54print 'Creating domain'
     55domain = Domain(mesh_points, vertices) #, boundary)
     56
    2657domain.check_integrity()
    2758domain.default_order = 2
     
    4778
    4879print 'Field values'
    49 domain.set_quantity('elevation', values)
     80domain.set_quantity('elevation', mesh_values)
    5081domain.set_quantity('friction', manning)
    5182
  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r485 r488  
    33   Implements a penalised least-squares fit and associated interpolations.
    44
    5    The panalty term (or smoothing term) is controlled by the smoothing
     5   The penalty term (or smoothing term) is controlled by the smoothing
    66   parameter alpha.
    77   With a value of alpha=0, the fit function will attempt
     
    2020
    2121#FIXME (Ole): Currently datapoints outside the triangular mesh are ignored.
    22 #             Is there a clean way of inlcuding them?
     22#             Is there a clean way of including them?
    2323
    2424
     
    152152        Inputs:
    153153       
    154           vertex_coordinates: List of coordinate pairs [xi, eta] of points
    155           constituting mesh (or a an m x 2 Numeric array)
     154          vertex_coordinates: List of coordinate pairs [xi, eta] of
     155          points constituting mesh (or a an m x 2 Numeric array)
    156156       
    157157          triangles: List of 3-tuples (or a Numeric array) of
    158158          integers representing indices of all vertices in the mesh.
    159159
    160           point_coordinates: List of coordinate pairs [x, y] of data points
    161           (or an nx2 Numeric array)
     160          point_coordinates: List of coordinate pairs [x, y] of
     161          data points (or an nx2 Numeric array)
     162          If point_coordinates is absent, only smoothing matrix will
     163          be built
    162164
    163165          alpha: Smoothing parameter
     
    504506    def fit_points(self, z):
    505507        """Like fit, but more robust when each point has two or more attributes
    506         FIXME(Ole): The name fit_points doesn't carry any meaning for me.
    507         How about something like fit_multiple or fit_columns?
     508        FIXME (Ole): The name fit_points doesn't carry any meaning
     509        for me. How about something like fit_multiple or fit_columns?
    508510        """
    509511       
     
    529531       
    530532    def interpolate(self, f):
    531         """Compute predicted values at data points implied in self.A.
     533        """Evaluate smooth surface f at data points implied in self.A.
    532534
    533535        The mesh values representing a smooth surface are
     
    546548       
    547549           
    548      
    549     #FIXME: We will need a method 'evaluate(self):' that will interpolate
    550     #a computed surface living on the mesh onto a collection of
    551     #arbitrary data points
    552     #
    553     #Precondition: self.fit(z) has stored its result in self.f.
    554     #
    555     #Input: data_points
    556     #Algorithm:
    557     #  1 Build a new temporary A matrix based on mesh and new data points
    558     #  2 Apply it to self.f (return A*self.f)
    559     #
    560     # ON
    561 
    562 
    563 
    564550#-------------------------------------------------------------
    565551if __name__ == "__main__":
  • inundation/ga/storm_surge/pyvolution/mesh_factory.py

    r486 r488  
    8686    points = []    #x, y
    8787    values = []    #z
    88     vertex_values = []    #Repeated z   
     88    ##vertex_values = []    #Repeated z   
    8989    triangles = [] #v0, v1, v2
    9090   
     
    172172               
    173173            triangles.append([j0, j1, j2])
    174             vertex_values.append([values[j0], values[j1], values[j2]])
     174            ##vertex_values.append([values[j0], values[j1], values[j2]])
    175175        else:
    176176            offending +=1
    177177
    178178    print 'Removed %d offending triangles out of %d' %(offending, len(lines))
    179     return points, triangles, vertex_values
     179    return points, triangles, values
Note: See TracChangeset for help on using the changeset viewer.