Changeset 1160


Ignore:
Timestamp:
Mar 29, 2005, 6:06:54 PM (19 years ago)
Author:
steve
Message:
 
Location:
inundation/ga/storm_surge
Files:
7 added
2 deleted
7 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/analytical solutions/Analytical solution_oblique_shock.py

    r774 r1160  
    1 """Example of shallow water wave equation
     1"""
     2Example of shallow water wave equation
    23consisting of an asymetrical converging channel.
    34
    4    Copyright 2004
    5    Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
    6    Geoscience Australia, ANU
    7    
     5Copyright 2004
     6Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
     7Geoscience Australia, ANU
     8
    89Specific methods pertaining to the 2D shallow water equation
    910are imported from shallow_water
     
    1314numerical vector named conserved_quantities.
    1415
    15 
    1616"""
    1717
    1818######################
    19 # Module imports 
     19# Module imports
    2020#
    2121
     
    2626import sys
    2727from os import sep
    28 sys.path.append('..'+s
    29 
    30 ep+'pyvolution')
     28sys.path.append('..'+sep+'pyvolution')
    3129
    3230
     
    5048
    5149points, elements, boundary = oblique(m, n, lenx, leny)
    52 domain = Domain(points, elements, boundary) 
     50domain = Domain(points, elements, boundary)
    5351
    5452# Order of solver
     
    9795for t in domain.evolve(yieldstep = 0.5, finaltime = 50):
    9896    domain.write_time()
    99    
     97
    10098print 'That took %.2f seconds' %(time.time()-t0)
    10199
     
    109107
    110108
    111    
     109
  • inundation/ga/storm_surge/pyvolution/cg_solve.py

    r1156 r1160  
    44import logging, logging.config
    55logger = logging.getLogger('cg_solve')
    6 logger.setLevel(logging.DEBUG)
     6logger.setLevel(logging.WARNING)
    77
    88try:
     
    1010except:
    1111    pass
    12 
    13 
    14 
    15 
    16 #
    17 
    18 
    19 
    20 
    2112
    2213def conjugate_gradient(A,b,x0=None,imax=10000,tol=1.0e-8,iprint=0):
     
    7970   #FIXME: Should this raise an exception?
    8071   if i==imax:
    81      print 'max number of iterations attained'
     72       logger.warning('max number of iterations attained')
    8273
    8374   return x
  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r1158 r1160  
    1212   A negative alpha is not allowed.
    1313   A typical value of alpha is 1.0e-6
    14          
     14
    1515
    1616   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
    17    Geoscience Australia, 2004.   
     17   Geoscience Australia, 2004.
    1818"""
    1919
     
    3232try:
    3333    from util import gradient
    34 except ImportError, e: 
     34except ImportError, e:
    3535    #FIXME reduce the dependency of modules in pyvolution
    3636    # Have util in a dir, working like load_mesh, and get rid of this
     
    3838        """
    3939        """
    40    
    41         det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)           
     40
     41        det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
    4242        a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
    4343        a /= det
    4444
    4545        b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
    46         b /= det           
     46        b /= det
    4747
    4848        return a, b
     
    5050
    5151DEFAULT_ALPHA = 0.001
    52    
     52
    5353def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
    5454                     alpha=DEFAULT_ALPHA, verbose= False,
     
    6161    point attributes to the mesh and write a mesh file with the
    6262    results.
    63    
     63
    6464
    6565    If data_origin is not None it is assumed to be
    6666    a 3-tuple with geo referenced
    6767    UTM coordinates (zone, easting, northing)
    68    
     68
    6969    mesh_origin is the same but refers to the input tsh file.
    7070    FIXME: When the tsh format contains it own origin, these parameters can go.
    7171    FIXME: And both origins should be obtained from the specified files.
    7272    """
    73    
     73
    7474    from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
    7575                 load_points_file, export_mesh_file, \
    7676                 concatinate_attributelist
    77    
     77
    7878    mesh_dict = mesh_file_to_mesh_dictionary(mesh_file)
    7979    vertex_coordinates = mesh_dict['vertices']
    8080    triangles = mesh_dict['triangles']
    81     if type(mesh_dict['vertex_attributes']) == ArrayType: 
     81    if type(mesh_dict['vertex_attributes']) == ArrayType:
    8282        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
    8383    else:
    84         old_point_attributes = mesh_dict['vertex_attributes'] 
    85        
     84        old_point_attributes = mesh_dict['vertex_attributes']
     85
    8686    if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
    8787        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
    88     else: 
     88    else:
    8989        old_title_list = mesh_dict['vertex_attribute_titles']
    90        
     90
    9191    if verbose: print 'tsh file %s loaded' %mesh_file
    92    
     92
    9393    # load in the .pts file
    9494    try:
     
    100100                                      delimiter = ' ',
    101101                                      verbose=verbose)
    102        
     102
    103103    point_coordinates = point_dict['pointlist']
    104104    title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
    105    
     105
    106106    if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
    107107        data_origin = point_dict['geo_reference'].get_origin()
    108108    else:
    109109        data_origin = (56, 0, 0) #FIXME(DSG-DSG)
    110        
     110
    111111    if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
    112112        mesh_origin = mesh_dict['geo_reference'].get_origin()
    113113    else:
    114114        mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
    115        
     115
    116116    if verbose: print "points file loaded"
    117117    if verbose:print "fitting to mesh"
     
    127127                    precrop = precrop)
    128128    if verbose: print "finished fitting to mesh"
    129    
     129
    130130    # convert array to list of lists
    131131    new_point_attributes = f.tolist()
     
    143143        mesh_dict['vertex_attribute_titles'] = title_list
    144144
    145     #FIXME (Ole): Remember to output mesh_origin as well   
     145    #FIXME (Ole): Remember to output mesh_origin as well
    146146    if verbose: print "exporting to file ",mesh_output_file
    147147    export_mesh_file(mesh_output_file, mesh_dict)
    148        
     148
    149149
    150150def fit_to_mesh(vertex_coordinates,
     
    162162    given data points with attributes.
    163163
    164          
     164
    165165        Inputs:
    166        
     166
    167167          vertex_coordinates: List of coordinate pairs [xi, eta] of points
    168168          constituting mesh (or a an m x 2 Numeric array)
    169        
     169
    170170          triangles: List of 3-tuples (or a Numeric array) of
    171171          integers representing indices of all vertices in the mesh.
     
    181181          UTM zone, easting and northing. If specified
    182182          point coordinates and vertex coordinates are assumed to be
    183           relative to their respective origins. 
    184          
     183          relative to their respective origins.
     184
    185185    """
    186186    interp = Interpolation(vertex_coordinates,
     
    193193                           mesh_origin = mesh_origin,
    194194                           precrop = precrop)
    195    
     195
    196196    vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
    197197    return vertex_attributes
    198198
    199    
    200    
     199
     200
    201201def pts2rectangular(pts_name, M, N, alpha = DEFAULT_ALPHA,
    202202                    verbose = False, reduction = 1, format = 'netcdf'):
    203     """Fits attributes from pts file to MxN rectangular mesh 
    204    
     203    """Fits attributes from pts file to MxN rectangular mesh
     204
    205205    Read pts file and create rectangular mesh of resolution MxN such that
    206     it covers all points specified in pts file. 
    207    
    208     FIXME: This may be a temporary function until we decide on 
     206    it covers all points specified in pts file.
     207
     208    FIXME: This may be a temporary function until we decide on
    209209    netcdf formats etc
    210210
    211211    FIXME: Uses elevation hardwired
    212     """   
    213    
     212    """
     213
    214214    import util, mesh_factory
    215215
     
    221221    elevation = attributes['elevation']  #Must be elevation
    222222    elevation = elevation[::reduction]
    223    
     223
    224224    if verbose: print 'Got %d data points' %len(points)
    225225
     
    229229    max_y = min_y = points[0][1]
    230230    for point in points[1:]:
    231         x = point[0] 
     231        x = point[0]
    232232        if x > max_x: max_x = x
    233         if x < min_x: min_x = x           
    234         y = point[1] 
     233        if x < min_x: min_x = x
     234        y = point[1]
    235235        if y > max_y: max_y = y
    236         if y < min_y: min_y = y           
    237    
     236        if y < min_y: min_y = y
     237
    238238    #Create appropriate mesh
    239239    vertex_coordinates, triangles, boundary =\
    240          mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y, 
     240         mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y,
    241241                                (min_x, min_y))
    242242
    243     #Fit attributes to mesh     
     243    #Fit attributes to mesh
    244244    vertex_attributes = fit_to_mesh(vertex_coordinates,
    245245                        triangles,
     
    248248
    249249
    250          
     250
    251251    return vertex_coordinates, triangles, boundary, vertex_attributes
    252                          
    253    
     252
     253
    254254
    255255class Interpolation:
     
    263263                 expand_search = True,
    264264                 max_points_per_cell = 30,
    265                  mesh_origin = None, 
     265                 mesh_origin = None,
    266266                 data_origin = None,
    267267                 precrop = False):
    268268
    269        
     269
    270270        """ Build interpolation matrix mapping from
    271271        function values at vertices to function values at data points
    272272
    273273        Inputs:
    274        
    275           vertex_coordinates: List of coordinate pairs [xi, eta] of 
     274
     275          vertex_coordinates: List of coordinate pairs [xi, eta] of
    276276          points constituting mesh (or a an m x 2 Numeric array)
    277        
     277
    278278          triangles: List of 3-tuples (or a Numeric array) of
    279279          integers representing indices of all vertices in the mesh.
    280280
    281           point_coordinates: List of coordinate pairs [x, y] of 
     281          point_coordinates: List of coordinate pairs [x, y] of
    282282          data points (or an nx2 Numeric array)
    283           If point_coordinates is absent, only smoothing matrix will 
    284           be built 
     283          If point_coordinates is absent, only smoothing matrix will
     284          be built
    285285
    286286          alpha: Smoothing parameter
     
    289289          UTM zone, easting and northing. If specified
    290290          point coordinates and vertex coordinates are assumed to be
    291           relative to their respective origins. 
    292          
     291          relative to their respective origins.
     292
    293293        """
    294294        from util import ensure_numeric
    295        
     295
    296296        #Convert input to Numeric arrays
    297297        triangles = ensure_numeric(triangles, Int)
     
    299299
    300300        #Build underlying mesh
    301         if verbose: print 'Building mesh' 
     301        if verbose: print 'Building mesh'
    302302        #self.mesh = General_mesh(vertex_coordinates, triangles,
    303303        #FIXME: Trying the normal mesh while testing precrop,
    304         #       The functionality of boundary_polygon is needed for that 
     304        #       The functionality of boundary_polygon is needed for that
    305305        self.mesh = Mesh(vertex_coordinates, triangles,
    306306                         origin = mesh_origin)
     
    322322                                        data_origin = data_origin,
    323323                                        precrop = precrop)
    324        
     324
    325325
    326326    def set_point_coordinates(self, point_coordinates,
     
    330330        """
    331331        self.build_coefficient_matrix_B(point_coordinates, data_origin)
    332        
     332
    333333    def build_coefficient_matrix_B(self, point_coordinates=None,
    334334                                   verbose = False, expand_search = True,
     
    337337                                   precrop = False):
    338338        """Build final coefficient matrix"""
    339        
     339
    340340
    341341        if self.alpha <> 0:
    342             if verbose: print 'Building smoothing matrix'         
     342            if verbose: print 'Building smoothing matrix'
    343343            self.build_smoothing_matrix_D()
    344        
     344
    345345        if point_coordinates is not None:
    346            
    347             if verbose: print 'Building interpolation matrix'         
     346
     347            if verbose: print 'Building interpolation matrix'
    348348            self.build_interpolation_matrix_A(point_coordinates,
    349349                                              verbose = verbose,
     
    361361            #Convert self.B matrix to CSR format for faster matrix vector
    362362            self.B = Sparse_CSR(self.B)
    363        
     363
    364364    def build_interpolation_matrix_A(self, point_coordinates,
    365365                                     verbose = False, expand_search = True,
     
    373373        This algorithm uses a quad tree data structure for fast binning of data points
    374374        origin is a 3-tuple consisting of UTM zone, easting and northing.
    375         If specified coordinates are assumed to be relative to this origin. 
     375        If specified coordinates are assumed to be relative to this origin.
    376376
    377377        This one will override any data_origin that may be specified in
    378         interpolation instance 
     378        interpolation instance
    379379
    380380        """
     
    386386            data_origin = self.data_origin #Use the one from
    387387                                           #interpolation instance
    388        
    389         #Convert input to Numeric arrays just in case. 
    390         point_coordinates = ensure_numeric(point_coordinates, Float)       
     388
     389        #Convert input to Numeric arrays just in case.
     390        point_coordinates = ensure_numeric(point_coordinates, Float)
    391391
    392392
     
    408408                    #
    409409                    #and similarly for eta
    410                    
     410
    411411                    x_offset = data_origin[1] - mesh_origin[1]
    412412                    y_offset = data_origin[2] - mesh_origin[2]
    413413                else: #Shift back to a zero origin
    414414                    x_offset = data_origin[1]
    415                     y_offset = data_origin[2]               
    416                    
     415                    y_offset = data_origin[2]
     416
    417417                point_coordinates[:,0] += x_offset
    418                 point_coordinates[:,1] += y_offset           
     418                point_coordinates[:,1] += y_offset
    419419            else:
    420420                if mesh_origin is not None:
    421421                    #Use mesh origin for data points
    422                     point_coordinates[:,0] -= mesh_origin[1] 
     422                    point_coordinates[:,0] -= mesh_origin[1]
    423423                    point_coordinates[:,1] -= mesh_origin[2]
    424424
     
    435435            P = self.mesh.get_boundary_polygon()
    436436
    437             if verbose: print 'Getting indices inside mesh boundary'           
     437            if verbose: print 'Getting indices inside mesh boundary'
    438438            indices = inside_polygon(point_coordinates, P, verbose = verbose)
    439            
     439
    440440            if verbose:
    441441                print 'Done'
     
    448448
    449449
    450        
    451         #Build n x m interpolation matrix       
     450
     451        #Build n x m interpolation matrix
    452452        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
    453453        n = point_coordinates.shape[0]     #Nbr of data points
     
    455455        if verbose: print 'Number of datapoints: %d' %n
    456456        if verbose: print 'Number of basis functions: %d' %m
    457        
     457
    458458        #FIXME (Ole): We should use CSR here since mat-mat mult is now OK.
    459459        #However, Sparse_CSR does not have the same methods as Sparse yet
     
    477477            candidate_vertices = root.search(x[0], x[1])
    478478            is_more_elements = True
    479            
     479
    480480            element_found, sigma0, sigma1, sigma2, k = \
    481481                self.search_triangles_of_vertices(candidate_vertices, x)
    482             while not element_found and is_more_elements and expand_search: 
    483                 #if verbose: print 'Expanding search' 
     482            while not element_found and is_more_elements and expand_search:
     483                #if verbose: print 'Expanding search'
    484484                candidate_vertices, branch = root.expand_search()
    485485                if branch == []:
     
    492492                    element_found, sigma0, sigma1, sigma2, k = \
    493493                      self.search_triangles_of_vertices(candidate_vertices, x)
    494                      
    495            
    496             #Update interpolation matrix A if necessary     
    497             if element_found is True:       
     494
     495
     496            #Update interpolation matrix A if necessary
     497            if element_found is True:
    498498                #Assign values to matrix A
    499499
     
    511511            else:
    512512                pass
    513                 #Ok if there is no triangle for datapoint 
     513                #Ok if there is no triangle for datapoint
    514514                #(as in brute force version)
    515515                #raise 'Could not find triangle for point', x
     
    529529            #For all vertices in same cell as point x
    530530            for v in candidate_vertices:
    531            
    532                 #for each triangle id (k) which has v as a vertex 
     531
     532                #for each triangle id (k) which has v as a vertex
    533533                for k, _ in self.mesh.vertexlist[v]:
    534                    
     534
    535535                    #Get the three vertex_points of candidate triangle
    536536                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
    537537                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
    538                     xi2 = self.mesh.get_vertex_coordinate(k, 2)     
     538                    xi2 = self.mesh.get_vertex_coordinate(k, 2)
    539539
    540540                    #print "PDSG - k", k
     
    544544                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
    545545                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
    546                    
     546
    547547                    #Get the three normals
    548548                    n0 = self.mesh.get_normal(k, 0)
    549549                    n1 = self.mesh.get_normal(k, 1)
    550                     n2 = self.mesh.get_normal(k, 2)               
    551 
    552                    
     550                    n2 = self.mesh.get_normal(k, 2)
     551
     552
    553553                    #Compute interpolation
    554554                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
     
    559559                    #print "PDSG - sigma1", sigma1
    560560                    #print "PDSG - sigma2", sigma2
    561                    
     561
    562562                    #FIXME: Maybe move out to test or something
    563563                    epsilon = 1.0e-6
    564564                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
    565                    
     565
    566566                    #Check that this triangle contains the data point
    567                    
     567
    568568                    #Sigmas can get negative within
    569569                    #machine precision on some machines (e.g nautilus)
    570                     #Hence the small eps                   
     570                    #Hence the small eps
    571571                    eps = 1.0e-15
    572572                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
     
    577577                    #Don't look for any other triangle
    578578                    break
    579             return element_found, sigma0, sigma1, sigma2, k     
    580                    
    581 
    582        
     579            return element_found, sigma0, sigma1, sigma2, k
     580
     581
     582
    583583    def build_interpolation_matrix_A_brute(self, point_coordinates):
    584584        """Build n x m interpolation matrix, where
     
    586586        m is the number of basis functions phi_k (one per vertex)
    587587
    588         This is the brute force which is too slow for large problems, 
     588        This is the brute force which is too slow for large problems,
    589589        but could be used for testing
    590590        """
    591591
    592592        from util import ensure_numeric
    593        
     593
    594594        #Convert input to Numeric arrays
    595         point_coordinates = ensure_numeric(point_coordinates, Float)         
    596        
    597         #Build n x m interpolation matrix       
     595        point_coordinates = ensure_numeric(point_coordinates, Float)
     596
     597        #Build n x m interpolation matrix
    598598        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
    599599        n = point_coordinates.shape[0]     #Nbr of data points
    600        
     600
    601601        self.A = Sparse(n,m)
    602602        self.AtA = Sparse(m,m)
     
    616616                xi0 = self.mesh.get_vertex_coordinate(k, 0)
    617617                xi1 = self.mesh.get_vertex_coordinate(k, 1)
    618                 xi2 = self.mesh.get_vertex_coordinate(k, 2)                 
     618                xi2 = self.mesh.get_vertex_coordinate(k, 2)
    619619
    620620                #Get the three normals
    621621                n0 = self.mesh.get_normal(k, 0)
    622622                n1 = self.mesh.get_normal(k, 1)
    623                 n2 = self.mesh.get_normal(k, 2)               
     623                n2 = self.mesh.get_normal(k, 2)
    624624
    625625                #Compute interpolation
     
    654654                            self.AtA[j,k] += sigmas[j]*sigmas[k]
    655655                k = k+1
    656        
    657 
    658        
     656
     657
     658
    659659    def get_A(self):
    660         return self.A.todense() 
     660        return self.A.todense()
    661661
    662662    def get_B(self):
    663663        return self.B.todense()
    664    
     664
    665665    def get_D(self):
    666666        return self.D.todense()
    667    
     667
    668668        #FIXME: Remember to re-introduce the 1/n factor in the
    669669        #interpolation term
    670        
     670
    671671    def build_smoothing_matrix_D(self):
    672672        """Build m x m smoothing matrix, where
     
    702702        self.D = Sparse(m,m)
    703703
    704         #For each triangle compute contributions to D = D1+D2       
     704        #For each triangle compute contributions to D = D1+D2
    705705        for i in range(len(self.mesh)):
    706706
     
    712712            v1 = self.mesh.triangles[i,1]
    713713            v2 = self.mesh.triangles[i,2]
    714            
     714
    715715            #Get the three vertex_points
    716716            xi0 = self.mesh.get_vertex_coordinate(i, 0)
    717717            xi1 = self.mesh.get_vertex_coordinate(i, 1)
    718             xi2 = self.mesh.get_vertex_coordinate(i, 2)                 
     718            xi2 = self.mesh.get_vertex_coordinate(i, 2)
    719719
    720720            #Compute gradients for each vertex
     
    726726
    727727            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
    728                               0, 0, 1)           
     728                              0, 0, 1)
    729729
    730730            #Compute diagonal contributions
    731731            self.D[v0,v0] += (a0*a0 + b0*b0)*area
    732732            self.D[v1,v1] += (a1*a1 + b1*b1)*area
    733             self.D[v2,v2] += (a2*a2 + b2*b2)*area           
     733            self.D[v2,v2] += (a2*a2 + b2*b2)*area
    734734
    735735            #Compute contributions for basis functions sharing edges
     
    744744            e20 = (a2*a0 + b2*b0)*area
    745745            self.D[v2,v0] += e20
    746             self.D[v0,v2] += e20             
    747 
    748            
     746            self.D[v0,v2] += e20
     747
     748
    749749    def fit(self, z):
    750750        """Fit a smooth surface to given 1d array of data points z.
    751751
    752752        The smooth surface is computed at each vertex in the underlying
    753         mesh using the formula given in the module doc string. 
     753        mesh using the formula given in the module doc string.
    754754
    755755        Pre Condition:
    756756          self.A, self.At and self.B have been initialised
    757          
     757
    758758        Inputs:
    759759          z: Single 1d vector or array of data at the point_coordinates.
    760760        """
    761        
     761
    762762        #Convert input to Numeric arrays
    763763        from util import ensure_numeric
    764764        z = ensure_numeric(z, Float)
    765        
     765
    766766        if len(z.shape) > 1 :
    767767            raise VectorShapeError, 'Can only deal with 1d data vector'
     
    769769        if self.point_indices is not None:
    770770            #Remove values for any points that were outside mesh
    771             z = take(z, self.point_indices) 
    772        
     771            z = take(z, self.point_indices)
     772
    773773        #Compute right hand side based on data
    774774        Atz = self.A.trans_mult(z)
    775775
    776        
     776
    777777        #Check sanity
    778778        n, m = self.A.shape
     
    788788
    789789        return conjugate_gradient(self.B, Atz, Atz,imax=2*len(Atz) )
    790         #FIXME: Should we store the result here for later use? (ON)       
    791 
    792            
     790        #FIXME: Should we store the result here for later use? (ON)
     791
     792
    793793    def fit_points(self, z, verbose=False):
    794794        """Like fit, but more robust when each point has two or more attributes
    795         FIXME (Ole): The name fit_points doesn't carry any meaning 
     795        FIXME (Ole): The name fit_points doesn't carry any meaning
    796796        for me. How about something like fit_multiple or fit_columns?
    797797        """
    798        
     798
    799799        try:
    800800            if verbose: print 'Solving penalised least_squares problem'
     
    805805            #Convert input to Numeric arrays
    806806            from util import ensure_numeric
    807             z = ensure_numeric(z, Float)           
    808            
    809             #Build n x m interpolation matrix       
     807            z = ensure_numeric(z, Float)
     808
     809            #Build n x m interpolation matrix
    810810            m = self.mesh.coordinates.shape[0] #Number of vertices
    811             n = z.shape[1]                     #Number of data points         
     811            n = z.shape[1]                     #Number of data points
    812812
    813813            f = zeros((m,n), Float) #Resulting columns
    814            
     814
    815815            for i in range(z.shape[1]):
    816816                f[:,i] = self.fit(z[:,i])
    817                
     817
    818818            return f
    819            
    820        
     819
     820
    821821    def interpolate(self, f):
    822         """Evaluate smooth surface f at data points implied in self.A. 
     822        """Evaluate smooth surface f at data points implied in self.A.
    823823
    824824        The mesh values representing a smooth surface are
    825825        assumed to be specified in f. This argument could,
    826826        for example have been obtained from the method self.fit()
    827        
     827
    828828        Pre Condition:
    829829          self.A has been initialised
    830        
     830
    831831        Inputs:
    832832          f: Vector or array of data at the mesh vertices.
    833833          If f is an array, interpolation will be done for each column as
    834834          per underlying matrix-matrix multiplication
    835          
     835
    836836        Output:
    837837          Interpolated values at data points implied in self.A
    838            
     838
    839839        """
    840840
    841841        return self.A * f
    842    
     842
    843843    def cull_outsiders(self, f):
    844844        pass
    845        
    846            
     845
     846
    847847#-------------------------------------------------------------
    848848if __name__ == "__main__":
     
    862862        point_file = sys.argv[2]
    863863        mesh_output_file = sys.argv[3]
    864        
     864
    865865        expand_search = False
    866866        if len(sys.argv) > 4:
    867867            if sys.argv[4][0] == "e" or sys.argv[4][0] == "E":
    868868                expand_search = True
    869             else:   
     869            else:
    870870                expand_search = False
    871                
     871
    872872        verbose = False
    873873        if len(sys.argv) > 5:
    874874            if sys.argv[5][0] == "n" or sys.argv[5][0] == "N":
    875875                verbose = False
    876             else:   
     876            else:
    877877                verbose = True
    878            
     878
    879879        if len(sys.argv) > 6:
    880880            alpha = sys.argv[6]
    881881        else:
    882882            alpha = DEFAULT_ALPHA
    883            
     883
    884884        t0 = time.time()
    885885        fit_to_mesh_file(mesh_file,
     
    889889                         verbose= verbose,
    890890                         expand_search = expand_search)
    891    
     891
    892892        print 'That took %.2f seconds' %(time.time()-t0)
    893        
     893
  • inundation/ga/storm_surge/pyvolution/mesh.py

    r1158 r1160  
    22
    33   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
    4    Geoscience Australia, 2004   
     4   Geoscience Australia, 2004
    55"""
    66
     
    1111
    1212    A triangular element is defined in terms of three vertex ids,
    13     ordered counter clock-wise, 
     13    ordered counter clock-wise,
    1414    each corresponding to a given coordinate set.
    1515    Vertices from different elements can point to the same
     
    1717
    1818    Coordinate sets are implemented as an N x 2 Numeric array containing
    19     x and y coordinates. 
    20    
     19    x and y coordinates.
     20
    2121
    2222    To instantiate:
     
    3838        c = [2.0,0.0]
    3939        e = [2.0, 2.0]
    40        
     40
    4141        points = [a, b, c, e]
    42         triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
     42        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
    4343        mesh = Mesh(points, triangles)
    44    
     44
    4545        #creates two triangles: bac and bce
    4646
     
    5656
    5757    #FIXME: Put in check for angles less than a set minimum
    58    
     58
    5959    def __init__(self, coordinates, triangles, boundary = None,
    6060                 tagged_elements = None, origin = None):
     
    7070
    7171        N = self.number_of_elements
    72        
     72
    7373        #Allocate space for geometric quantities
    7474        self.centroid_coordinates = zeros((N, 2), Float)
     
    7878        self.neighbours = zeros((N, 3), Int)
    7979        self.neighbour_edges = zeros((N, 3), Int)
    80         self.number_of_boundaries = zeros(N, Int) 
    81         self.surrogate_neighbours = zeros((N, 3), Int)       
    82        
     80        self.number_of_boundaries = zeros(N, Int)
     81        self.surrogate_neighbours = zeros((N, 3), Int)
     82
    8383        #Get x,y coordinates for all triangles and store
    8484        V = self.vertex_coordinates
    85        
     85
    8686        #Initialise each triangle
    8787        for i in range(N):
    8888            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
    89            
     89
    9090            x0 = V[i, 0]; y0 = V[i, 1]
    9191            x1 = V[i, 2]; y1 = V[i, 3]
    92             x2 = V[i, 4]; y2 = V[i, 5]           
     92            x2 = V[i, 4]; y2 = V[i, 5]
    9393
    9494            #Compute centroid
     
    9696            self.centroid_coordinates[i] = centroid
    9797
    98             #Midpoints       
     98            #Midpoints
    9999            m0 = array([(x1 + x2)/2, (y1 + y2)/2])
    100100            m1 = array([(x0 + x2)/2, (y0 + y2)/2])
     
    104104            #a triangle to the midpoint of the side of the triangle
    105105            #closest to the centroid
    106             d0 = sqrt(sum( (centroid-m0)**2 )) 
    107             d1 = sqrt(sum( (centroid-m1)**2 ))       
     106            d0 = sqrt(sum( (centroid-m0)**2 ))
     107            d1 = sqrt(sum( (centroid-m1)**2 ))
    108108            d2 = sqrt(sum( (centroid-m2)**2 ))
    109        
    110             self.radii[i] = min(d0, d1, d2)               
     109
     110            self.radii[i] = min(d0, d1, d2)
    111111
    112112            #Initialise Neighbours (-1 means that it is a boundary neighbour)
     
    118118
    119119
    120         #Build neighbour structure   
     120        #Build neighbour structure
    121121        self.build_neighbour_structure()
    122122
    123         #Build surrogate neighbour structure   
    124         self.build_surrogate_neighbour_structure()       
    125 
    126         #Build boundary dictionary mapping (id, edge) to symbolic tags       
     123        #Build surrogate neighbour structure
     124        self.build_surrogate_neighbour_structure()
     125
     126        #Build boundary dictionary mapping (id, edge) to symbolic tags
    127127        self.build_boundary_dictionary(boundary)
    128128
    129129        #Build tagged element  dictionary mapping (tag) to array of elements
    130130        self.build_tagged_elements_dictionary(tagged_elements)
    131        
     131
    132132        #Update boundary indices FIXME: OBSOLETE
    133         #self.build_boundary_structure()       
     133        #self.build_boundary_structure()
    134134
    135135        #FIXME check integrity?
    136        
     136
    137137    def __repr__(self):
    138138        return 'Mesh: %d triangles, %d elements, %d boundary segments'\
    139139               %(self.coordinates.shape[0], len(self), len(self.boundary))
    140140
    141            
     141
    142142    def build_neighbour_structure(self):
    143143        """Update all registered triangles to point to their neighbours.
    144        
     144
    145145        Also, keep a tally of the number of boundaries for each triangle
    146146
     
    149149          number_of_boundaries integer array is defined.
    150150        """
    151        
     151
    152152        #Step 1:
    153153        #Build dictionary mapping from segments (2-tuple of points)
    154154        #to left hand side edge (facing neighbouring triangle)
    155155
    156         N = self.number_of_elements       
     156        N = self.number_of_elements
    157157        neighbourdict = {}
    158158        for i in range(N):
     
    172172                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
    173173                    raise msg
    174                    
     174
    175175            neighbourdict[a,b] = (i, 2) #(id, edge)
    176176            neighbourdict[b,c] = (i, 0) #(id, edge)
     
    191191                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
    192192                self.number_of_boundaries[i] -= 1
    193                
     193
    194194            if neighbourdict.has_key((c,b)):
    195195                self.neighbours[i, 0] = neighbourdict[c,b][0]
    196196                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
    197                 self.number_of_boundaries[i] -= 1               
    198                
     197                self.number_of_boundaries[i] -= 1
     198
    199199            if neighbourdict.has_key((a,c)):
    200200                self.neighbours[i, 1] = neighbourdict[a,c][0]
    201201                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
    202                 self.number_of_boundaries[i] -= 1             
     202                self.number_of_boundaries[i] -= 1
    203203
    204204
     
    234234         keyed by volume id and edge:
    235235         { (id, edge): tag, ... }
    236          
     236
    237237         Postconditions:
    238238            self.boundary is defined.
     
    240240
    241241        from config import default_boundary_tag
    242        
     242
    243243        if boundary is None:
    244244            boundary = {}
     
    278278                            boundary[ (vol_id, edge_id) ] =\
    279279                                      default_boundary_tag
    280            
    281 
    282            
     280
     281
     282
    283283        self.boundary = boundary
    284            
     284
    285285
    286286    def build_tagged_elements_dictionary(self, tagged_elements = None):
     
    289289         keyed by tag:
    290290         { (tag): [e1, e2, e3..] }
    291          
     291
    292292         Postconditions:
    293293            self.element_tag is defined
    294294        """
    295295        from Numeric import array, Int
    296        
     296
    297297        if tagged_elements is None:
    298298            tagged_elements = {}
     
    301301            for tag in tagged_elements.keys():
    302302                tagged_elements[tag] = array(tagged_elements[tag]).astype(Int)
    303                
     303
    304304                msg = 'Not all elements exist. '
    305305                assert max(tagged_elements[tag]) < self.number_of_elements, msg
    306         #print "tagged_elements", tagged_elements               
     306        #print "tagged_elements", tagged_elements
    307307        self.tagged_elements = tagged_elements
    308            
     308
    309309    def build_boundary_structure(self):
    310         """Traverse boundary and 
     310        """Traverse boundary and
    311311        enumerate neighbour indices from -1 and
    312312        counting down.
     
    324324        """
    325325
    326         #FIXME: Now Obsolete - maybe use some comments from here in 
    327         #domain.set_boundary
     326        #FIXME: Now Obsolete - maybe use some comments from here in
     327        #domain.set_boundary
    328328
    329329        if self.boundary is None:
     
    332332            raise msg
    333333
    334        
     334
    335335        self.boundary_segments = self.boundary.keys()
    336336        self.boundary_segments.sort()
    337337
    338         index = -1       
     338        index = -1
    339339        for id, edge in self.boundary_segments:
    340340
     
    342342            #if self.neighbours[id, edge] > -1:
    343343            #    print 'Internal boundary'
    344            
     344
    345345            self.neighbours[id, edge] = index
    346346            index -= 1
    347347
    348    
     348
    349349    def get_boundary_tags(self):
    350350        """Return list of available boundary tags
    351351        """
    352        
     352
    353353        tags = {}
    354354        for v in self.boundary.values():
    355355            tags[v] = 1
    356356
    357         return tags.keys()   
    358            
     357        return tags.keys()
     358
    359359
    360360    def get_boundary_polygon(self):
     
    362362        """
    363363        from Numeric import allclose, sqrt, array, minimum, maximum
    364    
    365        
     364
     365
    366366
    367367        V = self.get_vertex_coordinates()
     
    377377
    378378        pmax = array( [max(self.coordinates[:,0]),
    379                        max(self.coordinates[:,1]) ] )       
    380        
     379                       max(self.coordinates[:,1]) ] )
     380
    381381        mindist = sqrt(sum( (pmax-pmin)**2 ))
    382382        for i, edge_id in self.boundary.keys():
     
    384384            if edge_id == 0: a = 1; b = 2
    385385            if edge_id == 1: a = 2; b = 0
    386             if edge_id == 2: a = 0; b = 1               
     386            if edge_id == 2: a = 0; b = 1
    387387
    388388            A = tuple(self.get_vertex_coordinate(i, a))
     
    393393            #a unique way of selecting
    394394            dist_A = sqrt(sum( (A-pmin)**2 ))
    395             dist_B = sqrt(sum( (B-pmin)**2 ))                         
     395            dist_B = sqrt(sum( (B-pmin)**2 ))
    396396
    397397            #Find minimal point
     
    415415            segments[A] = B
    416416
    417        
     417
    418418        #Start with smallest point and follow boundary (counter clock wise)
    419419        polygon = [p0]
     
    422422            polygon.append(p1)
    423423            p0 = p1
    424            
     424
    425425        return polygon
    426        
     426
    427427
    428428    def check_integrity(self):
     
    432432        Neighbour structure will be checked by class Mesh
    433433        """
    434        
     434
    435435        from config import epsilon
    436436        from math import pi
     
    438438
    439439        N = self.number_of_elements
    440        
     440
    441441        #Get x,y coordinates for all vertices for all triangles
    442442        V = self.get_vertex_coordinates()
     
    446446            x0 = V[i, 0]; y0 = V[i, 1]
    447447            x1 = V[i, 2]; y1 = V[i, 3]
    448             x2 = V[i, 4]; y2 = V[i, 5]           
     448            x2 = V[i, 4]; y2 = V[i, 5]
    449449
    450450            #Check that area hasn't been compromised
     
    454454                  %(x0,y0,x1,y1,x2,y2)
    455455            assert abs((area - ref)/area) < epsilon, msg
    456        
     456
    457457            #Check that points are arranged in counter clock-wise order
    458458            v0 = [x1-x0, y1-y0]
     
    462462            a0 = anglediff(v1, v0)
    463463            a1 = anglediff(v2, v1)
    464             a2 = anglediff(v0, v2)       
     464            a2 = anglediff(v0, v2)
    465465
    466466            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
     
    473473            normal0 = self.normals[i, 0:2]
    474474            normal1 = self.normals[i, 2:4]
    475             normal2 = self.normals[i, 4:6]           
     475            normal2 = self.normals[i, 4:6]
    476476
    477477            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
     
    480480
    481481        #Check that all vertices have been registered
    482         for v in self.vertexlist:
    483             msg = 'Some points do not belong to an element.\n'
    484             msg += 'Make sure all points appear as element vertices!'
    485             assert v is not None, msg
    486            
     482        for v_id, v in enumerate(self.vertexlist):
     483
     484            msg = 'Vertex %s does not belong to an element.'
     485            #assert v is not None, msg
     486            if v is None:
     487                print msg%v_id
     488
    487489        #Check integrity of neighbour structure
    488490        for i in range(N):
     
    495497                       (i, 1) in self.vertexlist[v] or\
    496498                       (i, 2) in self.vertexlist[v]
    497                
    498                
     499
     500
    499501
    500502            #Check neighbour structure
     
    505507                if neighbour_id >= 0:
    506508                    edge = self.neighbour_edges[i, k]
    507                     msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id) 
     509                    msg = 'Triangle %d has neighbour %d but it does not point back. \n' %(i,neighbour_id)
    508510                    msg += 'Only points to (%s)' %(self.neighbours[neighbour_id,:])
    509511                    assert self.neighbours[neighbour_id, edge] == i ,msg
     
    512514
    513515        #Check that all boundaries have
    514         # unique, consecutive, negative indices                   
     516        # unique, consecutive, negative indices
    515517
    516518        #L = len(self.boundary)
     
    531533        #
    532534        #See domain.set_boundary
    533        
    534            
     535
     536
    535537    def get_centroid_coordinates(self):
    536538        """Return all centroid coordinates.
    537539        Return all centroid coordinates for all triangles as an Nx2 array
    538         (ordered as x0, y0 for each triangle)       
     540        (ordered as x0, y0 for each triangle)
    539541        """
    540542        return self.centroid_coordinates
  • inundation/ga/storm_surge/pyvolution/test_cg_solve.py

    r1018 r1160  
    2828
    2929        assert allclose(x,xe)
     30
     31    def test_max_iter(self):
     32        """Small Sparse Matrix"""
     33
     34        A = [[2.0, -1.0, 0.0, 0.0 ],
     35             [-1.0, 2.0, -1.0, 0.0],
     36             [0.0, -1.0, 2.0, -1.0],
     37             [0.0,0.0, -1.0, 2.0]]
     38
     39        A = Sparse(A)
     40
     41        xe = [0.0, 1.0, 2.0, 3.0]
     42        b  = A*xe
     43        x =  [0.0, 0.0, 0.0, 0.0]
     44
     45        x = conjugate_gradient(A,b,x,iprint=0,imax=2)
    3046
    3147
     
    170186
    171187
    172 
  • inundation/ga/storm_surge/zeus/pyvolution.zpi

    r1156 r1160  
    3030    <Typedefs>Typedefs</Typedefs>
    3131    <Variables>Variables</Variables>
    32     <MakeFile></MakeFile>
     32    <MakeFile>makefile</MakeFile>
    3333    <SearchPath></SearchPath>
    3434    <RegexFile></RegexFile>
     
    4545    <ReleaseDirectory></ReleaseDirectory>
    4646    <ReleaseDebugger></ReleaseDebugger>
    47     <DebugMake></DebugMake>
    48     <DebugClean></DebugClean>
     47    <DebugMake>make $PF</DebugMake>
     48    <DebugClean>make clean</DebugClean>
    4949    <DebugCompile></DebugCompile>
    50     <DebugRebuild></DebugRebuild>
     50    <DebugRebuild>make $PF</DebugRebuild>
    5151    <DebugProjectOutput>Always</DebugProjectOutput>
    5252    <DebugCompilerOutput>Always</DebugCompilerOutput>
     
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
    75     <file>cg_solve.py</file>
    76     <file>compile.py</file>
    77     <file>config.py</file>
    78     <file>data_manager.py</file>
    79     <file>domain.py</file>
    80     <file>general_mesh.py</file>
    81     <file>interpolate_sww.py</file>
    82     <file>least_squares.py</file>
    83     <file>mesh.py</file>
    84     <file>mesh_factory.py</file>
    85     <file>netherlands.py</file>
    86     <file>pmesh2domain.py</file>
    87     <file>pressure_force.py</file>
    88     <file>quad.py</file>
    89     <file>quantity.py</file>
    90     <file>quantity_ext.c</file>
    91     <file>realtime_visualisation.py</file>
    92     <file>region.py</file>
    93     <file>shallow_water.py</file>
    94     <file>shallow_water_ext.c</file>
    95     <file>sparse.py</file>
    96     <file>sparse_ext.c</file>
    97     <file>test_advection.py</file>
    98     <file>test_all.py</file>
    99     <file>test_cg_solve.py</file>
    100     <file>test_data_manager.py</file>
    101     <file>test_domain.py</file>
    102     <file>test_general_mesh.py</file>
    103     <file>test_generic_boundary_conditions.py</file>
    104     <file>test_interpolate_sww.py</file>
    105     <file>test_least_squares.py</file>
    106     <file>test_mesh.py</file>
    107     <file>test_pmesh2domain.py</file>
    108     <file>test_quad.py</file>
    109     <file>test_quantity.py</file>
    110     <file>test_region.py</file>
    111     <file>test_shallow_water.py</file>
    112     <file>test_sparse.py</file>
    113     <file>test_util.py</file>
    114     <file>treenode.py</file>
    115     <file>util.py</file>
    116     <file>util_ext.c</file>
     75    <file>..\pyvolution\advection.py</file>
     76    <file>..\pyvolution\bed_w_eden_boundary.py</file>
     77    <file>..\pyvolution\bed_w_file_boundary.py</file>
     78    <file>..\pyvolution\cg_solve.py</file>
     79    <file>..\pyvolution\check_sww_tsh.py</file>
     80    <file>..\pyvolution\chris-dump.py</file>
     81    <file>..\pyvolution\combine_pts.py</file>
     82    <file>..\pyvolution\compile.py</file>
     83    <file>..\pyvolution\config.py</file>
     84    <file>..\pyvolution\cornell_room.py</file>
     85    <file>..\pyvolution\data_manager.py</file>
     86    <file>..\pyvolution\domain.py</file>
     87    <file>..\pyvolution\flatbed.py</file>
     88    <file>..\pyvolution\flatbed_compare.py</file>
     89    <file>..\pyvolution\general_mesh.py</file>
     90    <file>..\pyvolution\generic_boundary_conditions.py</file>
     91    <file>..\pyvolution\interpolate_sww.py</file>
     92    <file>..\pyvolution\island.py</file>
     93    <file>..\pyvolution\least_squares.py</file>
     94    <file>..\pyvolution\log.ini</file>
     95    <file>..\pyvolution\mesh.py</file>
     96    <file>..\pyvolution\mesh_factory.py</file>
     97    <file>..\pyvolution\most2nc.py</file>
     98    <file>..\pyvolution\netherlands.py</file>
     99    <file>..\pyvolution\netherlands-chris.py</file>
     100    <file>..\pyvolution\pmesh2domain.py</file>
     101    <file>..\pyvolution\pressure_force.py</file>
     102    <file>..\pyvolution\quad.py</file>
     103    <file>..\pyvolution\quantity.py</file>
     104    <file>..\pyvolution\quantity_ext.c</file>
     105    <file>..\pyvolution\realtime_visualisation.py</file>
     106    <file>..\pyvolution\region.py</file>
     107    <file>..\pyvolution\run_profile.py</file>
     108    <file>..\pyvolution\run_tsh_weir.py</file>
     109    <file>..\pyvolution\shallow_water.py</file>
     110    <file>..\pyvolution\shallow_water_ext.c</file>
     111    <file>..\pyvolution\show_balanced_limiters.py</file>
     112    <file>..\pyvolution\sparse.py</file>
     113    <file>..\pyvolution\sparse_ext.c</file>
     114    <file>..\pyvolution\test_advection.py</file>
     115    <file>..\pyvolution\test_all.py</file>
     116    <file>..\pyvolution\test_cg_solve.py</file>
     117    <file>..\pyvolution\test_combine_pts.py</file>
     118    <file>..\pyvolution\test_data_manager.py</file>
     119    <file>..\pyvolution\test_domain.py</file>
     120    <file>..\pyvolution\test_general_mesh.py</file>
     121    <file>..\pyvolution\test_generic_boundary_conditions.py</file>
     122    <file>..\pyvolution\test_interpolate_sww.py</file>
     123    <file>..\pyvolution\test_least_squares.py</file>
     124    <file>..\pyvolution\test_mesh.py</file>
     125    <file>..\pyvolution\test_most2nc.py</file>
     126    <file>..\pyvolution\test_pmesh2domain.py</file>
     127    <file>..\pyvolution\test_quad.py</file>
     128    <file>..\pyvolution\test_quantity.py</file>
     129    <file>..\pyvolution\test_region.py</file>
     130    <file>..\pyvolution\test_shallow_water.py</file>
     131    <file>..\pyvolution\test_sparse.py</file>
     132    <file>..\pyvolution\test_util.py</file>
     133    <file>..\pyvolution\treenode.py</file>
     134    <file>..\pyvolution\tsh2sww.py</file>
     135    <file>..\pyvolution\twolevels.py</file>
     136    <file>..\pyvolution\util.py</file>
     137    <file>..\pyvolution\util_ext.c</file>
     138    <file>..\pyvolution\util_ext.h</file>
     139    <file>..\pyvolution\wind_example_variable.py</file>
    117140    <folder name="Header Files" />
    118141    <folder name="Resource Files" />
Note: See TracChangeset for help on using the changeset viewer.