Changeset 6410


Ignore:
Timestamp:
Feb 25, 2009, 9:37:22 AM (14 years ago)
Author:
rwilson
Message:

numpy changes.

Location:
branches/numpy/anuga
Files:
47 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

    r6304 r6410  
    377377            #indices = range(M)
    378378
    379         return num.take(self.triangles, indices)
     379        return num.take(self.triangles, indices, axis=0)
    380380   
    381381

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6304 r6410  
    479479                                                     indices, verbose=verbose,
    480480                                                     use_cache=use_cache)
    481 #            if (isinstance(numeric, float) or isinstance(numeric, int)
    482 #                or isinstance(numeric, long)):
    483481            else:
    484482                try:
    485483                    numeric = float(numeric)
    486484                except ValueError:
    487                     msg = ("Illegal type for variable 'numeric': "
    488                            "%s, shape=%s\n%s"
    489                            % (type(numeric), numeric.shape, str(numeric)))
    490                     msg += ('\ntype(numeric)==FloatType is %s'
    491                             % str(type(numeric)==FloatType))
    492                     msg += ('\nisinstance(numeric, float)=%s'
    493                             % str(isinstance(numeric, float)))
    494                     msg += ('\nisinstance(numeric, num.ndarray)=%s'
    495                             % str(isinstance(numeric, num.ndarray)))
     485                    msg = ("Illegal type for variable 'numeric': %s"
     486                           % type(numeric))
    496487                    raise Exception, msg
    497488                self.set_values_from_constant(numeric, location,
    498489                                              indices, verbose)
    499         elif quantity is not None:
    500             if type(numeric) in [FloatType, IntType, LongType]:
    501                 self.set_values_from_constant(numeric, location,
    502                                               indices, verbose)
    503             else:
    504                 msg = ("Illegal type for variable 'numeric': %s, shape=%s\n%s"
    505                        % (type(numeric), numeric.shape, str(numeric)))
    506                 msg += ('\ntype(numeric)==FloatType is %s'
    507                         % str(type(numeric)==FloatType))
    508                 msg += ('\nisinstance(numeric, float)=%s'
    509                         % str(isinstance(numeric, float)))
    510                 msg += ('\nisinstance(numeric, num.ndarray)=%s'
    511                         % str(isinstance(numeric, num.ndarray)))
    512                 raise Exception, msg
    513490        elif quantity is not None:
    514491            self.set_values_from_quantity(quantity, location, indices, verbose)
    … …  
    745722                indices = range(len(self))
    746723
    747             V = num.take(self.domain.get_centroid_coordinates(), indices)
     724            V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0)
    748725            x = V[:,0]; y = V[:,1]
    749726            if use_cache is True:
    … …  
    12131190            if (indices ==  None):
    12141191                indices = range(len(self))
    1215             return num.take(self.centroid_values, indices)
     1192            return num.take(self.centroid_values, indices, axis=0)
    12161193        elif location == 'edges':
    12171194            if (indices ==  None):
    12181195                indices = range(len(self))
    1219             return num.take(self.edge_values, indices)
     1196            return num.take(self.edge_values, indices, axis=0)
    12201197        elif location == 'unique vertices':
    12211198            if (indices ==  None):
    … …  
    12431220            if (indices is None):
    12441221                indices = range(len(self))
    1245             return num.take(self.vertex_values, indices)
     1222            return num.take(self.vertex_values, indices, axis=0)
    12461223
    12471224    ##

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r6304 r6410  
    10301030        }
    10311031       
     1032        // check that numpy array objects arrays are C contiguous memory
     1033        CHECK_C_CONTIG(vertex_value_indices);
     1034        CHECK_C_CONTIG(number_of_triangles_per_node);
     1035        CHECK_C_CONTIG(vertex_values);
     1036        CHECK_C_CONTIG(A);
     1037
    10321038        N = vertex_value_indices -> dimensions[0];
    10331039        // printf("Got parameters, N=%d\n", N);

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6304 r6410  
    764764       
    765765        #print domain.quantities['friction'].get_values()
     766        msg = ("domain.quantities['friction'].get_values()=\n%s\n"
     767               'should equal\n'
     768               '[[ 0.09,  0.09,  0.09],\n'
     769               ' [ 0.09,  0.09,  0.09],\n'
     770               ' [ 0.07,  0.07,  0.07],\n'
     771               ' [ 0.07,  0.07,  0.07],\n'
     772               ' [ 1.0,  1.0,  1.0],\n'
     773               ' [ 1.0,  1.0,  1.0]]'
     774               % str(domain.quantities['friction'].get_values()))
    766775        assert num.allclose(domain.quantities['friction'].get_values(),
    767776                            [[ 0.09,  0.09,  0.09],
    … …  
    770779                             [ 0.07,  0.07,  0.07],
    771780                             [ 1.0,  1.0,  1.0],
    772                              [ 1.0,  1.0,  1.0]])
     781                             [ 1.0,  1.0,  1.0]]), msg
    773782       
    774783        domain.set_region([set_bottom_friction, set_top_friction])
    … …  
    793802
    794803#-------------------------------------------------------------
     804
    795805if __name__ == "__main__":
    796806    suite = unittest.makeSuite(Test_Domain,'test')

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ermapper.py

    r6304 r6410  
    183183       
    184184#-------------------------------------------------------------
     185
    185186if __name__ == "__main__":
    186187    suite = unittest.makeSuite(Test_ERMapper,'test')
    187188    runner = unittest.TextTestRunner()
    188189    runner.run(suite)
    189 
    190 
    191 
    192 

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6304 r6410  
    6363        domain = General_mesh(nodes, triangles,
    6464                       geo_reference = geo)
     65
    6566        verts = domain.get_vertex_coordinates(triangle_id=0)       
    66         self.assert_(num.allclose(num.array([b,a,c]), verts))
     67        msg = ("num.array([b,a,c])=\n%s\nshould be the same as 'verts'=\n%s"
     68               % (str(num.array([b,a,c])), str(verts)))
     69        #self.assert_(num.allclose(num.array([b,a,c]), verts))
     70        assert num.allclose(num.array([b,a,c]), verts), msg
     71
    6772        verts = domain.get_vertex_coordinates(triangle_id=0)       
    6873        self.assert_(num.allclose(num.array([b,a,c]), verts))
    … …  
    117122        domain = General_mesh(nodes, triangles)
    118123
    119         value = [7]
    120         assert num.allclose(domain.get_triangles(), triangles)
     124#        value = [7]
     125        msg = ("domain.get_triangles()=\n%s\nshould be the same as "
     126               "'triangles'=\n%s"
     127               % (str(domain.get_triangles()), str(triangles)))
     128        assert num.allclose(domain.get_triangles(), triangles), msg
     129        msg = ("domain.get_triangles([0,4])=\n%s\nshould be the same as "
     130               "'[triangles[0], triangles[4]]' which is\n%s"
     131               % (str(domain.get_triangles([0,4])),
     132                  str([triangles[0], triangles[4]])))
    121133        assert num.allclose(domain.get_triangles([0,4]),
    122                         [triangles[0], triangles[4]])
     134                            [triangles[0], triangles[4]]), msg
    123135       
    124136
    … …  
    278290        node = domain.get_node(2)       
    279291        #self.assertEqual(c, node)
    280         self.failUnless(num.alltrue(c == node))
     292        msg = ('\nc=%s\nmode=%s' % (str(c), str(node)))
     293        self.failUnless(num.alltrue(c == node), msg)
    281294       
    282295        node = domain.get_node(2, absolute=True)     
    … …  
    328341
    329342#-------------------------------------------------------------
     343
    330344if __name__ == "__main__":
    331345    suite = unittest.makeSuite(Test_General_Mesh,'test')
    332     #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')   
    333346    runner = unittest.TextTestRunner()
    334347    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r6304 r6410  
    270270
    271271
    272 
    273 
    274 
    275272#-------------------------------------------------------------
     273
    276274if __name__ == "__main__":
    277275    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test')
    278276    runner = unittest.TextTestRunner()
    279277    runner.run(suite)
    280 
    281 
    282 

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ghost.py

    r6304 r6410  
    4545
    4646
     47#-------------------------------------------------------------
    4748
    48 
    49 #-------------------------------------------------------------
    5049if __name__ == "__main__":
    5150    suite = unittest.makeSuite(Test_Domain,'test')

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r6304 r6410  
    18391839
    18401840
    1841 
    1842        
    18431841#-------------------------------------------------------------
     1842
    18441843if __name__ == "__main__":
    1845     #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
    1846     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
    1847     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
    18481844    suite = unittest.makeSuite(Test_Mesh,'test')
    18491845    runner = unittest.TextTestRunner()
    18501846    runner.run(suite)
    1851 
    1852 
    1853 
    1854 

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

    r6304 r6410  
    255255
    256256#-------------------------------------------------------------
     257
    257258if __name__ == "__main__":
    258259    suite = unittest.makeSuite(Test_pmesh2domain, 'test')
    259260    runner = unittest.TextTestRunner()
    260261    runner.run(suite)
    261 
    262 
    263 
    264 

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py

    r6304 r6410  
    504504        quantity.set_values(0.0)
    505505        quantity.set_values(3.14, polygon=polygon)
     506        msg = ('quantity.vertex_values=\n%s\nshould be close to\n'
     507               '[[0,0,0],\n'
     508               ' [3.14,3.14,3.14],\n'
     509               ' [3.14,3.14,3.14],\n'
     510               ' [0,0,0]]' % str(quantity.vertex_values))
    506511        assert num.allclose(quantity.vertex_values,
    507512                            [[0,0,0],
    508513                             [3.14,3.14,3.14],
    509514                             [3.14,3.14,3.14],                         
    510                              [0,0,0]])               
     515                             [0,0,0]]), msg
    511516
    512517
    … …  
    25082513
    25092514#-------------------------------------------------------------
     2515
    25102516if __name__ == "__main__":
    25112517    suite = unittest.makeSuite(Test_Quantity, 'test')   
    2512     #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
    25132518    runner = unittest.TextTestRunner()
    25142519    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py

    r6304 r6410  
    189189                         
    190190    def test_unique_vertices_average_loc_vert(self):
    191         """
    192         get values based on triangle lists.
    193         """
    194         from mesh_factory import rectangular
    195         from shallow_water import Domain
    196 
    197         #Create basic mesh
    198         points, vertices, boundary = rectangular(1, 3)
    199         #Create shallow water domain
    200         domain = Domain(points, vertices, boundary)
    201         domain.build_tagged_elements_dictionary({'bottom':[0,1],
    202                                                  'top':[4,5],
    203                                                  'not_bottom':[2,3,4,5]})
     191        """Get values based on triangle lists."""
     192
     193        from mesh_factory import rectangular
     194        from shallow_water import Domain
     195
     196        #Create basic mesh
     197        points, vertices, boundary = rectangular(1, 3)
     198
     199        #Create shallow water domain
     200        domain = Domain(points, vertices, boundary)
     201        domain.build_tagged_elements_dictionary({'bottom': [0, 1],
     202                                                 'top': [4, 5],
     203                                                 'not_bottom': [2, 3, 4, 5]})
    204204
    205205        #Set friction
    206206        domain.set_quantity('friction', add_x_y)
    207         av_bottom = 2.0/3.0
     207        av_bottom = 2.0 / 3.0
    208208        add = 60.0
    209209        calc_frict = av_bottom + add
    210210        domain.set_region(Add_value_to_region('bottom', 'friction', add,
    211                           initial_quantity='friction',
    212                            #location='unique vertices',
    213                            location='vertices',
    214                            average=True
    215                           ))
    216 
    217         #print domain.quantities['friction'].get_values()
     211                                              initial_quantity='friction',
     212                                              location='vertices',
     213                                              average=True))
     214
    218215        frict_points = domain.quantities['friction'].get_values()
     216
    219217        expected = [calc_frict, calc_frict, calc_frict]
    220         msg = ('frict_points[0]=%s\nexpected=%s' % (str(frict_points[0]),
    221                                                    str(expected)))
     218        msg = ('frict_points[0]=%s\nexpected=%s'
     219               % (str(frict_points[0]), str(expected)))
    222220        assert num.allclose(frict_points[0], expected), msg
    223         msg = ('frict_points[1]=%s\nexpected=%s' % (str(frict_points[1]),
    224                                                     str(expected)))
     221
     222        msg = ('frict_points[1]=%s\nexpected=%s'
     223               % (str(frict_points[1]), str(expected)))
    225224        assert num.allclose(frict_points[1], expected), msg
    226225 
    … …  
    264263                         
    265264#-------------------------------------------------------------
     265
    266266if __name__ == "__main__":
    267     suite = unittest.makeSuite(Test_Region, 'test')   
     267    #suite = unittest.makeSuite(Test_Region, 'test')   
     268    suite = unittest.makeSuite(Test_Region, 'test_unique_vertices_average_loc_vert')   
    268269    runner = unittest.TextTestRunner()
    269270    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_util.py

    r6304 r6410  
    183183
    184184        last_time_index = len(time)-1 #Last last_time_index
    185         d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1))
    186         d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
    187         d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
     185        d_stage = num.reshape(num.take(stage[last_time_index, :],
     186                                       [0,5,10,15],
     187                                       axis=0),
     188                              (4,1))
     189        d_uh = num.reshape(num.take(xmomentum[last_time_index, :],
     190                                    [0,5,10,15],
     191                                   axis=0),
     192                           (4,1))
     193        d_vh = num.reshape(num.take(ymomentum[last_time_index, :],
     194                                    [0,5,10,15],
     195                                   axis=0),
     196                           (4,1))
    188197        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    189198
    … …  
    195204
    196205        #And the midpoints are found now
    197         Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15])
    198         Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15])
     206        Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0)
     207        Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0)
    199208
    200209        diag = num.concatenate( (Dx, Dy), axis=1)
    … …  
    219228
    220229        timestep = 0 #First timestep
    221         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    222         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    223         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     230        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     231        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     232        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    224233        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    225234
    … …  
    241250
    242251        timestep = 33
    243         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    244         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    245         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     252        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     253        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     254        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    246255        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    247256
    … …  
    262271
    263272        timestep = 15
    264         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    265         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    266         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     273        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     274        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     275        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    267276        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    268277
    … …  
    275284        #
    276285        timestep = 16
    277         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    278         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    279         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     286        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     287        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     288        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    280289        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    281290
    … …  
    386395
    387396        last_time_index = len(time)-1 #Last last_time_index     
    388         d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1))
    389         d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
    390         d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
     397        d_stage = num.reshape(num.take(stage[last_time_index, :],
     398                                       [0,5,10,15],
     399                                       axis=0),
     400                              (4,1))
     401        d_uh = num.reshape(num.take(xmomentum[last_time_index, :],
     402                                    [0,5,10,15],
     403                                   axis=0),
     404                           (4,1))
     405        d_vh = num.reshape(num.take(ymomentum[last_time_index, :],
     406                                    [0,5,10,15],
     407                                    axis=0),
     408                           (4,1))
    391409        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    392410
    … …  
    398416
    399417        #And the midpoints are found now
    400         Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15])
    401         Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15])
     418        Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0)
     419        Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0)
    402420
    403421        diag = num.concatenate( (Dx, Dy), axis=1)
    … …  
    415433
    416434        t = time[last_time_index]                         
    417         q = f(t, point_id=0); assert num.allclose(r0, q)
    418         q = f(t, point_id=1); assert num.allclose(r1, q)
    419         q = f(t, point_id=2); assert num.allclose(r2, q)
     435
     436        q = f(t, point_id=0)
     437        msg = '\nr0=%s\nq=%s' % (str(r0), str(q))
     438        assert num.allclose(r0, q), msg
     439
     440        q = f(t, point_id=1)
     441        msg = '\nr1=%s\nq=%s' % (str(r1), str(q))
     442        assert num.allclose(r1, q), msg
     443
     444        q = f(t, point_id=2)
     445        msg = '\nr2=%s\nq=%s' % (str(r2), str(q))
     446        assert num.allclose(r2, q), msg
    420447
    421448
    … …  
    424451
    425452        timestep = 0 #First timestep
    426         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    427         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    428         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     453        d_stage = num.reshape(num.take(stage[timestep, :],
     454                                       [0,5,10,15],
     455                                       axis=0),
     456                              (4,1))
     457        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     458                                    [0,5,10,15],
     459                                    axis=0),
     460                           (4,1))
     461        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     462                                    [0,5,10,15],
     463                                    axis=0),
     464                           (4,1))
    429465        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    430466
    … …  
    446482
    447483        timestep = 33
    448         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    449         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    450         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     484        d_stage = num.reshape(num.take(stage[timestep, :],
     485                                       [0,5,10,15],
     486                                       axis=0),
     487                              (4,1))
     488        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     489                                    [0,5,10,15],
     490                                    axis=0),
     491                           (4,1))
     492        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     493                                    [0,5,10,15],
     494                                    axis=0),
     495                           (4,1))
    451496        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    452497
    … …  
    467512
    468513        timestep = 15
    469         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    470         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    471         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     514        d_stage = num.reshape(num.take(stage[timestep, :],
     515                                       [0,5,10,15],
     516                                       axis=0),
     517                              (4,1))
     518        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     519                                    [0,5,10,15],
     520                                    axis=0),
     521                           (4,1))
     522        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     523                                    [0,5,10,15],
     524                                    axis=0),
     525                           (4,1))
    472526        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    473527
    … …  
    480534        #
    481535        timestep = 16
    482         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    483         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    484         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     536        d_stage = num.reshape(num.take(stage[timestep, :],
     537                                       [0,5,10,15],
     538                                       axis=0),
     539                              (4,1))
     540        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     541                                    [0,5,10,15],
     542                                    axis=0),
     543                           (4,1))
     544        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     545                                    [0,5,10,15],
     546                                    axis=0),
     547                           (4,1))
    485548        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    486549
    … …  
    19131976        if 314 < angle < 316: v=1
    19141977        assert v==1
    1915 
    1916  
    1917 
    1918 
    19191978       
    19201979
    19211980#-------------------------------------------------------------
     1981
    19221982if __name__ == "__main__":
    19231983    suite = unittest.makeSuite(Test_Util,'test')
    1924 #    suite = unittest.makeSuite(Test_Util,'test_sww2csv_gauges')
    19251984#    runner = unittest.TextTestRunner(verbosity=2)
    19261985    runner = unittest.TextTestRunner(verbosity=1)
    19271986    runner.run(suite)
    1928 
    1929 
    1930 
    1931 

  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6304 r6410  
    476476        if boundary_polygon is not None:
    477477            #removes sts points that do not lie on boundary
    478             quantities[name] = num.take(quantities[name], gauge_id, 1)
     478            #quantities[name] = num.take(quantities[name], gauge_id, 1)  #was#
     479            quantities[name] = num.take(quantities[name], gauge_id, axis=1)
    479480           
    480481    # Close sww, tms or sts netcdf file         
    … …  
    17981799        # Remove the loners from verts
    17991800        # Could've used X=compress(less(loners,N),loners)
    1800         # verts=num.take(verts,X)  to Remove the loners from verts
     1801        # verts=num.take(verts,X,axis=0)  to Remove the loners from verts
    18011802        # but I think it would use more memory
    18021803        new_i = lone_start      # point at first loner - 'shuffle down' target

  • branches/numpy/anuga/alpha_shape/test_alpha_shape.py

    r6304 r6410  
    394394                  (46, 45)]
    395395        assert num.allclose(answer, result)   
     396
    396397#-------------------------------------------------------------
     398
    397399if __name__ == "__main__":
    398     #print "starting tests \n"
    399400    suite = unittest.makeSuite(TestCase,'test')
    400401    runner = unittest.TextTestRunner(verbosity=1)
    401402    runner.run(suite)
    402     #print "finished tests \n"
    403    
    404 
    405    
    406    
    407 
    408 
    409 

  • branches/numpy/anuga/caching/test_caching.py

    r6304 r6410  
    891891#-------------------------------------------------------------
    892892if __name__ == "__main__":
    893     #suite = unittest.makeSuite(Test_Caching, 'test_caching_of_callable_objects')
    894893    suite = unittest.makeSuite(Test_Caching, 'test')
    895894    runner = unittest.TextTestRunner()

  • branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

    r6360 r6410  
    505505       
    506506#-------------------------------------------------------------
     507
    507508if __name__ == "__main__":
    508 
    509509    suite = unittest.makeSuite(geo_referenceTestCase,'test')
    510     #suite = unittest.makeSuite(geo_referenceTestCase,'test_get_absolute')
    511510    runner = unittest.TextTestRunner() #verbosity=2)
    512511    runner.run(suite)

  • branches/numpy/anuga/coordinate_transforms/test_lat_long_UTM_conversion.py

    r6304 r6410  
    121121        assert num.allclose(lat,  lat_calced)
    122122        assert num.allclose(lon, long_calced)
     123
    123124#-------------------------------------------------------------
     125
    124126if __name__ == "__main__":
    125 
    126     #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
    127127    mysuite = unittest.makeSuite(TestCase,'test')
    128128    runner = unittest.TextTestRunner()

  • branches/numpy/anuga/coordinate_transforms/test_redfearn.py

    r6304 r6410  
    411411           
    412412#-------------------------------------------------------------
     413
    413414if __name__ == "__main__":
    414 
    415     #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
    416     #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')
    417415    mysuite = unittest.makeSuite(TestCase,'test')
    418416    runner = unittest.TextTestRunner()

  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6304 r6410  
    6868                max_vertices_per_cell=None,
    6969                start_blocking_len=500000,
    70                 use_cache=False,             
    71                 verbose=False):                 
     70                use_cache=False,
     71                verbose=False):
    7272    """Interpolate vertex_values to interpolation points.
    73    
     73
    7474    Inputs (mandatory):
    7575
    76    
     76
    7777    vertex_coordinates: List of coordinate pairs [xi, eta] of
    78                         points constituting a mesh 
     78                        points constituting a mesh
    7979                        (or an m x 2 numeric array or
    8080                        a geospatial object)
    … …  
    8383
    8484    triangles: List of 3-tuples (or a numeric array) of
    85                integers representing indices of all vertices 
     85               integers representing indices of all vertices
    8686               in the mesh.
    87                
     87
    8888    vertex_values: Vector or array of data at the mesh vertices.
    8989                   If array, interpolation will be done for each column as
    9090                   per underlying matrix-matrix multiplication
    91              
     91
    9292    interpolation_points: Interpolate mesh data to these positions.
    9393                          List of coordinate pairs [x, y] of
    94                           data points or an nx2 numeric array or a
     94                          data points or an nx2 numeric array or a
    9595                          Geospatial_data object
    96                
     96
    9797    Inputs (optional)
    98                
     98
    9999    mesh_origin: A geo_reference object or 3-tuples consisting of
    100100                 UTM zone, easting and northing.
    … …  
    108108                           object and a mesh origin, since geospatial has its
    109109                           own mesh origin.
    110                            
     110
    111111    start_blocking_len: If the # of points is more or greater than this,
    112                         start blocking 
    113                        
     112                        start blocking
     113
    114114    use_cache: True or False
    115                            
     115
    116116
    117117    Output:
    118    
     118
    119119    Interpolated values at specified point_coordinates
    120120
    121     Note: This function is a simple shortcut for case where 
     121    Note: This function is a simple shortcut for case where
    122122    interpolation matrix is unnecessary
    123     Note: This function does not take blocking into account, 
     123    Note: This function does not take blocking into account,
    124124    but allows caching.
    125    
     125
    126126    """
    127127
    128128    # FIXME(Ole): Probably obsolete since I is precomputed and
    129129    #             interpolate_block caches
    130        
     130
    131131    from anuga.caching import cache
    132132
    133133    # Create interpolation object with matrix
    134     args = (ensure_numeric(vertex_coordinates, num.float), 
     134    args = (ensure_numeric(vertex_coordinates, num.float),
    135135            ensure_numeric(triangles))
    136136    kwargs = {'mesh_origin': mesh_origin,
    137137              'max_vertices_per_cell': max_vertices_per_cell,
    138138              'verbose': verbose}
    139    
     139
    140140    if use_cache is True:
    141141        if sys.platform != 'win32':
    … …  
    149149    else:
    150150        I = apply(Interpolate, args, kwargs)
    151    
     151
    152152    # Call interpolate method with interpolation points
    153153    result = I.interpolate_block(vertex_values, interpolation_points,
    154154                                 use_cache=use_cache,
    155155                                 verbose=verbose)
    156                            
     156
    157157    return result
    158158
    159159
    160160##
    161 # @brief 
     161# @brief
    162162class Interpolate (FitInterpolate):
    163163
    … …  
    181181        Inputs:
    182182          vertex_coordinates: List of coordinate pairs [xi, eta] of
    183               points constituting a mesh (or an m x 2 numeric array or
     183              points constituting a mesh (or an m x 2 numeric array or
    184184              a geospatial object)
    185185              Points may appear multiple times
    … …  
    202202
    203203        # FIXME (Ole): Need an input check
    204        
     204
    205205        FitInterpolate.__init__(self,
    206206                                vertex_coordinates=vertex_coordinates,
    … …  
    209209                                verbose=verbose,
    210210                                max_vertices_per_cell=max_vertices_per_cell)
    211                                
     211
    212212        # Initialise variables
    213213        self._A_can_be_reused = False  # FIXME (Ole): Probably obsolete
    … …  
    215215        self.interpolation_matrices = {} # Store precomputed matrices
    216216
    217    
     217
    218218
    219219    ##
    … …  
    243243          point_coordinates: Interpolate mesh data to these positions.
    244244              List of coordinate pairs [x, y] of
    245               data points or an nx2 numeric array or a Geospatial_data object
    246              
    247               If point_coordinates is absent, the points inputted last time
     245              data points or an nx2 numeric array or a Geospatial_data object
     246
     247              If point_coordinates is absent, the points inputted last time
    248248              this method was called are used, if possible.
    249249
    250250          start_blocking_len: If the # of points is more or greater than this,
    251               start blocking 
    252 
    253         Output:
    254           Interpolated values at inputted points (z).
     251              start blocking
     252
     253        Output:
     254          Interpolated values at inputted points (z).
    255255        """
    256256
    … …  
    262262        # This has now been addressed through an attempt in interpolate_block
    263263
    264         if verbose: print 'Build intepolation object' 
     264        if verbose: print 'Build intepolation object'
    265265        if isinstance(point_coordinates, Geospatial_data):
    266266            point_coordinates = point_coordinates.get_data_points(absolute=True)
    … …  
    280280                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
    281281                raise Exception(msg)
    282            
     282
    283283        if point_coordinates is not None:
    284284            self._point_coordinates = point_coordinates
    … …  
    293293                start = 0
    294294                # creating a dummy array to concatenate to.
    295                
     295
    296296                f = ensure_numeric(f, num.float)
    297297                if len(f.shape) > 1:
    … …  
    299299                else:
    300300                    z = num.zeros((0,), num.int)        #array default#
    301                    
     301
    302302                for end in range(start_blocking_len,
    303303                                 len(point_coordinates),
    … …  
    313313                z = num.concatenate((z, t))
    314314        return z
    315    
     315
    316316
    317317    ##
    … …  
    322322    # @param verbose True if this function is verbose.
    323323    # @return ??
    324     def interpolate_block(self, f, point_coordinates, 
     324    def interpolate_block(self, f, point_coordinates,
    325325                          use_cache=False, verbose=False):
    326326        """
    … …  
    329329
    330330        Return the point data, z.
    331        
     331
    332332        See interpolate for doc info.
    333333        """
    334        
     334
    335335        # FIXME (Ole): I reckon we should change the interface so that
    336         # the user can specify the interpolation matrix instead of the 
     336        # the user can specify the interpolation matrix instead of the
    337337        # interpolation points to save time.
    338338
    … …  
    342342        # Convert lists to numeric arrays if necessary
    343343        point_coordinates = ensure_numeric(point_coordinates, num.float)
    344         f = ensure_numeric(f, num.float)               
     344        f = ensure_numeric(f, num.float)
    345345
    346346        from anuga.caching import myhash
    347         import sys 
    348          
     347        import sys
     348
    349349        if use_cache is True:
    350350            if sys.platform != 'win32':
    351351                # FIXME (Ole): (Why doesn't this work on windoze?)
    352352                # Still absolutely fails on Win 24 Oct 2008
    353            
     353
    354354                X = cache(self._build_interpolation_matrix_A,
    355355                          args=(point_coordinates,),
    356                           kwargs={'verbose': verbose},                       
    357                           verbose=verbose)       
     356                          kwargs={'verbose': verbose},
     357                          verbose=verbose)
    358358            else:
    359359                # FIXME
    … …  
    361361                # This will work on Linux as well if we want to use it there.
    362362                key = myhash(point_coordinates)
    363        
    364                 reuse_A = False 
    365            
     363
     364                reuse_A = False
     365
    366366                if self.interpolation_matrices.has_key(key):
    367                     X, stored_points = self.interpolation_matrices[key] 
     367                    X, stored_points = self.interpolation_matrices[key]
    368368                    if num.alltrue(stored_points == point_coordinates):
    369                         reuse_A = True          # Reuse interpolation matrix
    370                
     369                        reuse_A = True                # Reuse interpolation matrix
     370
    371371                if reuse_A is False:
    372372                    X = self._build_interpolation_matrix_A(point_coordinates,
    … …  
    375375        else:
    376376            X = self._build_interpolation_matrix_A(point_coordinates,
    377                                                    verbose=verbose)           
    378 
    379         # Unpack result                                       
     377                                                   verbose=verbose)
     378
     379        # Unpack result
    380380        self._A, self.inside_poly_indices, self.outside_poly_indices = X
    381381
    … …  
    389389        msg += ' I got %d points and %d matrix rows.' \
    390390               % (point_coordinates.shape[0], self._A.shape[0])
    391         assert point_coordinates.shape[0] == self._A.shape[0], msg       
     391        assert point_coordinates.shape[0] == self._A.shape[0], msg
    392392
    393393        msg = 'The number of columns in matrix A must be the same as the '
    394394        msg += 'number of mesh vertices.'
    395395        msg += ' I got %d vertices and %d matrix columns.' \
    396                % (f.shape[0], self._A.shape[1])       
     396               % (f.shape[0], self._A.shape[1])
    397397        assert self._A.shape[1] == f.shape[0], msg
    398398
    … …  
    409409        """
    410410        Return the point data, z.
    411        
     411
    412412        Precondition: The _A matrix has been created
    413413        """
    … …  
    416416
    417417        # Taking into account points outside the mesh.
    418         for i in self.outside_poly_indices: 
     418        for i in self.outside_poly_indices:
    419419            z[i] = NAN
    420420        return z
    … …  
    456456                                   self.mesh.get_boundary_polygon(),
    457457                                   closed=True, verbose=verbose)
    458        
     458
    459459        # Build n x m interpolation matrix
    460460        if verbose and len(outside_poly_indices) > 0:
    461461            print '\n WARNING: Points outside mesh boundary. \n'
    462            
     462
    463463        # Since you can block, throw a warning, not an error.
    464464        if verbose and 0 == len(inside_poly_indices):
    465465            print '\n WARNING: No points within the mesh! \n'
    466            
     466
    467467        m = self.mesh.number_of_nodes  # Nbr of basis functions (1/vertex)
    468468        n = point_coordinates.shape[0] # Nbr of data points
    … …  
    474474
    475475        n = len(inside_poly_indices)
    476        
     476
    477477        # Compute matrix elements for points inside the mesh
    478478        if verbose: print 'Building interpolation matrix from %d points' %n
    … …  
    485485            element_found, sigma0, sigma1, sigma2, k = \
    486486                           search_tree_of_vertices(self.root, self.mesh, x)
    487            
    488             # Update interpolation matrix A if necessary
     487
     488            # Update interpolation matrix A if necessary
    489489            if element_found is True:
    490490                # Assign values to matrix A
    … …  
    499499                    A[i, j] = sigmas[j]
    500500            else:
    501                 msg = 'Could not find triangle for point', x 
     501                msg = 'Could not find triangle for point', x
    502502                raise Exception(msg)
    503503        return A, inside_poly_indices, outside_poly_indices
    504504
    505505
    506        
    507        
    508        
    509        
     506
     507
     508
     509
    510510##
    511511# @brief ??
    … …  
    527527            List of coordinate pairs [x, y] of
    528528            data points or an nx2 numeric array or a Geospatial_data object
    529              
     529
    530530    No test for this yet.
    531531    Note, this has no time the input data has no time dimension.  Which is
    532532    different from most of the data we interpolate, eg sww info.
    533  
     533
    534534    Output:
    535535        Interpolated values at inputted points.
    … …  
    537537
    538538    interp = Interpolate(vertices,
    539                          triangles, 
     539                         triangles,
    540540                         max_vertices_per_cell=max_points_per_cell,
    541541                         mesh_origin=mesh_origin)
    542            
     542
    543543    calc = interp.interpolate(vertex_attributes,
    544544                              points,
    … …  
    554554# @param velocity_y_file Name of the output y velocity  file.
    555555# @param stage_file Name of the output stage file.
    556 # @param froude_file 
     556# @param froude_file
    557557# @param time_thinning Time thinning step to use.
    558558# @param verbose True if this function is to be verbose.
    … …  
    604604                                 time_thinning=time_thinning,
    605605                                 use_cache=use_cache)
    606    
     606
    607607    depth_writer = writer(file(depth_file, "wb"))
    608608    velocity_x_writer = writer(file(velocity_x_file, "wb"))
    … …  
    620620    velocity_y_writer.writerow(heading)
    621621    if stage_file is not None:
    622         stage_writer.writerow(heading) 
     622        stage_writer.writerow(heading)
    623623    if froude_file is not None:
    624         froude_writer.writerow(heading)     
    625    
     624        froude_writer.writerow(heading)
     625
    626626    for time in callable_sww.get_time():
    627627        depths = [time]
    628628        velocity_xs = [time]
    629629        velocity_ys = [time]
    630         if stage_file is not None: 
    631             stages = [time] 
    632         if froude_file is not None: 
    633             froudes = [time] 
     630        if stage_file is not None:
     631            stages = [time]
     632        if froude_file is not None:
     633            froudes = [time]
    634634        for point_i, point in enumerate(points):
    635635            quantities = callable_sww(time,point_i)
    636            
     636
    637637            w = quantities[0]
    638638            z = quantities[1]
    639639            momentum_x = quantities[2]
    640640            momentum_y = quantities[3]
    641             depth = w - z 
    642              
     641            depth = w - z
     642
    643643            if w == NAN or z == NAN or momentum_x == NAN:
    644644                velocity_x = NAN
    … …  
    677677
    678678        if stage_file is not None:
    679             stage_writer.writerow(stages)   
     679            stage_writer.writerow(stages)
    680680        if froude_file is not None:
    681             froude_writer.writerow(froudes)           
     681            froude_writer.writerow(froudes)
    682682
    683683
    684684##
    685 # @brief 
     685# @brief
    686686class Interpolation_function:
    687687    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
    … …  
    695695    Mandatory input
    696696        time:                 px1 array of monotonously increasing times (float)
    697         quantities:           Dictionary of arrays or 1 array (float) 
     697        quantities:           Dictionary of arrays or 1 array (float)
    698698                              The arrays must either have dimensions pxm or mx1.
    699699                              The resulting function will be time dependent in
    700700                              the former case while it will be constant with
    701701                              respect to time in the latter case.
    702        
     702
    703703    Optional input:
    704704        quantity_names:       List of keys into the quantities dictionary for
    … …  
    706706        vertex_coordinates:   mx2 array of coordinates (float)
    707707        triangles:            nx3 array of indices into vertex_coordinates (Int)
    708         interpolation_points: Nx2 array of coordinates to be interpolated to 
     708        interpolation_points: Nx2 array of coordinates to be interpolated to
    709709        verbose:              Level of reporting
    710710
    … …  
    742742                 time,
    743743                 quantities,
    744                  quantity_names=None, 
     744                 quantity_names=None,
    745745                 vertex_coordinates=None,
    746746                 triangles=None,
    … …  
    762762            print 'Interpolation_function: input checks'
    763763
    764         # Check temporal info
     764        # Check temporal info
    765765        time = ensure_numeric(time)
    766766        if not num.alltrue(time[1:] - time[:-1] >= 0):
    … …  
    791791            if triangles is not None:
    792792                triangles = ensure_numeric(triangles)
    793             self.spatial = True         
     793            self.spatial = True
    794794
    795795        if verbose is True:
    796796            print 'Interpolation_function: thinning by %d' % time_thinning
    797797
    798            
     798
    799799        # Thin timesteps if needed
    800800        # Note array() is used to make the thinned arrays contiguous in memory
    801         self.time = num.array(time[::time_thinning])         
     801        self.time = num.array(time[::time_thinning])
    802802        for name in quantity_names:
    803803            if len(quantities[name].shape) == 2:
    … …  
    806806        if verbose is True:
    807807            print 'Interpolation_function: precomputing'
    808            
     808
    809809        # Save for use with statistics
    810810        self.quantities_range = {}
    … …  
    812812            q = quantities[name][:].flatten()
    813813            self.quantities_range[name] = [min(q), max(q)]
    814        
    815         self.quantity_names = quantity_names       
    816         self.vertex_coordinates = vertex_coordinates 
     814
     815        self.quantity_names = quantity_names
     816        self.vertex_coordinates = vertex_coordinates
    817817        self.interpolation_points = interpolation_points
    818818
    … …  
    825825            #if self.spatial is False:
    826826            #    raise 'Triangles and vertex_coordinates must be specified'
    827             # 
     827            #
    828828            try:
    829                 self.interpolation_points = \
     829                self.interpolation_points = \
    830830                    interpolation_points = ensure_numeric(interpolation_points)
    831831            except:
    832                 msg = 'Interpolation points must be an N x 2 numeric array ' \
     832                msg = 'Interpolation points must be an N x 2 numeric array ' \
    833833                      'or a list of points\n'
    834                 msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
     834                msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
    835835                raise msg
    836836
    … …  
    839839                # mesh boundary as defined by triangles and vertex_coordinates.
    840840                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
    841                 from anuga.utilities.polygon import outside_polygon           
     841                from anuga.utilities.polygon import outside_polygon
    842842
    843843                # Create temporary mesh object from mesh info passed
    844                 # into this function. 
     844                # into this function.
    845845                mesh = Mesh(vertex_coordinates, triangles)
    846846                mesh_boundary_polygon = mesh.get_boundary_polygon()
    … …  
    868868                            # FIXME (Ole): Why only Windoze?
    869869                            from anuga.utilities.polygon import plot_polygons
    870                             #out_interp_pts = num.take(interpolation_points,[indices])
    871870                            title = ('Interpolation points fall '
    872871                                     'outside specified mesh')
    … …  
    886885                        print msg
    887886                    #raise Exception(msg)
    888                    
     887
    889888            elif triangles is None and vertex_coordinates is not None:    #jj
    890889                #Dealing with sts file
    … …  
    911910            m = len(self.interpolation_points)
    912911            p = len(self.time)
    913            
    914             for name in quantity_names:
     912
     913            for name in quantity_names:
    915914                self.precomputed_values[name] = num.zeros((p, m), num.float)
    916915
    … …  
    918917                print 'Build interpolator'
    919918
    920                
     919
    921920            # Build interpolator
    922921            if triangles is not None and vertex_coordinates is not None:
    923                 if verbose:               
     922                if verbose:
    924923                    msg = 'Building interpolation matrix from source mesh '
    925924                    msg += '(%d vertices, %d triangles)' \
    … …  
    927926                              triangles.shape[0])
    928927                    print msg
    929                
     928
    930929                # This one is no longer needed for STS files
    931930                interpol = Interpolate(vertex_coordinates,
    932931                                       triangles,
    933                                        verbose=verbose)               
    934                
     932                                       verbose=verbose)
     933
    935934            elif triangles is None and vertex_coordinates is not None:
    936935                if verbose:
    … …  
    943942                print 'Interpolating (%d interpolation points, %d timesteps).' \
    944943                      %(self.interpolation_points.shape[0], self.time.shape[0]),
    945            
     944
    946945                if time_thinning > 1:
    947946                    print 'Timesteps were thinned by a factor of %d' \
    … …  
    950949                    print
    951950
    952             for i, t in enumerate(self.time):
     951            for i, t in enumerate(self.time):
    953952                # Interpolate quantities at this timestep
    954953                if verbose and i%((p+10)/10) == 0:
    955954                    print '  time step %d of %d' %(i, p)
    956                    
     955
    957956                for name in quantity_names:
    958957                    if len(quantities[name].shape) == 2:
    … …  
    963962                    if verbose and i%((p+10)/10) == 0:
    964963                        print '    quantity %s, size=%d' %(name, len(Q))
    965                        
    966                     # Interpolate 
     964
     965                    # Interpolate
    967966                    if triangles is not None and vertex_coordinates is not None:
    968967                        result = interpol.interpolate(Q,
    … …  
    976975                                                      interpolation_points=\
    977976                                                          self.interpolation_points)
    978                        
     977
    979978                    #assert len(result), len(interpolation_points)
    980979                    self.precomputed_values[name][i, :] = result
    … …  
    10031002    def __call__(self, t, point_id=None, x=None, y=None):
    10041003        """Evaluate f(t) or f(t, point_id)
    1005        
    1006         Inputs:
    1007           t:        time - Model time. Must lie within existing timesteps
    1008           point_id: index of one of the preprocessed points.
    1009 
    1010           If spatial info is present and all of point_id
    1011           are None an exception is raised 
    1012                    
     1004
     1005        Inputs:
     1006          t:        time - Model time. Must lie within existing timesteps
     1007          point_id: index of one of the preprocessed points.
     1008
     1009          If spatial info is present and all of point_id
     1010          are None an exception is raised
     1011
    10131012          If no spatial info is present, point_id arguments are ignored
    10141013          making f a function of time only.
    10151014
    1016           FIXME: f(t, x, y) x, y could overrided location, point_id ignored 
    1017           FIXME: point_id could also be a slice
    1018           FIXME: What if x and y are vectors?
     1015          FIXME: f(t, x, y) x, y could overrided location, point_id ignored
     1016          FIXME: point_id could also be a slice
     1017          FIXME: What if x and y are vectors?
    10191018          FIXME: What about f(x,y) without t?
    10201019        """
    10211020
    10221021        from math import pi, cos, sin, sqrt
    1023         from anuga.abstract_2d_finite_volumes.util import mean       
     1022        from anuga.abstract_2d_finite_volumes.util import mean
    10241023
    10251024        if self.spatial is True:
    … …  
    10571056        # Compute interpolated values
    10581057        q = num.zeros(len(self.quantity_names), num.float)
    1059         for i, name in enumerate(self.quantity_names):
     1058        for i, name in enumerate(self.quantity_names):
    10601059            Q = self.precomputed_values[name]
    10611060
    10621061            if self.spatial is False:
    1063                 # If there is no spatial info               
     1062                # If there is no spatial info
    10641063                assert len(Q.shape) == 1
    10651064
    … …  
    10761075                        Q1 = Q[self.index+1, point_id]
    10771076
    1078             # Linear temporal interpolation   
     1077            # Linear temporal interpolation
    10791078            if ratio > 0:
    10801079                if Q0 == NAN and Q1 == NAN:
    … …  
    10921091            # Replicate q according to x and y
    10931092            # This is e.g used for Wind_stress
    1094             if x is None or y is None: 
     1093            if x is None or y is None:
    10951094                return q
    10961095            else:
    … …  
    11061105                    for col in q:
    11071106                        res.append(col*num.ones(N, num.float))
    1108                        
     1107
    11091108                return res
    11101109
    … …  
    11221121        """Output statistics about interpolation_function
    11231122        """
    1124        
     1123
    11251124        vertex_coordinates = self.vertex_coordinates
    1126         interpolation_points = self.interpolation_points               
     1125        interpolation_points = self.interpolation_points
    11271126        quantity_names = self.quantity_names
    11281127        #quantities = self.quantities
    1129         precomputed_values = self.precomputed_values                 
    1130                
     1128        precomputed_values = self.precomputed_values
     1129
    11311130        x = vertex_coordinates[:,0]
    1132         y = vertex_coordinates[:,1]               
     1131        y = vertex_coordinates[:,1]
    11331132
    11341133        str =  '------------------------------------------------\n'
    … …  
    11441143        for name in quantity_names:
    11451144            minq, maxq = self.quantities_range[name]
    1146             str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
     1145            str += '    %s in [%f, %f]\n' %(name, minq, maxq)
    11471146            #q = quantities[name][:].flatten()
    11481147            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
    11491148
    1150         if interpolation_points is not None:   
     1149        if interpolation_points is not None:
    11511150            str += '  Interpolation points (xi, eta):'\
    11521151                   ' number of points == %d\n' %interpolation_points.shape[0]
    … …  
    11561155                                             max(interpolation_points[:,1]))
    11571156            str += '  Interpolated quantities (over all timesteps):\n'
    1158        
     1157
    11591158            for name in quantity_names:
    11601159                q = precomputed_values[name][:].flatten()
    … …  
    11851184    print "x",x
    11861185    print "y",y
    1187    
     1186
    11881187    print "time", time
    11891188    print "quantities", quantities
    … …  
    11951194    interp = Interpolation_interface(time,
    11961195                                     quantities,
    1197                                      quantity_names=quantity_names, 
     1196                                     quantity_names=quantity_names,
    11981197                                     vertex_coordinates=vertex_coordinates,
    11991198                                     triangles=volumes,
    … …  
    12081207    """
    12091208    obsolete - Nothing should be calling this
    1210    
     1209
    12111210    Read in an sww file.
    1212    
     1211
    12131212    Input;
    12141213    file_name - the sww file
    1215    
     1214
    12161215    Output;
    12171216    x - Vector of x values
    … …  
    12261225    msg = 'Function read_sww in interpolat.py is obsolete'
    12271226    raise Exception, msg
    1228    
     1227
    12291228    #FIXME Have this reader as part of data_manager?
    1230    
    1231     from Scientific.IO.NetCDF import NetCDFFile     
     1229
     1230    from Scientific.IO.NetCDF import NetCDFFile
    12321231    import tempfile
    12331232    import sys
    12341233    import os
    1235    
     1234
    12361235    #Check contents
    12371236    #Get NetCDF
    1238    
     1237
    12391238    # see if the file is there.  Throw a QUIET IO error if it isn't
    12401239    # I don't think I could implement the above
    1241    
     1240
    12421241    #throws prints to screen if file not present
    12431242    junk = tempfile.mktemp(".txt")
    … …  
    12451244    stdout = sys.stdout
    12461245    sys.stdout = fd
    1247     fid = NetCDFFile(file_name, netcdf_mode_r) 
     1246    fid = NetCDFFile(file_name, netcdf_mode_r)
    12481247    sys.stdout = stdout
    12491248    fd.close()
    12501249    #clean up
    1251     os.remove(junk)     
    1252      
     1250    os.remove(junk)
     1251
    12531252    # Get the variables
    12541253    x = fid.variables['x'][:]
    12551254    y = fid.variables['y'][:]
    1256     volumes = fid.variables['volumes'][:] 
     1255    volumes = fid.variables['volumes'][:]
    12571256    time = fid.variables['time'][:]
    12581257
    … …  
    12661265    for name in keys:
    12671266        quantities[name] = fid.variables[name][:]
    1268            
     1267
    12691268    fid.close()
    12701269    return x, y, volumes, time, quantities

  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6360 r6410  
    216216        else:
    217217            self.data_points = ensure_numeric(data_points)
    218             return
    219 
    220             print 'self.data_points=%s' % str(self.data_points)
    221             print 'self.data_points.shape=%s' % str(self.data_points.shape)
    222218            if not (0,) == self.data_points.shape:
    223219                assert len(self.data_points.shape) == 2
    … …  
    551547
    552548        attributes = {}
    553         if file_name[-4:]== ".pts":
     549        if file_name[-4:] == ".pts":
    554550            try:
    555551                data_points, attributes, geo_reference = \
    … …  
    558554                msg = 'Could not open file %s ' % file_name
    559555                raise IOError, msg
    560         elif file_name[-4:]== ".txt" or file_name[-4:]== ".csv":
     556        elif file_name[-4:] == ".txt" or file_name[-4:]== ".csv":
    561557            try:
    562558                data_points, attributes, geo_reference = \
    … …  
    588584    def export_points_file(self, file_name, absolute=True,
    589585                           as_lat_long=False, isSouthHemisphere=True):
    590 
    591586        """write a points file as a text (.csv) or binary (.pts) file
    592587
    … …  
    659654        # FIXME: add the geo_reference to this
    660655        points = self.get_data_points()
    661         sampled_points = num.take(points, indices)
     656        sampled_points = num.take(points, indices, axis=0)
    662657
    663658        attributes = self.get_all_attributes()
    … …  
    666661        if attributes is not None:
    667662            for key, att in attributes.items():
    668                 sampled_attributes[key] = num.take(att, indices)
     663                sampled_attributes[key] = num.take(att, indices, axis=0)
    669664
    670665        return Geospatial_data(sampled_points, sampled_attributes)
    … …  
    695690        """
    696691
    697         i=0
     692        i = 0
    698693        self_size = len(self)
    699694        random_list = []
    … …  
    704699        if verbose: print "make unique random number list and get indices"
    705700
    706         total=num.array(range(self_size))
     701        total = num.array(range(self_size), num.int)    #array default#
    707702        total_list = total.tolist()
    708703
    … …  
    734729        if verbose: print "make random number list and get indices"
    735730
    736         j=0
    737         k=1
     731        j = 0
     732        k = 1
    738733        remainder_list = total_list[:]
    739734
    … …  
    10771072    """
    10781073
    1079     line = file_pointer.readline().strip()
     1074    line = file_pointer.readline()
    10801075    header = clean_line(line, delimiter)
    10811076
    … …  
    12521247    # Create new file
    12531248    outfile.institution = 'Geoscience Australia'
    1254     outfile.description = 'NetCDF format for compact and portable storage ' \
    1255                           'of spatial point data'
     1249    outfile.description = ('NetCDF format for compact and portable storage '
     1250                           'of spatial point data')
    12561251
    12571252    # Dimension definitions
    … …  
    15781573
    15791574
    1580     attribute_smoothed='elevation'
     1575    attribute_smoothed = 'elevation'
    15811576
    15821577    if mesh_file is None:
    … …  
    16371632    data = num.array([], dtype=num.float)
    16381633
    1639     data=num.resize(data, (len(points), 3+len(alphas)))
     1634    data = num.resize(data, (len(points), 3+len(alphas)))
    16401635
    16411636    # gets relative point from sample
    … …  
    16451640    data[:,2] = elevation_sample
    16461641
    1647     normal_cov=num.array(num.zeros([len(alphas), 2]), dtype=num.float)
     1642    normal_cov = num.array(num.zeros([len(alphas), 2]), dtype=num.float)
    16481643
    16491644    if verbose: print 'Setup computational domains with different alphas'
    … …  
    16841679
    16851680    normal_cov0 = normal_cov[:,0]
    1686     normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0))
     1681    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)
    16871682
    16881683    if plot_name is not None:
    … …  
    19031898
    19041899    normal_cov0 = normal_cov[:,0]
    1905     normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0))
     1900    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)
    19061901
    19071902    if plot_name is not None:

  • branches/numpy/anuga/geospatial_data/test_geospatial_data.py

    r6360 r6410  
    207207        new_geo_ref = Geo_reference(56, x_p, y_p)
    208208        G.set_geo_reference(new_geo_ref)
    209         msg = ('points_ab=%s, but\nG.get_data_points(absolute=True)=%s'
     209        msg = ('points_ab=\n%s\nbut G.get_data_points(absolute=True)=\n%s'
    210210               % (str(points_ab), str(G.get_data_points(absolute=True))))
    211211        assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg
    … …  
    342342
    343343        assert num.allclose(P, num.concatenate( (P1,P2) ))
    344         msg = 'P=%s' % str(P)
     344        msg = ('P=\n%s\nshould be close to\n'
     345               '[[2.0, 4.1], [4.0, 7.3],\n'
     346               ' [5.1, 9.1], [6.1, 6.3]]'
     347               % str(P))
    345348        assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3],
    346349                                [5.1, 9.1], [6.1, 6.3]]), msg
    … …  
    369372
    370373        P1 = G1.get_data_points(absolute=False)
    371         msg = ('P1=%s, but\npoints1 - <...>=%s'
     374        msg = ('P1=\n%s\nbut points1 - <...>=\n%s'
    372375               % (str(P1),
    373376                  str(points1 - [geo_ref1.get_xllcorner(),
    … …  
    439442        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
    440443        P = G.get_data_points(absolute=True)
    441         msg = 'P=%s' % str(P)
     444        msg = 'P=\n%s' % str(P)
    442445        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg
    443446
    … …  
    10221025        os.remove(fileName)
    10231026
    1024         msg = ('results.get_data_points()=%s, but\nG.get_data_points(True)=%s'
     1027        msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s'
    10251028               % (str(results.get_data_points()), str(G.get_data_points(True))))
    10261029        assert num.allclose(results.get_data_points(),

  • branches/numpy/anuga/load_mesh/loadASCII.py

    r6360 r6410  
    626626        num.array(mesh['vertex_attribute_titles'], num.character)
    627627    mesh['segments'] = num.array(mesh['segments'], IntType)
    628     print ("Before num.array(), mesh['segment_tags']=%s"
    629            % str(mesh['segment_tags']))
    630628    mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character)
    631     print ("After num.array(), mesh['segment_tags']=%s"
    632            % str(mesh['segment_tags']))
    633     print "mesh['segment_tags'].shape=%s" % str(mesh['segment_tags'].shape)
    634629    mesh['triangles'] = num.array(mesh['triangles'], IntType)
    635630    mesh['triangle_tags'] = num.array(mesh['triangle_tags'])
    … …  
    11061101def take_points(dict, indices_to_keep):
    11071102    dict = point_atts2array(dict)
    1108     dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep)
     1103    dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0)
    11091104
    11101105    for key in dict['attributelist'].keys():
    11111106        dict['attributelist'][key] = num.take(dict['attributelist'][key],
    1112                                               indices_to_keep)
     1107                                              indices_to_keep, axis=0)
    11131108
    11141109    return dict

  • branches/numpy/anuga/load_mesh/test_loadASCII.py

    r6360 r6410  
    472472if __name__ == '__main__':
    473473    suite = unittest.makeSuite(loadASCIITestCase,'test')
    474     #suite = unittest.makeSuite(loadASCIITestCase,'test_read_write_msh_file6')
    475474    runner = unittest.TextTestRunner() #verbosity=2)
    476475    runner.run(suite)

  • branches/numpy/anuga/mesh_engine/mesh_engine_c_layer.c

    r6304 r6410  
    6060#include "Python.h"
    6161
     62#include "util_ext.h"
     63
     64
    6265//#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH
    6366
    … …  
    119122    return NULL;
    120123  }
     124
     125  // check that numpy array objects arrays are C contiguous memory
     126  CHECK_C_CONTIG(pointlist);
     127  CHECK_C_CONTIG(seglist);
     128  CHECK_C_CONTIG(holelist);
     129  CHECK_C_CONTIG(regionlist);
     130  CHECK_C_CONTIG(pointattributelist);
     131  CHECK_C_CONTIG(segmarkerlist);
    121132 
    122133  /* Initialize  points */

  • branches/numpy/anuga/mesh_engine/test_all.py

    r4898 r6410  
    7676    return unittest.TestSuite(map(load, modules))
    7777
     78
    7879if __name__ == '__main__':   
    7980    suite = regressionTest()

  • branches/numpy/anuga/mesh_engine/test_generate_mesh.py

    r6304 r6410  
    463463                                     correct.flat),
    464464                        'Failed')
    465      
     465   
     466
    466467if __name__ == "__main__":
    467 
    468468    suite = unittest.makeSuite(triangTestCase,'test')
    469469    runner = unittest.TextTestRunner()  #verbosity=2)

  • branches/numpy/anuga/pmesh/test_mesh.py

    r6360 r6410  
    23162316if __name__ == "__main__":
    23172317    suite = unittest.makeSuite(meshTestCase,'test')
    2318     #suite = unittest.makeSuite(meshTestCase,'test_import_mesh')
    23192318    runner = unittest.TextTestRunner() #verbosity=2)
    23202319    runner.run(suite)

  • branches/numpy/anuga/shallow_water/data_manager.py

    r6304 r6410  
    26692669    # Interpolate using quantity values
    26702670    if verbose: print 'Interpolating'
    2671     interpolated_values = interp.interpolate(q, data_points).flatten
     2671    interpolated_values = interp.interpolate(q, data_points).flatten()    #????#
    26722672
    26732673    if verbose:
    … …  
    73087308            # FIXME (Ole): Make a generic polygon input check in polygon.py
    73097309            # and call it here
    7310 
    7311             points = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
    7312 
     7310            points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
     7311                                                            y[:,num.newaxis]),
     7312                                                           axis=1))
    73137313            point_indices = inside_polygon(points, polygon)
    73147314
    73157315            # Restrict quantities to polygon
    7316             elevation = num.take(elevation, point_indices)
     7316            elevation = num.take(elevation, point_indices, axis=0)
    73177317            stage = num.take(stage, point_indices, axis=1)
    73187318
    73197319            # Get info for location of maximal runup
    7320             points_in_polygon = num.take(points, point_indices)
     7320            points_in_polygon = num.take(points, point_indices, axis=0)
     7321
    73217322            x = points_in_polygon[:,0]
    73227323            y = points_in_polygon[:,1]
    … …  
    73767377            else:
    73777378                # Find maximum elevation among wet nodes
    7378                 wet_elevation = num.take(elevation, wet_nodes)
     7379                wet_elevation = num.take(elevation, wet_nodes, axis=0)
    73797380                runup_index = num.argmax(wet_elevation)
    73807381                runup = max(wet_elevation)
    … …  
    73857386
    73867387                # Record location
    7387                 wet_x = num.take(x, wet_nodes)
    7388                 wet_y = num.take(y, wet_nodes)
     7388                wet_x = num.take(x, wet_nodes, axis=0)
     7389                wet_y = num.take(y, wet_nodes, axis=0)
    73897390                maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]
    73907391

  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r6304 r6410  
    1 """Finite-volume computations of the shallow water wave equation.
     1"""
     2Finite-volume computations of the shallow water wave equation.
    23
    34Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
    … …  
    5354    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
    5455
    55     Hydrodynamic modelling of coastal inundation. 
     56    Hydrodynamic modelling of coastal inundation.
    5657    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
    5758    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
    … …  
    7172    $Author$
    7273    $Date$
    73 
    7474"""
    7575
    … …  
    7979# $LastChangedRevision$
    8080# $LastChangedBy$
     81
    8182
    8283import numpy as num
    … …  
    109110from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    110111
    111 from anuga.fit_interpolate.interpolate import Modeltime_too_late, Modeltime_too_early
    112 
    113 from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
     112from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
     113                                              Modeltime_too_early
     114
     115from anuga.utilities.polygon import inside_polygon, polygon_area, \
     116                                    is_inside_polygon
    114117
    115118from types import IntType, FloatType
    … …  
    117120
    118121
    119 
    120 #
     122################################################################################
    121123# Shallow water domain
    122 #---------------------
     124################################################################################
     125
     126##
     127# @brief Class for a shallow water domain.
    123128class Domain(Generic_Domain):
    124129
    125130    conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
    126131    other_quantities = ['elevation', 'friction']
    127    
     132
     133    ##
     134    # @brief Instantiate a shallow water domain.
     135    # @param coordinates
     136    # @param vertices
     137    # @param boundary
     138    # @param tagged_elements
     139    # @param geo_reference
     140    # @param use_inscribed_circle
     141    # @param mesh_filename
     142    # @param use_cache
     143    # @param verbose
     144    # @param full_send_dict
     145    # @param ghost_recv_dict
     146    # @param processor
     147    # @param numproc
     148    # @param number_of_full_nodes
     149    # @param number_of_full_triangles
    128150    def __init__(self,
    129151                 coordinates=None,
    … …  
    142164                 number_of_full_nodes=None,
    143165                 number_of_full_triangles=None):
    144 
    145166
    146167        other_quantities = ['elevation', 'friction']
    … …  
    162183                                numproc,
    163184                                number_of_full_nodes=number_of_full_nodes,
    164                                 number_of_full_triangles=number_of_full_triangles)
    165 
     185                                number_of_full_triangles=number_of_full_triangles)
    166186
    167187        self.set_minimum_allowed_height(minimum_allowed_height)
    168        
     188
    169189        self.maximum_allowed_speed = maximum_allowed_speed
    170190        self.g = g
    171         self.beta_w      = beta_w
    172         self.beta_w_dry  = beta_w_dry
    173         self.beta_uh     = beta_uh
     191        self.beta_w = beta_w
     192        self.beta_w_dry = beta_w_dry
     193        self.beta_uh = beta_uh
    174194        self.beta_uh_dry = beta_uh_dry
    175         self.beta_vh     = beta_vh
     195        self.beta_vh = beta_vh
    176196        self.beta_vh_dry = beta_vh_dry
    177197        self.alpha_balance = alpha_balance
    … …  
    188208        self.set_store_vertices_uniquely(False)
    189209        self.minimum_storable_height = minimum_storable_height
    190         self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
     210        self.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    191211
    192212        # Limiters
    … …  
    195215        self.use_centroid_velocities = use_centroid_velocities
    196216
    197 
     217    ##
     218    # @brief
     219    # @param beta
    198220    def set_all_limiters(self, beta):
    199         """Shorthand to assign one constant value [0,1[ to all limiters.
     221        """Shorthand to assign one constant value [0,1] to all limiters.
    200222        0 Corresponds to first order, where as larger values make use of
    201         the second order scheme. 
    202         """
    203 
    204         self.beta_w      = beta
    205         self.beta_w_dry  = beta
     223        the second order scheme.
     224        """
     225
     226        self.beta_w = beta
     227        self.beta_w_dry = beta
    206228        self.quantities['stage'].beta = beta
    207        
    208         self.beta_uh     = beta
     229
     230        self.beta_uh = beta
    209231        self.beta_uh_dry = beta
    210232        self.quantities['xmomentum'].beta = beta
    211        
    212         self.beta_vh     = beta
     233
     234        self.beta_vh = beta
    213235        self.beta_vh_dry = beta
    214236        self.quantities['ymomentum'].beta = beta
    215        
    216        
    217 
     237
     238    ##
     239    # @brief
     240    # @param flag
     241    # @param reduction
    218242    def set_store_vertices_uniquely(self, flag, reduction=None):
    219243        """Decide whether vertex values should be stored uniquely as
    … …  
    230254            #self.reduction = min  #Looks better near steep slopes
    231255
    232 
     256    ##
     257    # @brief Set the minimum depth that will be written to an SWW file.
     258    # @param minimum_storable_height The minimum stored height (in m).
    233259    def set_minimum_storable_height(self, minimum_storable_height):
    234         """
    235         Set the minimum depth that will be recognised when writing
     260        """Set the minimum depth that will be recognised when writing
    236261        to an sww file. This is useful for removing thin water layers
    237262        that seems to be caused by friction creep.
    … …  
    239264        The minimum allowed sww depth is in meters.
    240265        """
     266
    241267        self.minimum_storable_height = minimum_storable_height
    242268
    243 
     269    ##
     270    # @brief
     271    # @param minimum_allowed_height
    244272    def set_minimum_allowed_height(self, minimum_allowed_height):
    245         """
    246         Set the minimum depth that will be recognised in the numerical
    247         scheme
     273        """Set the minimum depth recognised in the numerical scheme.
    248274
    249275        The minimum allowed depth is in meters.
    … …  
    258284        #and deal with them adaptively similarly to how we used to use 1 order
    259285        #steps to recover.
     286
    260287        self.minimum_allowed_height = minimum_allowed_height
    261         self.H0 = minimum_allowed_height   
    262        
    263 
     288        self.H0 = minimum_allowed_height
     289
     290    ##
     291    # @brief
     292    # @param maximum_allowed_speed
    264293    def set_maximum_allowed_speed(self, maximum_allowed_speed):
    265         """
    266         Set the maximum particle speed that is allowed in water
     294        """Set the maximum particle speed that is allowed in water
    267295        shallower than minimum_allowed_height. This is useful for
    268296        controlling speeds in very thin layers of water and at the same time
    269297        allow some movement avoiding pooling of water.
    270 
    271         """
     298        """
     299
    272300        self.maximum_allowed_speed = maximum_allowed_speed
    273301
    274 
    275     def set_points_file_block_line_size(self,points_file_block_line_size):
    276         """
    277         Set the minimum depth that will be recognised when writing
     302    ##
     303    # @brief
     304    # @param points_file_block_line_size
     305    def set_points_file_block_line_size(self, points_file_block_line_size):
     306        """Set the minimum depth that will be recognised when writing
    278307        to an sww file. This is useful for removing thin water layers
    279308        that seems to be caused by friction creep.
    … …  
    282311        """
    283312        self.points_file_block_line_size = points_file_block_line_size
    284        
    285        
     313
     314    ##
     315    # @brief Set the quantities that will be written to an SWW file.
     316    # @param q The quantities to be written.
     317    # @note Param 'q' may be None, single quantity or list of quantity strings.
     318    # @note If 'q' is None, no quantities will be stored in the SWW file.
    286319    def set_quantities_to_be_stored(self, q):
    287320        """Specify which quantities will be stored in the sww file.
    … …  
    295328        each yieldstep (This is in addition to the quantities elevation
    296329        and friction)
    297        
     330
    298331        If q is None, storage will be switched off altogether.
    299332        """
    300 
    301333
    302334        if q is None:
    … …  
    310342        # Check correcness
    311343        for quantity_name in q:
    312             msg = 'Quantity %s is not a valid conserved quantity'\
    313                   %quantity_name
    314            
     344            msg = ('Quantity %s is not a valid conserved quantity'
     345                   % quantity_name)
    315346            assert quantity_name in self.conserved_quantities, msg
    316347
    317348        self.quantities_to_be_stored = q
    318349
    319 
    320 
     350    ##
     351    # @brief
     352    # @param indices
    321353    def get_wet_elements(self, indices=None):
    322354        """Return indices for elements where h > minimum_allowed_height
    … …  
    328360            indices = get_wet_elements()
    329361
    330         Note, centroid values are used for this operation           
     362        Note, centroid values are used for this operation
    331363        """
    332364
    … …  
    335367        from anuga.config import minimum_allowed_height
    336368
    337        
    338369        elevation = self.get_quantity('elevation').\
     370                        get_values(location='centroids', indices=indices)
     371        stage = self.get_quantity('stage').\
    339372                    get_values(location='centroids', indices=indices)
    340         stage = self.get_quantity('stage').\
    341                 get_values(location='centroids', indices=indices)       
    342373        depth = stage - elevation
    343374
    … …  
    345376        wet_indices = num.compress(depth > minimum_allowed_height,
    346377                                   num.arange(len(depth)))
    347         return wet_indices
    348 
    349 
     378        return wet_indices
     379
     380    ##
     381    # @brief
     382    # @param indices
    350383    def get_maximum_inundation_elevation(self, indices=None):
    351384        """Return highest elevation where h > 0
    … …  
    357390            q = get_maximum_inundation_elevation()
    358391
    359         Note, centroid values are used for this operation           
     392        Note, centroid values are used for this operation
    360393        """
    361394
    362395        wet_elements = self.get_wet_elements(indices)
    363396        return self.get_quantity('elevation').\
    364                get_maximum_value(indices=wet_elements)
    365 
    366 
     397                   get_maximum_value(indices=wet_elements)
     398
     399    ##
     400    # @brief
     401    # @param indices
    367402    def get_maximum_inundation_location(self, indices=None):
    368403        """Return location of highest elevation where h > 0
    … …  
    374409            q = get_maximum_inundation_location()
    375410
    376         Note, centroid values are used for this operation           
     411        Note, centroid values are used for this operation
    377412        """
    378413
    379414        wet_elements = self.get_wet_elements(indices)
    380415        return self.get_quantity('elevation').\
    381                get_maximum_location(indices=wet_elements)   
    382                
    383                
    384                
    385                
    386     def get_flow_through_cross_section(self, polyline,
    387                                        verbose=False):               
    388         """Get the total flow through an arbitrary poly line.       
    389        
    390         This is a run-time equivalent of the function with same name
     416                   get_maximum_location(indices=wet_elements)
     417
     418    ##
     419    # @brief Get the total flow through an arbitrary poly line.
     420    # @param polyline Representation of desired cross section.
     421    # @param verbose True if this method is to be verbose.
     422    # @note 'polyline' may contain multiple sections allowing complex shapes.
     423    # @note Assume absolute UTM coordinates.
     424    def get_flow_through_cross_section(self, polyline, verbose=False):
     425        """Get the total flow through an arbitrary poly line.
     426
     427        This is a run-time equivalent of the function with same name
    391428        in data_manager.py
    392        
     429
    393430        Input:
    394             polyline: Representation of desired cross section - it may contain 
    395                       multiple sections allowing for complex shapes. Assume 
     431            polyline: Representation of desired cross section - it may contain
     432                      multiple sections allowing for complex shapes. Assume
    396433                      absolute UTM coordinates.
    397                       Format [[x0, y0], [x1, y1], ...]       
    398                  
    399         Output:       
     434                      Format [[x0, y0], [x1, y1], ...]
     435
     436        Output:
    400437            Q: Total flow [m^3/s] across given segments.
    401        
    402          
    403         """       
    404        
     438        """
     439
    405440        # Find all intersections and associated triangles.
    406         segments = self.get_intersecting_segments(polyline,
    407                                                   use_cache=True,
     441        segments = self.get_intersecting_segments(polyline, use_cache=True,
    408442                                                  verbose=verbose)
    409443
    410444        # Get midpoints
    411         midpoints = segment_midpoints(segments)       
    412        
     445        midpoints = segment_midpoints(segments)
     446
    413447        # Make midpoints Geospatial instances
    414         midpoints = ensure_geospatial(midpoints, self.geo_reference)       
    415        
    416         # Compute flow       
     448        midpoints = ensure_geospatial(midpoints, self.geo_reference)
     449
     450        # Compute flow
    417451        if verbose: print 'Computing flow through specified cross section'
    418        
     452
    419453        # Get interpolated values
    420454        xmomentum = self.get_quantity('xmomentum')
    421         ymomentum = self.get_quantity('ymomentum')       
    422        
    423         uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
    424         vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
    425        
     455        ymomentum = self.get_quantity('ymomentum')
     456
     457        uh = xmomentum.get_values(interpolation_points=midpoints,
     458                                  use_cache=True)
     459        vh = ymomentum.get_values(interpolation_points=midpoints,
     460                                  use_cache=True)
     461
    426462        # Compute and sum flows across each segment
    427         total_flow=0
     463        total_flow = 0
    428464        for i in range(len(uh)):
    429            
    430             # Inner product of momentum vector with segment normal [m^2/s]
     465            # Inner product of momentum vector with segment normal [m^2/s]
    431466            normal = segments[i].normal
    432             normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 
    433                
     467            normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
     468
    434469            # Flow across this segment [m^3/s]
    435470            segment_flow = normal_momentum*segments[i].length
    … …  
    437472            # Accumulate
    438473            total_flow += segment_flow
    439            
     474
    440475        return total_flow
    441476
    442        
    443                
     477    ##
     478    # @brief
     479    # @param polyline Representation of desired cross section.
     480    # @param kind Select energy type to compute ('specific' or 'total').
     481    # @param verbose True if this method is to be verbose.
     482    # @note 'polyline' may contain multiple sections allowing complex shapes.
     483    # @note Assume absolute UTM coordinates.
    444484    def get_energy_through_cross_section(self, polyline,
    445485                                         kind='total',
    446                                          verbose=False):               
     486                                         verbose=False):
    447487        """Obtain average energy head [m] across specified cross section.
    448488
    449489        Inputs:
    450             polyline: Representation of desired cross section - it may contain 
    451                       multiple sections allowing for complex shapes. Assume 
     490            polyline: Representation of desired cross section - it may contain
     491                      multiple sections allowing for complex shapes. Assume
    452492                      absolute UTM coordinates.
    453493                      Format [[x0, y0], [x1, y1], ...]
    454             kind:     Select which energy to compute. 
     494            kind:     Select which energy to compute.
    455495                      Options are 'specific' and 'total' (default)
    456496
    … …  
    458498            E: Average energy [m] across given segments for all stored times.
    459499
    460         The average velocity is computed for each triangle intersected by 
    461         the polyline and averaged weighted by segment lengths. 
    462 
    463         The typical usage of this function would be to get average energy of 
    464         flow in a channel, and the polyline would then be a cross section 
     500        The average velocity is computed for each triangle intersected by
     501        the polyline and averaged weighted by segment lengths.
     502
     503        The typical usage of this function would be to get average energy of
     504        flow in a channel, and the polyline would then be a cross section
    465505        perpendicular to the flow.
    466506
    467         #FIXME (Ole) - need name for this energy reflecting that its dimension 
     507        #FIXME (Ole) - need name for this energy reflecting that its dimension
    468508        is [m].
    469509        """
    470510
    471         from anuga.config import g, epsilon, velocity_protection as h0       
    472                                          
    473        
     511        from anuga.config import g, epsilon, velocity_protection as h0
     512
    474513        # Find all intersections and associated triangles.
    475         segments = self.get_intersecting_segments(polyline,
    476                                                   use_cache=True,
     514        segments = self.get_intersecting_segments(polyline, use_cache=True,
    477515                                                  verbose=verbose)
    478516
    479517        # Get midpoints
    480         midpoints = segment_midpoints(segments)       
    481        
     518        midpoints = segment_midpoints(segments)
     519
    482520        # Make midpoints Geospatial instances
    483         midpoints = ensure_geospatial(midpoints, self.geo_reference)       
    484        
    485         # Compute energy       
    486         if verbose: print 'Computing %s energy' %kind       
    487        
     521        midpoints = ensure_geospatial(midpoints, self.geo_reference)
     522
     523        # Compute energy
     524        if verbose: print 'Computing %s energy' %kind
     525
    488526        # Get interpolated values
    489         stage = self.get_quantity('stage')       
    490         elevation = self.get_quantity('elevation')               
     527        stage = self.get_quantity('stage')
     528        elevation = self.get_quantity('elevation')
    491529        xmomentum = self.get_quantity('xmomentum')
    492         ymomentum = self.get_quantity('ymomentum')       
     530        ymomentum = self.get_quantity('ymomentum')
    493531
    494532        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
    495         z = elevation.get_values(interpolation_points=midpoints, use_cache=True)       
    496         uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True)
    497         vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True)
    498         h = w-z # Depth
    499        
     533        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
     534        uh = xmomentum.get_values(interpolation_points=midpoints,
     535                                  use_cache=True)
     536        vh = ymomentum.get_values(interpolation_points=midpoints,
     537                                  use_cache=True)
     538        h = w-z                # Depth
     539
    500540        # Compute total length of polyline for use with weighted averages
    501541        total_line_length = 0.0
    502542        for segment in segments:
    503543            total_line_length += segment.length
    504        
     544
    505545        # Compute and sum flows across each segment
    506         average_energy=0.0
     546        average_energy = 0.0
    507547        for i in range(len(w)):
    508            
    509548            # Average velocity across this segment
    510549            if h[i] > epsilon:
    … …  
    514553            else:
    515554                u = v = 0.0
    516                
    517             speed_squared = u*u + v*v   
     555
     556            speed_squared = u*u + v*v
    518557            kinetic_energy = 0.5*speed_squared/g
    519            
     558
    520559            if kind == 'specific':
    521560                segment_energy = h[i] + kinetic_energy
    522561            elif kind == 'total':
    523                 segment_energy = w[i] + kinetic_energy               
     562                segment_energy = w[i] + kinetic_energy
    524563            else:
    525564                msg = 'Energy kind must be either "specific" or "total".'
    526565                msg += ' I got %s' %kind
    527                
    528566
    529567            # Add to weighted average
    530568            weigth = segments[i].length/total_line_length
    531569            average_energy += segment_energy*weigth
    532              
    533            
     570
    534571        return average_energy
    535        
    536 
    537                        
    538 
     572
     573    ##
     574    # @brief Run integrity checks on shallow water domain.
    539575    def check_integrity(self):
    540         Generic_Domain.check_integrity(self)
     576        Generic_Domain.check_integrity(self)    # super integrity check
    541577
    542578        #Check that we are solving the shallow water wave equation
    543 
    544579        msg = 'First conserved quantity must be "stage"'
    545580        assert self.conserved_quantities[0] == 'stage', msg
    … …  
    549584        assert self.conserved_quantities[2] == 'ymomentum', msg
    550585
     586    ##
     587    # @brief
    551588    def extrapolate_second_order_sw(self):
    552         #Call correct module function
    553         #(either from this module or C-extension)
     589        #Call correct module function (either from this module or C-extension)
    554590        extrapolate_second_order_sw(self)
    555591
     592    ##
     593    # @brief
    556594    def compute_fluxes(self):
    557         #Call correct module function
    558         #(either from this module or C-extension)
     595        #Call correct module function (either from this module or C-extension)
    559596        compute_fluxes(self)
    560597
     598    ##
     599    # @brief
    561600    def distribute_to_vertices_and_edges(self):
    562601        # Call correct module function
    563602        if self.use_edge_limiter:
    564             distribute_using_edge_limiter(self)           
     603            distribute_using_edge_limiter(self)
    565604        else:
    566605            distribute_using_vertex_limiter(self)
    567606
    568 
    569 
    570 
     607    ##
     608    # @brief Evolve the model by one step.
     609    # @param yieldstep
     610    # @param finaltime
     611    # @param duration
     612    # @param skip_initial_step
    571613    def evolve(self,
    572                yieldstep = None,
    573                finaltime = None,
    574                duration = None,
    575                skip_initial_step = False):
    576         """Specialisation of basic evolve method from parent class
    577         """
     614               yieldstep=None,
     615               finaltime=None,
     616               duration=None,
     617               skip_initial_step=False):
     618        """Specialisation of basic evolve method from parent class"""
    578619
    579620        # Call check integrity here rather than from user scripts
    580621        # self.check_integrity()
    581622
    582         msg = 'Parameter beta_w must be in the interval [0, 2['
     623        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
    583624        assert 0 <= self.beta_w <= 2.0, msg
    584625
    585 
    586626        # Initial update of vertex and edge values before any STORAGE
    587         # and or visualisation
     627        # and or visualisation.
    588628        # This is done again in the initialisation of the Generic_Domain
    589         # evolve loop but we do it here to ensure the values are ok for storage
     629        # evolve loop but we do it here to ensure the values are ok for storage.
    590630        self.distribute_to_vertices_and_edges()
    591631
    592632        if self.store is True and self.time == 0.0:
    593633            self.initialise_storage()
    594             # print 'Storing results in ' + self.writer.filename
    595634        else:
    596635            pass
    … …  
    601640
    602641        # Call basic machinery from parent class
    603         for t in Generic_Domain.evolve(self,
    604                                        yieldstep=yieldstep,
    605                                        finaltime=finaltime,
    606                                        duration=duration,
     642        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
     643                                       finaltime=finaltime, duration=duration,
    607644                                       skip_initial_step=skip_initial_step):
    608 
    609645            # Store model data, e.g. for subsequent visualisation
    610646            if self.store is True:
    … …  
    617653            yield(t)
    618654
    619 
     655    ##
     656    # @brief
    620657    def initialise_storage(self):
    621658        """Create and initialise self.writer object for storing data.
    … …  
    631668        self.writer.store_connectivity()
    632669
    633 
     670    ##
     671    # @brief
     672    # @param name
    634673    def store_timestep(self, name):
    635674        """Store named quantity and time.
    … …  
    638677           self.write has been initialised
    639678        """
     679
    640680        self.writer.store_timestep(name)
    641681
    642        
     682    ##
     683    # @brief Get time stepping statistics string for printing.
     684    # @param track_speeds
     685    # @param triangle_id
    643686    def timestepping_statistics(self,
    644687                                track_speeds=False,
    645                                 triangle_id=None):       
     688                                triangle_id=None):
    646689        """Return string with time stepping statistics for printing or logging
    647690
    … …  
    651694        """
    652695
    653         from anuga.config import epsilon, g               
    654 
     696        from anuga.config import epsilon, g
    655697
    656698        # Call basic machinery from parent class
    657         msg = Generic_Domain.timestepping_statistics(self,
    658                                                      track_speeds,
     699        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
    659700                                                     triangle_id)
    660701
    661702        if track_speeds is True:
    662 
    663703            # qwidth determines the text field used for quantities
    664704            qwidth = self.qwidth
    665        
     705
    666706            # Selected triangle
    667707            k = self.k
    … …  
    669709            # Report some derived quantities at vertices, edges and centroid
    670710            # specific to the shallow water wave equation
    671 
    672711            z = self.quantities['elevation']
    673             w = self.quantities['stage']           
     712            w = self.quantities['stage']
    674713
    675714            Vw = w.get_values(location='vertices', indices=[k])[0]
    … …  
    679718            Vz = z.get_values(location='vertices', indices=[k])[0]
    680719            Ez = z.get_values(location='edges', indices=[k])[0]
    681             Cz = z.get_values(location='centroids', indices=[k])                             
    682                
     720            Cz = z.get_values(location='centroids', indices=[k])
    683721
    684722            name = 'depth'
    … …  
    686724            Eh = Ew-Ez
    687725            Ch = Cw-Cz
    688            
     726
    689727            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    690728                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
    691            
     729
    692730            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    693731                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
    694            
     732
    695733            s += '    %s: centroid_value = %.4f\n'\
    696                  %(name.ljust(qwidth), Ch[0])                               
    697            
     734                 %(name.ljust(qwidth), Ch[0])
     735
    698736            msg += s
    699737
    … …  
    704742            Euh = uh.get_values(location='edges', indices=[k])[0]
    705743            Cuh = uh.get_values(location='centroids', indices=[k])
    706            
     744
    707745            Vvh = vh.get_values(location='vertices', indices=[k])[0]
    708746            Evh = vh.get_values(location='edges', indices=[k])[0]
    … …  
    712750            Vu = Vuh/(Vh + epsilon)
    713751            Eu = Euh/(Eh + epsilon)
    714             Cu = Cuh/(Ch + epsilon)           
     752            Cu = Cuh/(Ch + epsilon)
    715753            name = 'U'
    716754            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    717755                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
    718            
     756
    719757            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    720758                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
    721            
     759
    722760            s += '    %s: centroid_value = %.4f\n'\
    723                  %(name.ljust(qwidth), Cu[0])                               
    724            
     761                 %(name.ljust(qwidth), Cu[0])
     762
    725763            msg += s
    726764
    727765            Vv = Vvh/(Vh + epsilon)
    728766            Ev = Evh/(Eh + epsilon)
    729             Cv = Cvh/(Ch + epsilon)           
     767            Cv = Cvh/(Ch + epsilon)
    730768            name = 'V'
    731769            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    732770                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
    733            
     771
    734772            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    735773                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
    736            
     774
    737775            s += '    %s: centroid_value = %.4f\n'\
    738                  %(name.ljust(qwidth), Cv[0])                               
    739            
     776                 %(name.ljust(qwidth), Cv[0])
     777
    740778            msg += s
    741 
    742779
    743780            # Froude number in each direction
    … …  
    746783            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
    747784            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
    748            
     785
    749786            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    750787                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
    751            
     788
    752789            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    753790                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
    754            
     791
    755792            s += '    %s: centroid_value = %.4f\n'\
    756                  %(name.ljust(qwidth), Cfx[0])                               
    757            
     793                 %(name.ljust(qwidth), Cfx[0])
     794
    758795            msg += s
    759 
    760796
    761797            name = 'Froude (y)'
    … …  
    763799            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
    764800            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
    765            
     801
    766802            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
    767803                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
    768            
     804
    769805            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
    770806                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
    771            
     807
    772808            s += '    %s: centroid_value = %.4f\n'\
    773                  %(name.ljust(qwidth), Cfy[0])                               
    774            
    775             msg += s           
    776 
    777                
     809                 %(name.ljust(qwidth), Cfy[0])
     810
     811            msg += s
    778812
    779813        return msg
    780        
    781        
    782 
    783 #=============== End of class Shallow Water Domain ===============================
    784 
     814
     815################################################################################
     816# End of class Shallow Water Domain
     817################################################################################
    785818
    786819#-----------------
    … …  
    788821#-----------------
    789822
     823## @brief Compute fluxes and timestep suitable for all volumes in domain.
     824# @param domain The domain to calculate fluxes for.
    790825def compute_fluxes(domain):
    791     """Compute all fluxes and the timestep suitable for all volumes
    792     in domain.
     826    """Compute fluxes and timestep suitable for all volumes in domain.
    793827
    794828    Compute total flux for each conserved quantity using "flux_function"
    … …  
    806840      domain.explicit_update is reset to computed flux values
    807841      domain.timestep is set to the largest step satisfying all volumes.
    808    
    809842
    810843    This wrapper calls the underlying C version of compute fluxes
    … …  
    812845
    813846    import sys
    814 
    815     N = len(domain) # number_of_triangles
     847    from shallow_water_ext import compute_fluxes_ext_central \
     848                                  as compute_fluxes_ext
     849
     850    N = len(domain)    # number_of_triangles
    816851
    817852    # Shortcuts
    … …  
    822857
    823858    timestep = float(sys.maxint)
    824     from shallow_water_ext import\
    825          compute_fluxes_ext_central as compute_fluxes_ext
    826 
    827859
    828860    flux_timestep = compute_fluxes_ext(timestep,
    … …  
    853885    domain.flux_timestep = flux_timestep
    854886
    855 
    856 
    857 #---------------------------------------
     887################################################################################
    858888# Module functions for gradient limiting
    859 #---------------------------------------
    860 
    861 
    862 # MH090605 The following method belongs to the shallow_water domain class
    863 # see comments in the corresponding method in shallow_water_ext.c
     889################################################################################
     890
     891##
     892# @brief Wrapper for C version of extrapolate_second_order_sw.
     893# @param domain The domain to operate on.
     894# @note MH090605 The following method belongs to the shallow_water domain class
     895#       see comments in the corresponding method in shallow_water_ext.c
    864896def extrapolate_second_order_sw(domain):
    865     """Wrapper calling C version of extrapolate_second_order_sw
    866     """
     897    """Wrapper calling C version of extrapolate_second_order_sw"""
     898
    867899    import sys
     900    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
    868901
    869902    N = len(domain) # number_of_triangles
    … …  
    875908    Elevation = domain.quantities['elevation']
    876909
    877     from shallow_water_ext import extrapolate_second_order_sw as extrapol2
    878910    extrapol2(domain,
    879911              domain.surrogate_neighbours,
    … …  
    891923              int(domain.optimise_dry_cells))
    892924
    893 
     925##
     926# @brief Distribution from centroids to vertices specific to the SWW eqn.
     927# @param domain The domain to operate on.
    894928def distribute_using_vertex_limiter(domain):
    895     """Distribution from centroids to vertices specific to the
    896     shallow water wave
    897     equation.
     929    """Distribution from centroids to vertices specific to the SWW equation.
    898930
    899931    It will ensure that h (w-z) is always non-negative even in the
    … …  
    912944    Postcondition
    913945      Conserved quantities defined at vertices
    914 
    915946    """
    916 
    917    
    918947
    919948    # Remove very thin layers of water
    … …  
    946975                raise 'Unknown order'
    947976
    948 
    949977    # Take bed elevation into account when water heights are small
    950978    balance_deep_and_shallow(domain)
    … …  
    955983        Q.interpolate_from_vertices_to_edges()
    956984
    957 
    958 
     985##
     986# @brief Distribution from centroids to edges specific to the SWW eqn.
     987# @param domain The domain to operate on.
    959988def distribute_using_edge_limiter(domain):
    960     """Distribution from centroids to edges specific to the
    961     shallow water wave
    962     equation.
     989    """Distribution from centroids to edges specific to the SWW eqn.
    963990
    964991    It will ensure that h (w-z) is always non-negative even in the
    … …  
    9761003    Postcondition
    9771004      Conserved quantities defined at vertices
    978 
    9791005    """
    9801006
    9811007    # Remove very thin layers of water
    9821008    protect_against_infinitesimal_and_negative_heights(domain)
    983 
    9841009
    9851010    for name in domain.conserved_quantities:
    … …  
    9991024        Q.interpolate_from_vertices_to_edges()
    10001025
    1001 
     1026##
     1027# @brief  Protect against infinitesimal heights and associated high velocities.
     1028# @param domain The domain to operate on.
    10021029def protect_against_infinitesimal_and_negative_heights(domain):
    1003     """Protect against infinitesimal heights and associated high velocities
    1004     """
     1030    """Protect against infinitesimal heights and associated high velocities"""
     1031
     1032    from shallow_water_ext import protect
    10051033
    10061034    # Shortcuts
    … …  
    10101038    ymomc = domain.quantities['ymomentum'].centroid_values
    10111039
    1012     from shallow_water_ext import protect
    1013 
    10141040    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
    10151041            domain.epsilon, wc, zc, xmomc, ymomc)
    10161042
    1017 
    1018 
     1043##
     1044# @brief Wrapper for C function balance_deep_and_shallow_c().
     1045# @param domain The domain to operate on.
    10191046def balance_deep_and_shallow(domain):
    10201047    """Compute linear combination between stage as computed by
    … …  
    10311058    """
    10321059
    1033     from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c
    1034 
     1060    from shallow_water_ext import balance_deep_and_shallow \
     1061                                  as balance_deep_and_shallow_c
    10351062
    10361063    # Shortcuts
    10371064    wc = domain.quantities['stage'].centroid_values
    10381065    zc = domain.quantities['elevation'].centroid_values
    1039 
    10401066    wv = domain.quantities['stage'].vertex_values
    10411067    zv = domain.quantities['elevation'].vertex_values
    … …  
    10501076
    10511077    balance_deep_and_shallow_c(domain,
    1052                                wc, zc, wv, zv, wc, 
     1078                               wc, zc, wv, zv, wc,
    10531079                               xmomc, ymomc, xmomv, ymomv)
    10541080
    10551081
    1056 
    1057 
    1058 #------------------------------------------------------------------
     1082################################################################################
    10591083# Boundary conditions - specific to the shallow water wave equation
    1060 #------------------------------------------------------------------
     1084################################################################################
     1085
     1086##
     1087# @brief Class for a reflective boundary.
     1088# @note Inherits from Boundary.
    10611089class Reflective_boundary(Boundary):
    10621090    """Reflective boundary returns same conserved quantities as
    … …  
    10681096    """
    10691097
    1070     def __init__(self, domain = None):
     1098    ##
     1099    # @brief Instantiate a Reflective_boundary.
     1100    # @param domain
     1101    def __init__(self, domain=None):
    10711102        Boundary.__init__(self)
    10721103
    10731104        if domain is None:
    10741105            msg = 'Domain must be specified for reflective boundary'
    1075             raise msg
     1106            raise Exception, msg
    10761107
    10771108        # Handy shorthands
    1078         self.stage   = domain.quantities['stage'].edge_values
    1079         self.xmom    = domain.quantities['xmomentum'].edge_values
    1080         self.ymom    = domain.quantities['ymomentum'].edge_values
     1109        self.stage = domain.quantities['stage'].edge_values
     1110        self.xmom = domain.quantities['xmomentum'].edge_values
     1111        self.ymom = domain.quantities['ymomentum'].edge_values
    10811112        self.normals = domain.normals
    10821113
    10831114        self.conserved_quantities = num.zeros(3, num.float)
    10841115
     1116    ##
     1117    # @brief Return a representation of this instance.
    10851118    def __repr__(self):
    10861119        return 'Reflective_boundary'
    10871120
    1088 
     1121    ##
     1122    # @brief Calculate reflections (reverse outward momentum).
     1123    # @param vol_id
     1124    # @param edge_id
    10891125    def evaluate(self, vol_id, edge_id):
    10901126        """Reflective boundaries reverses the outward momentum
    … …  
    10991135        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
    11001136
    1101 
    11021137        r = rotate(q, normal, direction = 1)
    11031138        r[1] = -r[1]
    … …  
    11071142
    11081143
    1109 
    1110 
     1144##
     1145# @brief Class for a transmissive boundary.
     1146# @note Inherits from Boundary.
    11111147class Transmissive_momentum_set_stage_boundary(Boundary):
    11121148    """Returns same momentum conserved quantities as
    … …  
    11171153    Example:
    11181154
    1119     def waveform(t): 
     1155    def waveform(t):
    11201156        return sea_level + normalized_amplitude/cosh(t-25)**2
    11211157
    11221158    Bts = Transmissive_momentum_set_stage_boundary(domain, waveform)
    1123    
    11241159
    11251160    Underlying domain must be specified when boundary is instantiated
    11261161    """
    11271162
    1128     def __init__(self, domain = None, function=None):
     1163    ##
     1164    # @brief Instantiate a Reflective_boundary.
     1165    # @param domain
     1166    # @param function
     1167    def __init__(self, domain=None, function=None):
    11291168        Boundary.__init__(self)
    11301169
    11311170        if domain is None:
    11321171            msg = 'Domain must be specified for this type boundary'
    1133             raise msg
     1172            raise Exception, msg
    11341173
    11351174        if function is None:
    11361175            msg = 'Function must be specified for this type boundary'
    1137             raise msg
    1138 
    1139         self.domain   = domain
     1176            raise Exception, msg
     1177
     1178        self.domain = domain
    11401179        self.function = function
    11411180
     1181    ##
     1182    # @brief Return a representation of this instance.
    11421183    def __repr__(self):
    11431184        return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain
    11441185
     1186    ##
     1187    # @brief Calculate transmissive results.
     1188    # @param vol_id
     1189    # @param edge_id
    11451190    def evaluate(self, vol_id, edge_id):
    11461191        """Transmissive momentum set stage boundaries return the edge momentum
    … …  
    11501195        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
    11511196
    1152 
    11531197        t = self.domain.get_time()
    11541198
    11551199        if hasattr(self.function, 'time'):
    1156             # Roll boundary over if time exceeds           
     1200            # Roll boundary over if time exceeds
    11571201            while t > self.function.time[-1]:
    1158                 msg = 'WARNING: domain time %.2f has exceeded' %t
     1202                msg = 'WARNING: domain time %.2f has exceeded' % t
    11591203                msg += 'time provided in '
    11601204                msg += 'transmissive_momentum_set_stage boundary object.\n'
    … …  
    11631207                t -= self.function.time[-1]
    11641208
    1165 
    11661209        value = self.function(t)
    11671210        try:
    … …  
    11691212        except:
    11701213            x = float(value[0])
    1171            
     1214
    11721215        q[0] = x
    11731216        return q
    1174 
    11751217
    11761218        # FIXME: Consider this (taken from File_boundary) to allow
    … …  
    11831225
    11841226
    1185 # Backward compatibility       
    1186 # FIXME(Ole): Deprecate
     1227##
     1228# @brief Deprecated boundary class.
    11871229class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary):
    11881230    pass
    11891231
    1190      
    1191 
     1232
     1233##
     1234# @brief A transmissive boundary, momentum set to zero.
     1235# @note Inherits from Bouondary.
    11921236class Transmissive_stage_zero_momentum_boundary(Boundary):
    1193     """Return same stage as those present in its neighbour volume. Set momentum to zero.
     1237    """Return same stage as those present in its neighbour volume.
     1238    Set momentum to zero.
    11941239
    11951240    Underlying domain must be specified when boundary is instantiated
    11961241    """
    11971242
     1243    ##
     1244    # @brief Instantiate a Transmissive (zero momentum) boundary.
     1245    # @param domain
    11981246    def __init__(self, domain=None):
    11991247        Boundary.__init__(self)
    12001248
    12011249        if domain is None:
    1202             msg = 'Domain must be specified for '
    1203             msg += 'Transmissive_stage_zero_momentum boundary'
     1250            msg = ('Domain must be specified for '
     1251                   'Transmissive_stage_zero_momentum boundary')
    12041252            raise Exception, msg
    12051253
    12061254        self.domain = domain
    12071255
     1256    ##
     1257    # @brief Return a representation of this instance.
    12081258    def __repr__(self):
    1209         return 'Transmissive_stage_zero_momentum_boundary(%s)' %self.domain
    1210 
     1259        return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain
     1260
     1261    ##
     1262    # @brief Calculate transmissive (zero momentum) results.
     1263    # @param vol_id
     1264    # @param edge_id
    12111265    def evaluate(self, vol_id, edge_id):
    12121266        """Transmissive boundaries return the edge values
    … …  
    12151269
    12161270        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
    1217        
     1271
    12181272        q[1] = q[2] = 0.0
    12191273        return q
    12201274
    12211275
    1222        
     1276##
     1277# @brief Class for a Dirichlet discharge boundary.
     1278# @note Inherits from Boundary.
    12231279class Dirichlet_discharge_boundary(Boundary):
    12241280    """
    … …  
    12291285    """
    12301286
     1287    ##
     1288    # @brief Instantiate a Dirichlet discharge boundary.
     1289    # @param domain
     1290    # @param stage0
     1291    # @param wh0
    12311292    def __init__(self, domain=None, stage0=None, wh0=None):
    12321293        Boundary.__init__(self)
    12331294
    12341295        if domain is None:
    1235             msg = 'Domain must be specified for this type boundary'
    1236             raise msg
     1296            msg = 'Domain must be specified for this type of boundary'
     1297            raise Exception, msg
    12371298
    12381299        if stage0 is None:
    1239             raise 'set stage'
     1300            raise Exception, 'Stage must be specified for this type of boundary'
    12401301
    12411302        if wh0 is None:
    12421303            wh0 = 0.0
    12431304
    1244         self.domain   = domain
    1245         self.stage0  = stage0
     1305        self.domain = domain
     1306        self.stage0 = stage0
    12461307        self.wh0 = wh0
    12471308
     1309    ##
     1310    # @brief Return a representation of this instance.
    12481311    def __repr__(self):
    1249         return 'Dirichlet_Discharge_boundary(%s)' %self.domain
    1250 
     1312        return 'Dirichlet_Discharge_boundary(%s)' % self.domain
     1313
     1314    ##
     1315    # @brief Calculate Dirichlet discharge boundary results.
     1316    # @param vol_id
     1317    # @param edge_id
    12511318    def evaluate(self, vol_id, edge_id):
    1252         """Set discharge in the (inward) normal direction
    1253         """
     1319        """Set discharge in the (inward) normal direction"""
    12541320
    12551321        normal = self.domain.get_normal(vol_id,edge_id)
    12561322        q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]]
    12571323        return q
    1258 
    12591324
    12601325        # FIXME: Consider this (taken from File_boundary) to allow
    … …  
    12671332
    12681333
    1269        
    1270 # Backward compatibility       
     1334# Backward compatibility
    12711335# FIXME(Ole): Deprecate
     1336##
     1337# @brief Deprecated
    12721338class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary):
    12731339    pass
    1274                                                    
    1275    
    1276 
    1277        
    1278        
     1340
     1341
     1342##
     1343# @brief A class for a 'field' boundary.
     1344# @note Inherits from Boundary.
    12791345class Field_boundary(Boundary):
    1280     """Set boundary from given field represented in an sww file containing values
    1281     for stage, xmomentum and ymomentum.
    1282     Optionally, the user can specify mean_stage to offset the stage provided in the
    1283     sww file.
    1284 
    1285     This function is a thin wrapper around the generic File_boundary. The
     1346    """Set boundary from given field represented in an sww file containing
     1347    values for stage, xmomentum and ymomentum.
     1348
     1349    Optionally, the user can specify mean_stage to offset the stage provided
     1350    in the sww file.
     1351
     1352    This function is a thin wrapper around the generic File_boundary. The
    12861353    difference between the file_boundary and field_boundary is only that the
    12871354    field_boundary will allow you to change the level of the stage height when
    1288     you read in the boundary condition. This is very useful when running 
    1289     different tide heights in the same area as you need only to convert one 
    1290     boundary condition to a SWW file, ideally for tide height of 0 m 
     1355    you read in the boundary condition. This is very useful when running
     1356    different tide heights in the same area as you need only to convert one
     1357    boundary condition to a SWW file, ideally for tide height of 0 m
    12911358    (saving disk space). Then you can use field_boundary to read this SWW file
    12921359    and change the stage height (tide) on the fly depending on the scenario.
    1293    
    12941360    """
    12951361
    1296 
    1297     def __init__(self, filename, domain,
     1362    ##
     1363    # @brief Construct an instance of a 'field' boundary.
     1364    # @param filename Name of SWW file containing stage, x and ymomentum.
     1365    # @param domain Shallow water domain for which the boundary applies.
     1366    # @param mean_stage Mean water level added to stage derived from SWW file.
     1367    # @param time_thinning Time step thinning factor.
     1368    # @param time_limit
     1369    # @param boundary_polygon
     1370    # @param default_boundary None or an instance of Boundary.
     1371    # @param use_cache True if caching is to be used.
     1372    # @param verbose True if this method is to be verbose.
     1373    def __init__(self,
     1374                 filename,
     1375                 domain,
    12981376                 mean_stage=0.0,
    12991377                 time_thinning=1,
    13001378                 time_limit=None,
    1301                  boundary_polygon=None,   
    1302                  default_boundary=None,                 
     1379                 boundary_polygon=None,
     1380                 default_boundary=None,
    13031381                 use_cache=False,
    13041382                 verbose=False):
    … …  
    13101388                    from the sww file
    13111389        time_thinning: Will set how many time steps from the sww file read in
    1312                        will be interpolated to the boundary. For example if 
     1390                       will be interpolated to the boundary. For example if
    13131391                       the sww file has 1 second time steps and is 24 hours
    1314                        in length it has 86400 time steps. If you set 
    1315                        time_thinning to 1 it will read all these steps. 
     1392                       in length it has 86400 time steps. If you set
     1393                       time_thinning to 1 it will read all these steps.
    13161394                       If you set it to 100 it will read every 100th step eg
    13171395                       only 864 step. This parameter is very useful to increase
    1318                        the speed of a model run that you are setting up 
     1396                       the speed of a model run that you are setting up
    13191397                       and testing.
    1320                        
    1321         default_boundary: Must be either None or an instance of a 
     1398
     1399        default_boundary: Must be either None or an instance of a
    13221400                          class descending from class Boundary.
    1323                           This will be used in case model time exceeds 
     1401                          This will be used in case model time exceeds
    13241402                          that available in the underlying data.
    1325                                                
     1403
    13261404        use_cache:
    13271405        verbose:
    1328        
    13291406        """
    13301407
    13311408        # Create generic file_boundary object
    1332         self.file_boundary = File_boundary(filename, domain,
     1409        self.file_boundary = File_boundary(filename,
     1410                                           domain,
    13331411                                           time_thinning=time_thinning,
    13341412                                           time_limit=time_limit,
    … …  
    13381416                                           verbose=verbose)
    13391417
    1340        
    13411418        # Record information from File_boundary
    13421419        self.F = self.file_boundary.F
    13431420        self.domain = self.file_boundary.domain
    1344        
     1421
    13451422        # Record mean stage
    13461423        self.mean_stage = mean_stage
    13471424
    1348 
     1425    ##
     1426    # @note Generate a string representation of this instance.
    13491427    def __repr__(self):
    13501428        return 'Field boundary'
    13511429
    1352 
     1430    ##
     1431    # @brief Calculate 'field' boundary results.
     1432    # @param vol_id
     1433    # @param edge_id
    13531434    def evaluate(self, vol_id=None, edge_id=None):
    13541435        """Return linearly interpolated values based on domain.time
    … …  
    13561437        vol_id and edge_id are ignored
    13571438        """
    1358        
     1439
    13591440        # Evaluate file boundary
    13601441        q = self.file_boundary.evaluate(vol_id, edge_id)
    … …  
    13661447        return q
    13671448
    1368    
    1369 
    1370 #-----------------------
     1449
     1450################################################################################
    13711451# Standard forcing terms
    1372 #-----------------------
    1373 
     1452################################################################################
     1453
     1454##
     1455# @brief Apply gravitational pull in the presence of bed slope.
     1456# @param domain The domain to apply gravity to.
     1457# @note Wrapper for C function gravity_c().
    13741458def gravity(domain):
    13751459    """Apply gravitational pull in the presence of bed slope
    … …  
    13771461    """
    13781462
     1463    from shallow_water_ext import gravity as gravity_c
     1464
    13791465    xmom = domain.quantities['xmomentum'].explicit_update
    13801466    ymom = domain.quantities['ymomentum'].explicit_update
    … …  
    13881474    x = domain.get_vertex_coordinates()
    13891475    g = domain.g
    1390    
    1391 
    1392     from shallow_water_ext import gravity as gravity_c
    1393     gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6)
    1394 
    1395 
    1396 
     1476
     1477    gravity_c(g, h, z, x, xmom, ymom)    #, 1.0e-6)
     1478
     1479##
     1480# @brief Apply friction to a surface (implicit).
     1481# @param domain The domain to apply Manning friction to.
     1482# @note Wrapper for C function manning_friction_c().
    13971483def manning_friction_implicit(domain):
    1398     """Apply (Manning) friction to water momentum   
     1484    """Apply (Manning) friction to water momentum
    13991485    Wrapper for c version
    14001486    """
    14011487
    1402 
    1403     #print 'Implicit friction'
     1488    from shallow_water_ext import manning_friction as manning_friction_c
    14041489
    14051490    xmom = domain.quantities['xmomentum']
    … …  
    14201505    g = domain.g
    14211506
    1422     from shallow_water_ext import manning_friction as manning_friction_c
    14231507    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
    14241508
    1425 
     1509##
     1510# @brief Apply friction to a surface (explicit).
     1511# @param domain The domain to apply Manning friction to.
     1512# @note Wrapper for C function manning_friction_c().
    14261513def manning_friction_explicit(domain):
    1427     """Apply (Manning) friction to water momentum   
     1514    """Apply (Manning) friction to water momentum
    14281515    Wrapper for c version
    14291516    """
    14301517
    1431     # print 'Explicit friction'
     1518    from shallow_water_ext import manning_friction as manning_friction_c
    14321519
    14331520    xmom = domain.quantities['xmomentum']
    … …  
    14481535    g = domain.g
    14491536
    1450     from shallow_water_ext import manning_friction as manning_friction_c
    14511537    manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
    14521538
    14531539
    14541540# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
    1455 # Is it still needed (30 Oct 2007)?
     1541##
     1542# @brief Apply linear friction to a surface.
     1543# @param domain The domain to apply Manning friction to.
     1544# @note Is this still used (30 Oct 2007)?
    14561545def linear_friction(domain):
    14571546    """Apply linear friction to water momentum
    … …  
    14861575                ymom_update[k] += S*vh[k]
    14871576
    1488 
    1489 
    1490 #---------------------------------
     1577################################################################################
    14911578# Experimental auxiliary functions
    1492 #---------------------------------
     1579################################################################################
     1580
     1581##
     1582# @brief Check forcefield parameter.
     1583# @param f Object to check.
     1584# @note 'f' may be a callable object or a scalar value.
    14931585def check_forcefield(f):
    1494     """Check that f is either
     1586    """Check that force object is as expected.
     1587   
     1588    Check that f is either:
    14951589    1: a callable object f(t,x,y), where x and y are vectors
    14961590       and that it returns an array or a list of same length
    … …  
    15161610            msg += 'not be converted into a numeric array of floats.\n'
    15171611            msg += 'Specified function should return either list or array.'
    1518             raise msg
     1612            raise Exception, msg
    15191613
    15201614        # Is this really what we want?
    1521         msg = 'Return vector from function %s ' %f
    1522         msg += 'must have same length as input vectors'
    1523         msg += ' (type(q)=%s' % type(q)
     1615        msg = 'Function %s must return vector' % str(f)
     1616        assert hasattr(q, 'len'), msg
     1617
     1618        msg = ('Return vector from function %s must have same '
     1619               'length as input vectors' % f)
    15241620        assert len(q) == N, msg
    1525 
    15261621    else:
    15271622        try:
    15281623            f = float(f)
    15291624        except:
    1530             msg = 'Force field %s must be either a scalar' %f
    1531             msg += ' or a vector function'
    1532             raise Exception(msg)
     1625            msg = ('Force field %s must be either a vector function or a '
     1626                   'scalar value (coercible to float).' % str(f))
     1627            raise Exception, msg
     1628
    15331629    return f
    15341630
    15351631
     1632##
     1633# Class to apply a wind stress to a domain.
    15361634class Wind_stress:
    15371635    """Apply wind stress to water momentum in terms of
    … …  
    15391637    """
    15401638
     1639    ##
     1640    # @brief Create an instance of Wind_stress.
     1641    # @param *args
     1642    # @param **kwargs
    15411643    def __init__(self, *args, **kwargs):
    15421644        """Initialise windfield from wind speed s [m/s]
    … …  
    15911693        else:
    15921694           # Assume info is in 2 keyword arguments
    1593 
    15941695           if len(kwargs) == 2:
    15951696               s = kwargs['s']
    15961697               phi = kwargs['phi']
    15971698           else:
    1598                raise 'Assumes two keyword arguments: s=..., phi=....'
     1699               raise Exception, 'Assumes two keyword arguments: s=..., phi=....'
    15991700
    16001701        self.speed = check_forcefield(s)
    … …  
    16031704        self.const = eta_w*rho_a/rho_w
    16041705
    1605 
     1706    ##
     1707    # @brief 'execute' this class instance.
     1708    # @param domain
    16061709    def __call__(self, domain):
    1607         """Evaluate windfield based on values found in domain
    1608         """
     1710        """Evaluate windfield based on values found in domain"""
    16091711
    16101712        from math import pi, cos, sin, sqrt
    … …  
    16131715        ymom_update = domain.quantities['ymomentum'].explicit_update
    16141716
    1615         N = len(domain) # number_of_triangles
     1717        N = len(domain)    # number_of_triangles
    16161718        t = domain.time
    16171719
    … …  
    16211723        else:
    16221724            # Assume s is a scalar
    1623 
    16241725            try:
    16251726                s_vec = self.speed * num.ones(N, num.float)
    … …  
    16281729                raise msg
    16291730
    1630 
    16311731        if callable(self.phi):
    16321732            xc = domain.get_centroid_coordinates()
    … …  
    16451745
    16461746
     1747##
     1748# @brief Assign wind field values
     1749# @param xmom_update
     1750# @param ymom_update
     1751# @param s_vec
     1752# @param phi_vec
     1753# @param const
    16471754def assign_windfield_values(xmom_update, ymom_update,
    16481755                            s_vec, phi_vec, const):
    … …  
    16501757    A c version also exists (for speed)
    16511758    """
     1759
    16521760    from math import pi, cos, sin, sqrt
    16531761
    … …  
    16701778
    16711779
    1672 
    1673 
    1674 
     1780##
     1781# @brief A class for a general explicit forcing term.
    16751782class General_forcing:
    16761783    """General explicit forcing term for update of quantity
    1677    
     1784
    16781785    This is used by Inflow and Rainfall for instance
    1679    
     1786
    16801787
    16811788    General_forcing(quantity_name, rate, center, radius, polygon)
    16821789
    16831790    domain:     ANUGA computational domain
    1684     quantity_name: Name of quantity to update. 
     1791    quantity_name: Name of quantity to update.
    16851792                   It must be a known conserved quantity.
    1686                    
    1687     rate [?/s]: Total rate of change over the specified area. 
     1793
     1794    rate [?/s]: Total rate of change over the specified area.
    16881795                This parameter can be either a constant or a
    1689                 function of time. Positive values indicate increases, 
     1796                function of time. Positive values indicate increases,
    16901797                negative values indicate decreases.
    16911798                Rate can be None at initialisation but must be specified
    … …  
    17011808    Either center, radius or polygon can be specified but not both.
    17021809    If neither are specified the entire domain gets updated.
    1703     All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.   
    1704    
     1810    All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain.
     1811
    17051812    Inflow or Rainfall for examples of use
    17061813    """
    … …  
    17091816    # FIXME (AnyOne) : Add various methods to allow spatial variations
    17101817
     1818    ##
     1819    # @brief Create an instance of this forcing term.
     1820    # @param domain
     1821    # @param quantity_name
     1822    # @param rate
     1823    # @param center
     1824    # @param radius
     1825    # @param polygon
     1826    # @param default_rate
     1827    # @param verbose
    17111828    def __init__(self,
    17121829                 domain,
    17131830                 quantity_name,
    17141831                 rate=0.0,
    1715                  center=None, radius=None,
     1832                 center=None,
     1833                 radius=None,
    17161834                 polygon=None,
    17171835                 default_rate=None,
    17181836                 verbose=False):
    1719                      
     1837
     1838        from math import pi, cos, sin
     1839
    17201840        if center is None:
    1721             msg = 'I got radius but no center.'       
     1841            msg = 'I got radius but no center.'
    17221842            assert radius is None, msg
    1723            
     1843
    17241844        if radius is None:
    1725             msg += 'I got center but no radius.'       
     1845            msg += 'I got center but no radius.'
    17261846            assert center is None, msg
    1727            
    1728            
    1729                      
    1730         from math import pi, cos, sin
    17311847
    17321848        self.domain = domain
    … …  
    17351851        self.center = ensure_numeric(center)
    17361852        self.radius = radius
    1737         self.polygon = polygon       
     1853        self.polygon = polygon
    17381854        self.verbose = verbose
    1739         self.value = 0.0 # Can be used to remember value at
    1740                          # previous timestep in order to obtain rate
     1855        self.value = 0.0    # Can be used to remember value at
     1856                            # previous timestep in order to obtain rate
    17411857
    17421858        # Get boundary (in absolute coordinates)
    17431859        bounding_polygon = domain.get_boundary_polygon()
    17441860
    1745 
    17461861        # Update area if applicable
    1747         self.exchange_area = None       
     1862        self.exchange_area = None
    17481863        if center is not None and radius is not None:
    17491864            assert len(center) == 2
    … …  
    17541869
    17551870            # Check that circle center lies within the mesh.
    1756             msg = 'Center %s specified for forcing term did not' %(str(center))
     1871            msg = 'Center %s specified for forcing term did not' % str(center)
    17571872            msg += 'fall within the domain boundary.'
    17581873            assert is_inside_polygon(center, bounding_polygon), msg
    … …  
    17621877            periphery_points = []
    17631878            for i in range(N):
    1764 
    17651879                theta = 2*pi*i/100
    1766                
     1880
    17671881                x = center[0] + radius*cos(theta)
    17681882                y = center[1] + radius*sin(theta)
    … …  
    17701884                periphery_points.append([x,y])
    17711885
    1772 
    17731886            for point in periphery_points:
    1774                 msg = 'Point %s on periphery for forcing term' %(str(point))
     1887                msg = 'Point %s on periphery for forcing term' % str(point)
    17751888                msg += ' did not fall within the domain boundary.'
    17761889                assert is_inside_polygon(point, bounding_polygon), msg
    17771890
    17781891        if polygon is not None:
    1779 
    17801892            # Check that polygon lies within the mesh.
    17811893            for point in self.polygon:
    1782                 msg = 'Point %s in polygon for forcing term' %(point)
     1894                msg = 'Point %s in polygon for forcing term' % str(point)
    17831895                msg += ' did not fall within the domain boundary.'
    17841896                assert is_inside_polygon(point, bounding_polygon), msg
    1785                
    1786             # Compute area and check that it is greater than 0   
     1897
     1898            # Compute area and check that it is greater than 0
    17871899            self.exchange_area = polygon_area(self.polygon)
    1788            
    1789             msg = 'Polygon %s in forcing term' %(self.polygon)
    1790             msg += ' has area = %f' %self.exchange_area
    1791             assert self.exchange_area > 0.0           
    1792            
    1793                
    1794                
     1900
     1901            msg = 'Polygon %s in forcing term' % str(self.polygon)
     1902            msg += ' has area = %f' % self.exchange_area
     1903            assert self.exchange_area > 0.0
    17951904
    17961905        # Pointer to update vector
    … …  
    17981907
    17991908        # Determine indices in flow area
    1800         N = len(domain)   
     1909        N = len(domain)
    18011910        points = domain.get_centroid_coordinates(absolute=True)
    18021911
    … …  
    18051914        if self.center is not None and self.radius is not None:
    18061915            # Inlet is circular
    1807            
    1808             inlet_region = 'center=%s, radius=%s' %(self.center, self.radius)
    1809            
     1916            inlet_region = 'center=%s, radius=%s' % (self.center, self.radius)
     1917
    18101918            self.exchange_indices = []
    18111919            for k in range(N):
    1812                 x, y = points[k,:] # Centroid
    1813                
     1920                x, y = points[k,:]    # Centroid
     1921
    18141922                c = self.center
    18151923                if ((x-c[0])**2+(y-c[1])**2) < self.radius**2:
    18161924                    self.exchange_indices.append(k)
    1817                    
    1818         if self.polygon is not None:                   
     1925
     1926        if self.polygon is not None:
    18191927            # Inlet is polygon
    1820            
    1821             inlet_region = 'polygon=%s, area=%f m^2' %(self.polygon,
    1822                                                        self.exchange_area)
    1823                        
     1928            inlet_region = 'polygon=%s, area=%f m^2' % (self.polygon,
     1929                                                        self.exchange_area)
     1930
    18241931            self.exchange_indices = inside_polygon(points, self.polygon)
    1825            
    1826            
    1827            
     1932
    18281933        if self.exchange_indices is not None:
    1829             #print inlet_region
    1830        
    18311934            if len(self.exchange_indices) == 0:
    18321935                msg = 'No triangles have been identified in '
    1833                 msg += 'specified region: %s' %inlet_region
     1936                msg += 'specified region: %s' % inlet_region
    18341937                raise Exception, msg
    1835            
     1938
    18361939        # Check and store default_rate
    1837         msg = 'Keyword argument default_rate must be either None '
    1838         msg += 'or a function of time.\n I got %s' %(str(default_rate))
    1839         assert default_rate is None or \
    1840                type(default_rate) in [IntType, FloatType] or \
    1841                callable(default_rate), msg
    1842        
     1940        msg = ('Keyword argument default_rate must be either None '
     1941               'or a function of time.\nI got %s.' % str(default_rate))
     1942        assert (default_rate is None or
     1943                type(default_rate) in [IntType, FloatType] or
     1944                callable(default_rate)), msg
     1945
    18431946        if default_rate is not None:
    1844 
    18451947            # If it is a constant, make it a function
    18461948            if not callable(default_rate):
    … …  
    18481950                default_rate = lambda t: tmp
    18491951
    1850            
    18511952            # Check that default_rate is a function of one argument
    18521953            try:
    … …  
    18561957
    18571958        self.default_rate = default_rate
    1858         self.default_rate_invoked = False    # Flag       
    1859        
    1860 
     1959        self.default_rate_invoked = False    # Flag
     1960
     1961    ##
     1962    # @brief Execute this instance.
     1963    # @param domain
    18611964    def __call__(self, domain):
    1862         """Apply inflow function at time specified in domain and update stage
    1863         """
     1965        """Apply inflow function at time specified in domain, update stage"""
    18641966
    18651967        # Call virtual method allowing local modifications
    1866 
    18671968        t = domain.get_time()
    18681969        try:
    … …  
    18721973        except Modeltime_too_late, e:
    18731974            if self.default_rate is None:
    1874                 raise Exception, e # Reraise exception
     1975                raise Exception, e    # Reraise exception
    18751976            else:
    18761977                # Pass control to default rate function
    18771978                rate = self.default_rate(t)
    1878                
     1979
    18791980                if self.default_rate_invoked is False:
    18801981                    # Issue warning the first time
    1881                     msg = '%s' %str(e)
    1882                     msg += 'Instead I will use the default rate: %s\n'\
    1883                         %str(self.default_rate)
    1884                     msg += 'Note: Further warnings will be supressed'
     1982                    msg = ('%s\n'
     1983                           'Instead I will use the default rate: %s\n'
     1984                           'Note: Further warnings will be supressed'
     1985                           % (str(e), str(self.default_rate)))
    18851986                    warn(msg)
    1886                    
     1987
    18871988                    # FIXME (Ole): Replace this crude flag with
    18881989                    # Python's ability to print warnings only once.
    18891990                    # See http://docs.python.org/lib/warning-filter.html
    18901991                    self.default_rate_invoked = True
    1891                    
    1892 
    1893            
    1894                
    18951992
    18961993        if rate is None:
    1897             msg = 'Attribute rate must be specified in General_forcing'
    1898             msg += ' or its descendants before attempting to call it'
     1994            msg = ('Attribute rate must be specified in General_forcing '
     1995                   'or its descendants before attempting to call it')
    18991996            raise Exception, msg
    1900        
    19011997
    19021998        # Now rate is a number
    19031999        if self.verbose is True:
    1904             print 'Rate of %s at time = %.2f = %f' %(self.quantity_name,
    1905                                                      domain.get_time(),
    1906                                                      rate)
    1907 
     2000            print 'Rate of %s at time = %.2f = %f' % (self.quantity_name,
     2001                                                      domain.get_time(),
     2002                                                      rate)
    19082003
    19092004        if self.exchange_indices is None:
    … …  
    19142009                self.update[k] += rate
    19152010
    1916 
     2011    ##
     2012    # @brief Update the internal rate.
     2013    # @param t A callable or scalar used to set the rate.
     2014    # @return The new rate.
    19172015    def update_rate(self, t):
    19182016        """Virtual method allowing local modifications by writing an
    19192017        overriding version in descendant
    1920        
    1921         """
     2018        """
     2019
    19222020        if callable(self.rate):
    19232021            rate = self.rate(t)
    … …  
    19272025        return rate
    19282026
    1929 
     2027    ##
     2028    # @brief Get values for the specified quantity.
     2029    # @param quantity_name Name of the quantity of interest.
     2030    # @return The value(s) of the quantity.
     2031    # @note If 'quantity_name' is None, use self.quantity_name.
    19302032    def get_quantity_values(self, quantity_name=None):
    1931         """Return values for specified quantity restricted to opening 
    1932        
     2033        """Return values for specified quantity restricted to opening
     2034
    19332035        Optionally a quantity name can be specified if values from another
    19342036        quantity is sought
    19352037        """
    1936        
     2038
    19372039        if quantity_name is None:
    19382040            quantity_name = self.quantity_name
    1939            
     2041
    19402042        q = self.domain.quantities[quantity_name]
    1941         return q.get_values(location='centroids', 
     2043        return q.get_values(location='centroids',
    19422044                            indices=self.exchange_indices)
    1943    
    1944 
     2045
     2046    ##
     2047    # @brief Set value for the specified quantity.
     2048    # @param val The value object used to set value.
     2049    # @param quantity_name Name of the quantity of interest.
     2050    # @note If 'quantity_name' is None, use self.quantity_name.
    19452051    def set_quantity_values(self, val, quantity_name=None):
    1946         """Set values for specified quantity restricted to opening 
    1947        
     2052        """Set values for specified quantity restricted to opening
     2053
    19482054        Optionally a quantity name can be specified if values from another
    1949         quantity is sought       
     2055        quantity is sought
    19502056        """
    19512057
    19522058        if quantity_name is None:
    19532059            quantity_name = self.quantity_name
    1954                    
    1955         q = self.domain.quantities[self.quantity_name]               
    1956         q.set_values(val,
    1957                      location='centroids',
    1958                      indices=self.exchange_indices)   
    1959 
    1960 
    1961 
     2060
     2061        q = self.domain.quantities[self.quantity_name]
     2062        q.set_values(val,
     2063                     location='centroids',
     2064                     indices=self.exchange_indices)
     2065
     2066
     2067##
     2068# @brief A class for rainfall forcing function.
     2069# @note Inherits from General_forcing.
    19622070class Rainfall(General_forcing):
    19632071    """Class Rainfall - general 'rain over entire domain' forcing term.
    1964    
     2072
    19652073    Used for implementing Rainfall over the entire domain.
    1966        
     2074
    19672075        Current Limited to only One Gauge..
    1968        
    1969         Need to add Spatial Varying Capability 
     2076
     2077        Need to add Spatial Varying Capability
    19702078        (This module came from copying and amending the Inflow Code)
    1971    
     2079
    19722080    Rainfall(rain)
    19732081
    1974     domain   
    1975     rain [mm/s]:  Total rain rate over the specified domain. 
     2082    domain
     2083    rain [mm/s]:  Total rain rate over the specified domain.
    19762084                  NOTE: Raingauge Data needs to reflect the time step.
    19772085                  IE: if Gauge is mm read at a time step, then the input
    … …  
    19802088
    19812089                  This parameter can be either a constant or a
    1982                   function of time. Positive values indicate inflow, 
     2090                  function of time. Positive values indicate inflow,
    19832091                  negative values indicate outflow.
    19842092                  (and be used for Infiltration - Write Seperate Module)
    … …  
    19882096
    19892097    polygon: Specifies a polygon to restrict the rainfall.
    1990    
     2098
    19912099    Examples
    19922100    How to put them in a run File...
    … …  
    19982106    # input of Rainfall in mm/s
    19992107
    2000     catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 
     2108    catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms'))
    20012109                        # Note need path to File in String.
    20022110                        # Else assumed in same directory
    … …  
    20052113    """
    20062114
    2007    
     2115    ##
     2116    # @brief Create an instance of the class.
     2117    # @param domain Domain of interest.
     2118    # @param rate Total rain rate over the specified domain (mm/s).
     2119    # @param center
     2120    # @param radius
     2121    # @param polygon Polygon  to restrict rainfall.
     2122    # @param default_rate
     2123    # @param verbose True if this instance is to be verbose.
    20082124    def __init__(self,
    20092125                 domain,
    20102126                 rate=0.0,
    2011                  center=None, radius=None,
     2127                 center=None,
     2128                 radius=None,
    20122129                 polygon=None,
    2013                  default_rate=None,                 
     2130                 default_rate=None,
    20142131                 verbose=False):
    20152132
    … …  
    20202137            rain = rate/1000.0
    20212138
    2022         if default_rate is not None:   
     2139        if default_rate is not None:
    20232140            if callable(default_rate):
    20242141                default_rain = lambda t: default_rate(t)/1000.0
    … …  
    20282145            default_rain = None
    20292146
    2030 
    2031            
    2032            
     2147        # pass to parent class
    20332148        General_forcing.__init__(self,
    20342149                                 domain,
    20352150                                 'stage',
    20362151                                 rate=rain,
    2037                                  center=center, radius=radius,
     2152                                 center=center,
     2153                                 radius=radius,
    20382154                                 polygon=polygon,
    20392155                                 default_rate=default_rain,
    20402156                                 verbose=verbose)
    20412157
    2042        
    2043 
    2044 
    2045 
    2046 
     2158
     2159##
     2160# @brief A class for inflow (rain and drain) forcing function.
     2161# @note Inherits from General_forcing.
    20472162class Inflow(General_forcing):
    20482163    """Class Inflow - general 'rain and drain' forcing term.
    2049    
     2164
    20502165    Useful for implementing flows in and out of the domain.
    2051    
     2166
    20522167    Inflow(flow, center, radius, polygon)
    20532168
    20542169    domain
    2055     rate [m^3/s]: Total flow rate over the specified area. 
     2170    rate [m^3/s]: Total flow rate over the specified area.
    20562171                  This parameter can be either a constant or a
    2057                   function of time. Positive values indicate inflow, 
     2172                  function of time. Positive values indicate inflow,
    20582173                  negative values indicate outflow.
    20592174                  The specified flow will be divided by the area of
    2060                   the inflow region and then applied to update stage.     
     2175                  the inflow region and then applied to update stage.
    20612176    center [m]: Coordinates at center of flow point
    20622177    radius [m]: Size of circular area
    … …  
    20642179
    20652180    Either center, radius or polygon must be specified
    2066    
     2181
    20672182    Examples
    20682183
    20692184    # Constant drain at 0.003 m^3/s.
    20702185    # The outflow area is 0.07**2*pi=0.0154 m^2
    2071     # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s 
    2072     #         
     2186    # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s
     2187    #
    20732188    Inflow((0.7, 0.4), 0.07, -0.003)
    20742189
    … …  
    20762191    # Tap turning up to a maximum inflow of 0.0142 m^3/s.
    20772192    # The inflow area is 0.03**2*pi = 0.00283 m^2
    2078     # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s     
     2193    # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s
    20792194    # over the specified area
    20802195    Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142))
    … …  
    20912206
    20922207    domain.forcing_terms.append(hydrograph)
    2093    
    20942208    """
    20952209
    2096 
     2210    ##
     2211    # @brief Create an instance of the class.
     2212    # @param domain Domain of interest.
     2213    # @param rate Total rain rate over the specified domain (mm/s).
     2214    # @param center
     2215    # @param radius
     2216    # @param polygon Polygon  to restrict rainfall.
     2217    # @param default_rate
     2218    # @param verbose True if this instance is to be verbose.
    20972219    def __init__(self,
    20982220                 domain,
    20992221                 rate=0.0,
    2100                  center=None, radius=None,
     2222                 center=None,
     2223                 radius=None,
    21012224                 polygon=None,
    21022225                 default_rate=None,
    2103                  verbose=False):                 
    2104 
    2105 
     2226                 verbose=False):
    21062227        # Create object first to make area is available
    21072228        General_forcing.__init__(self,
    … …  
    21092230                                 'stage',
    21102231                                 rate=rate,
    2111                                  center=center, radius=radius,
     2232                                 center=center,
     2233                                 radius=radius,
    21122234                                 polygon=polygon,
    21132235                                 default_rate=default_rate,
    21142236                                 verbose=verbose)
    21152237
     2238    ##
     2239    # @brief Update the instance rate.
     2240    # @param t New rate object.
    21162241    def update_rate(self, t):
    21172242        """Virtual method allowing local modifications by writing an
    21182243        overriding version in descendant
    21192244
    2120         This one converts m^3/s to m/s which can be added directly 
     2245        This one converts m^3/s to m/s which can be added directly
    21212246        to 'stage' in ANUGA
    21222247        """
    … …  
    21302255
    21312256
    2132 
    2133 
    2134 #------------------
     2257################################################################################
    21352258# Initialise module
    2136 #------------------
    2137 
     2259################################################################################
    21382260
    21392261from anuga.utilities import compile
    21402262if compile.can_use_C_extension('shallow_water_ext.c'):
    2141     # Underlying C implementations can be accessed
    2142 
     2263    # Underlying C implementations can be accessed
    21432264    from shallow_water_ext import rotate, assign_windfield_values
    21442265else:
    2145     msg = 'C implementations could not be accessed by %s.\n ' %__file__
     2266    msg = 'C implementations could not be accessed by %s.\n ' % __file__
    21462267    msg += 'Make sure compile_all.py has been run as described in '
    21472268    msg += 'the ANUGA installation guide.'
    21482269    raise Exception, msg
    2149 
    21502270
    21512271# Optimisation with psyco
    … …  
    21602280            #Psyco isn't supported on 64 bit systems, but it doesn't matter
    21612281        else:
    2162             msg = 'WARNING: psyco (speedup) could not import'+\
    2163                   ', you may want to consider installing it'
     2282            msg = ('WARNING: psyco (speedup) could not be imported, '
     2283                   'you may want to consider installing it')
    21642284            print msg
    21652285    else:
    … …  
    21672287        psyco.bind(Domain.compute_fluxes)
    21682288
     2289
    21692290if __name__ == "__main__":
    21702291    pass
    2171 
    2172 

  • branches/numpy/anuga/shallow_water/shallow_water_ext.c

    r6304 r6410  
    879879  }
    880880
     881  // check that numpy array objects arrays are C contiguous memory
     882  CHECK_C_CONTIG(h);
     883  CHECK_C_CONTIG(v);
     884  CHECK_C_CONTIG(x);
     885  CHECK_C_CONTIG(xmom);
     886  CHECK_C_CONTIG(ymom);
     887
    881888  N = h -> dimensions[0];
    882889  for (k=0; k<N; k++) {
    … …  
    937944  }
    938945
     946  // check that numpy array objects arrays are C contiguous memory
     947  CHECK_C_CONTIG(w);
     948  CHECK_C_CONTIG(z);
     949  CHECK_C_CONTIG(uh);
     950  CHECK_C_CONTIG(vh);
     951  CHECK_C_CONTIG(eta);
     952  CHECK_C_CONTIG(xmom);
     953  CHECK_C_CONTIG(ymom);
    939954
    940955  N = w -> dimensions[0];
    … …  
    967982                        &xmom, &ymom))
    968983    return NULL;
     984
     985  // check that numpy array objects arrays are C contiguous memory
     986  CHECK_C_CONTIG(w);
     987  CHECK_C_CONTIG(z);
     988  CHECK_C_CONTIG(uh);
     989  CHECK_C_CONTIG(vh);
     990  CHECK_C_CONTIG(eta);
     991  CHECK_C_CONTIG(xmom);
     992  CHECK_C_CONTIG(ymom);
    969993
    970994  N = w -> dimensions[0];
    … …  
    14991523  }
    15001524
     1525  // check that numpy array objects arrays are C contiguous memory
     1526  CHECK_C_CONTIG(surrogate_neighbours);
     1527  CHECK_C_CONTIG(number_of_boundaries);
     1528  CHECK_C_CONTIG(centroid_coordinates);
     1529  CHECK_C_CONTIG(stage_centroid_values);
     1530  CHECK_C_CONTIG(xmom_centroid_values);
     1531  CHECK_C_CONTIG(ymom_centroid_values);
     1532  CHECK_C_CONTIG(elevation_centroid_values);
     1533  CHECK_C_CONTIG(vertex_coordinates);
     1534  CHECK_C_CONTIG(stage_vertex_values);
     1535  CHECK_C_CONTIG(xmom_vertex_values);
     1536  CHECK_C_CONTIG(ymom_vertex_values);
     1537  CHECK_C_CONTIG(elevation_vertex_values);
     1538 
    15011539  // Get the safety factor beta_w, set in the config.py file.
    15021540  // This is used in the limiting process
    … …  
    20522090  }
    20532091
    2054  
     2092  // check that numpy array objects arrays are C contiguous memory
     2093  CHECK_C_CONTIG(neighbours);
     2094  CHECK_C_CONTIG(neighbour_edges);
     2095  CHECK_C_CONTIG(normals);
     2096  CHECK_C_CONTIG(edgelengths);
     2097  CHECK_C_CONTIG(radii);
     2098  CHECK_C_CONTIG(areas);
     2099  CHECK_C_CONTIG(tri_full_flag);
     2100  CHECK_C_CONTIG(stage_edge_values);
     2101  CHECK_C_CONTIG(xmom_edge_values);
     2102  CHECK_C_CONTIG(ymom_edge_values);
     2103  CHECK_C_CONTIG(bed_edge_values);
     2104  CHECK_C_CONTIG(stage_boundary_values);
     2105  CHECK_C_CONTIG(xmom_boundary_values);
     2106  CHECK_C_CONTIG(ymom_boundary_values);
     2107  CHECK_C_CONTIG(stage_explicit_update);
     2108  CHECK_C_CONTIG(xmom_explicit_update);
     2109  CHECK_C_CONTIG(ymom_explicit_update);
     2110  CHECK_C_CONTIG(already_computed_flux);
     2111  CHECK_C_CONTIG(max_speed_array);
     2112
    20552113  int number_of_elements = stage_edge_values -> dimensions[0];
    20562114

  • branches/numpy/anuga/shallow_water/test_all.py

    r6304 r6410  
    7676    return unittest.TestSuite(map(load, modules))
    7777
     78################################################################################
     79
    7880if __name__ == '__main__':   
     81
    7982    suite = regressionTest()
    8083    runner = unittest.TextTestRunner() #verbosity=2)

  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6304 r6410  
    75177517        #print 2.0*num.transpose(ha) - stage
    75187518       
    7519         ha_permutation = num.take(ha, permutation)
    7520         ua_permutation = num.take(ua, permutation)         
    7521         va_permutation = num.take(va, permutation)                 
    7522         gauge_depth_permutation = num.take(gauge_depth, permutation)                         
     7519        ha_permutation = num.take(ha, permutation, axis=0)
     7520        ua_permutation = num.take(ua, permutation, axis=0)
     7521        va_permutation = num.take(va, permutation, axis=0)
     7522        gauge_depth_permutation = num.take(gauge_depth, permutation, axis=0)
    75237523
    75247524       
    … …  
    75357535            #2.0*ha necessary because using two files with weights=1 are used
    75367536           
    7537         depth_permutation = num.take(depth, permutation)                     
     7537        depth_permutation = num.take(depth, permutation, axis=0)                     
    75387538       
    75397539
    … …  
    78587858        # quantities written to the mux2 files subject to the permutation vector.
    78597859       
    7860         ha_ref = num.take(ha0, permutation)
    7861         ua_ref = num.take(ua0, permutation)       
    7862         va_ref = num.take(va0, permutation)               
    7863 
    7864         gauge_depth_ref = num.take(gauge_depth, permutation)                     
     7860        ha_ref = num.take(ha0, permutation, axis=0)
     7861        ua_ref = num.take(ua0, permutation, axis=0)       
     7862        va_ref = num.take(va0, permutation, axis=0)               
     7863
     7864        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                     
    78657865       
    78667866        assert num.allclose(num.transpose(ha_ref)+tide, stage0)  # Meters
    … …  
    79667966        # quantities written to the mux2 files subject to the permutation vector.
    79677967       
    7968         ha_ref = num.take(ha1, permutation)
    7969         ua_ref = num.take(ua1, permutation)       
    7970         va_ref = num.take(va1, permutation)               
    7971 
    7972         gauge_depth_ref = num.take(gauge_depth, permutation)                         
     7968        ha_ref = num.take(ha1, permutation, axis=0)
     7969        ua_ref = num.take(ua1, permutation, axis=0)       
     7970        va_ref = num.take(va1, permutation, axis=0)               
     7971
     7972        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
    79737973
    79747974
    … …  
    80748074        # quantities written to the mux2 files subject to the permutation vector.
    80758075       
    8076         ha_ref = weights[0]*num.take(ha0, permutation) + weights[1]*num.take(ha1, permutation)
    8077         ua_ref = weights[0]*num.take(ua0, permutation) + weights[1]*num.take(ua1, permutation)       
    8078         va_ref = weights[0]*num.take(va0, permutation) + weights[1]*num.take(va1, permutation)               
    8079 
    8080         gauge_depth_ref = num.take(gauge_depth, permutation)                         
     8076        ha_ref = (weights[0]*num.take(ha0, permutation, axis=0)
     8077                  + weights[1]*num.take(ha1, permutation, axis=0))
     8078        ua_ref = (weights[0]*num.take(ua0, permutation, axis=0)
     8079                  + weights[1]*num.take(ua1, permutation, axis=0))
     8080        va_ref = (weights[0]*num.take(va0, permutation, axis=0)
     8081                  + weights[1]*num.take(va1, permutation, axis=0))
     8082
     8083        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
    80818084
    80828085
    … …  
    1058810591        h = stage-z
    1058910592        for i in range(len(stage)):
    10590             if h[i] == 0.0:
    10591                 assert xmomentum[i] == 0.0
    10592                 assert ymomentum[i] == 0.0               
     10593            if num.alltrue(h[i] == 0.0):
     10594                assert num.alltrue(xmomentum[i] == 0.0)
     10595                assert num.alltrue(ymomentum[i] == 0.0)
    1059310596            else:
    10594                 assert h[i] >= domain.minimum_storable_height
     10597                msg = ('h[i]=\n%s\ndomain.minimum_storable_height=%s'
     10598                       % (str(h[i]), str(domain.minimum_storable_height)))
     10599                assert num.alltrue(h[i] >= domain.minimum_storable_height), msg
    1059510600       
    1059610601        fid.close()
    … …  
    1132911334           
    1133011335#-------------------------------------------------------------
     11336
    1133111337if __name__ == "__main__":
    11332 
    1133311338    suite = unittest.makeSuite(Test_Data_Manager,'test')
    11334     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts')
    11335     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    11336     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
    11337     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources')
    11338     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources')
     11339    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_maximum_inundation')
    1133911340
    1134011341    # FIXME (Ole): This is the test that fails under Windows
    1134111342    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2')
    1134211343    #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV')
    11343 
    1134411344   
    1134511345    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':

  • branches/numpy/anuga/shallow_water/test_eq.py

    r6304 r6410  
    6060
    6161#-------------------------------------------------------------
     62
    6263if __name__ == "__main__":
    6364    suite = unittest.makeSuite(Test_eq,'test_Okada_func')

  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6304 r6410  
    199199
    200200    return 17.7
     201
     202
     203def scalar_func_list(t,x,y):
     204    """Function that returns a scalar.
     205    Used to test error message when numeric array is expected
     206    """
     207
     208    return [17.7]
    201209
    202210
    … …  
    24622470
    24632471        try:
    2464             domain.forcing_terms.append(Wind_stress(s = scalar_func,
     2472            domain.forcing_terms.append(Wind_stress(s = scalar_func_list,
    24652473                                                    phi = angle))
    24662474        except AssertionError:
    … …  
    24742482            domain.forcing_terms.append(Wind_stress(s = speed,
    24752483                                                    phi = scalar_func))
    2476         except AssertionError:
     2484        except Exception:
    24772485            pass
    24782486        else:
    … …  
    55825590        #print take(cv2, (0,3,8))
    55835591
    5584         assert num.allclose( num.take(cv1, (0,8,16)), num.take(cv2, (0,3,8))) #Diag
    5585         assert num.allclose( num.take(cv1, (0,6,12)), num.take(cv2, (0,1,4))) #Bottom
    5586         assert num.allclose( num.take(cv1, (12,14,16)), num.take(cv2, (4,6,8))) #RHS
     5592        assert num.allclose(num.take(cv1, (0,8,16), axis=0),
     5593                            num.take(cv2, (0,3,8), axis=0))    #Diag
     5594        assert num.allclose(num.take(cv1, (0,6,12), axis=0),
     5595                            num.take(cv2, (0,1,4), axis=0))    #Bottom
     5596        assert num.allclose(num.take(cv1, (12,14,16), axis=0),
     5597                            num.take(cv2, (4,6,8), axis=0))    #RHS
    55875598
    55885599        #Cleanup
    … …  
    56905701
    56915702        #print points[0], points[5], points[10], points[15]
    5692         msg = ('value was %s,\nshould be [[0,0], [1.0/3, 1.0/3], '
    5693                '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15])))
    5694         assert num.allclose(num.take(points, [0,5,10,15]),
     5703        msg = ('value was\n%s\nshould be\n'
     5704               '[[0,0], [1.0/3, 1.0/3],\n'
     5705               '[2.0/3, 2.0/3], [1,1]]'
     5706               % str(num.take(points, [0,5,10,15], axis=0)))
     5707        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
    56955708                            [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg
    56965709
    … …  
    58735886
    58745887        #print points[0], points[5], points[10], points[15]
    5875         msg = ('values was %s,\nshould be [[0,0], [1.0/3, 1.0/3], '
    5876                '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15])))
    5877         assert num.allclose(num.take(points, [0,5,10,15]),
     5888        msg = ('values was\n%s\nshould be\n'
     5889               '[[0,0], [1.0/3, 1.0/3],\n'
     5890               '[2.0/3, 2.0/3], [1,1]]'
     5891               % str(num.take(points, [0,5,10,15], axis=0)))
     5892        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
    58785893                            [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg
    58795894
    … …  
    60486063
    60496064        #print points[0], points[5], points[10], points[15]
    6050         msg = ('value was %s,\nshould be [[0,0], [1.0/3, 1.0/3], '
    6051                '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15])))
    6052         assert num.allclose(num.take(points, [0,5,10,15]),
     6065        msg = ('value was\n%s\nshould be\n'
     6066               '[[0,0], [1.0/3, 1.0/3],\n'
     6067               '[2.0/3, 2.0/3], [1,1]]'
     6068               % str(num.take(points, [0,5,10,15], axis=0)))
     6069        assert num.allclose(num.take(points, [0,5,10,15], axis=0),
    60536070                            [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg
    60546071
    … …  
    66266643       
    66276644if __name__ == "__main__":
    6628 
    6629     suite = unittest.makeSuite(Test_Shallow_Water, 'test')
    6630     #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
    6631     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g')   
    6632     #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')   
    6633     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
    6634     #suite = unittest.makeSuite(Test_Shallow_Water,'test_inflow_outflow_conservation')
    6635     #suite = unittest.makeSuite(Test_Shallow_Water,'test_outflow_conservation_problem_temp')   
    6636    
    6637 
    6638    
     6645    #suite = unittest.makeSuite(Test_Shallow_Water, 'test')
     6646    suite = unittest.makeSuite(Test_Shallow_Water, 'test_get_maximum_inundation_from_sww')
    66396647    runner = unittest.TextTestRunner(verbosity=1)   
    66406648    runner.run(suite)

  • branches/numpy/anuga/shallow_water/test_smf.py

    r6304 r6410  
    132132
    133133#-------------------------------------------------------------
     134
    134135if __name__ == "__main__":
    135     #suite = unittest.makeSuite(Test_smf,'test_Double_gaussian')
    136136    suite = unittest.makeSuite(Test_smf,'test')
    137137    runner = unittest.TextTestRunner()

  • branches/numpy/anuga/shallow_water/test_system.py

    r6304 r6410  
    187187       
    188188#-------------------------------------------------------------
     189
    189190if __name__ == "__main__":
    190191    suite = unittest.makeSuite(Test_system,'test')
    191     #suite = unittest.makeSuite(Test_system,'test_boundary_timeII')
    192192    runner = unittest.TextTestRunner()
    193193    runner.run(suite)

  • branches/numpy/anuga/shallow_water/test_tsunami_okada.py

    r6304 r6410  
    291291
    292292#-------------------------------------------------------------
     293
    293294if __name__ == "__main__":
    294295    suite = unittest.makeSuite(Test_eq,'test')

  • branches/numpy/anuga/shallow_water/tsunami_okada.py

    r6304 r6410  
    215215            yrec=x
    216216            for i in range(0,len(zrec[0])):
    217                 if zrec[0][i]==yrec and zrec[1][i]==xrec:
     217                if (num.alltrue(zrec[0][i]==yrec) and
     218                    num.alltrue(zrec[1][i]==xrec)):
    218219                    Z=zrec[2][i]
    219220                    Z=0.001*Z

  • branches/numpy/anuga/utilities/polygon.py

    r6360 r6410  
    11#!/usr/bin/env python
    2 """Polygon manipulations
    3 """
     2
     3"""Polygon manipulations"""
    44
    55import numpy as num
    … …  
    1111
    1212
     13##
     14# @brief Determine whether a point is on a line segment.
     15# @param point (x, y) of point in question (tuple, list or array).
     16# @param line ((x1,y1), (x2,y2)) for line (tuple, list or array).
     17# @param rtol Relative error for 'close'.
     18# @param atol Absolute error for 'close'.
     19# @return True or False.
    1320def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8):
    1421    """Determine whether a point is on a line segment
    1522
    16     Input: 
     23    Input:
    1724        point is given by [x, y]
    1825        line is given by [x0, y0], [x1, y1]] or
    … …  
    2128    Output:
    2229
    23     Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding.
     30    Note: Line can be degenerate and function still works to discern coinciding
     31          points from non-coinciding.
    2432    """
    2533
    … …  
    3139                         line[1,0], line[1,1],
    3240                         rtol, atol)
    33    
     41
    3442    return bool(res)
    3543
    … …  
    4250# result functions for possible states
    4351def lines_dont_coincide(p0,p1,p2,p3):               return (3, None)
    44 def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2, num.array([p0,p1]))
    45 def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2, num.array([p2,p3]))
    46 def lines_overlap_same_direction(p0,p1,p2,p3):      return (2, num.array([p0,p3]))
    47 def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2, num.array([p2,p1]))
    48 def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2, num.array([p0,p2]))
    49 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, num.array([p3,p1]))
     52def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2,
     53                                                            num.array([p0,p1]))
     54def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2,
     55                                                            num.array([p2,p3]))
     56def lines_overlap_same_direction(p0,p1,p2,p3):      return (2,
     57                                                            num.array([p0,p3]))
     58def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2,
     59                                                            num.array([p2,p1]))
     60def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2,
     61                                                            num.array([p0,p2]))
     62def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2,
     63                                                            num.array([p3,p1]))
    5064
    5165# this function called when an impossible state is found
    52 def lines_error(p1, p2, p3, p4): raise RuntimeError, ("INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s" %
    53                                                       (str(p1), str(p2), str(p3), str(p4)))
     66def lines_error(p1, p2, p3, p4):
     67    raise RuntimeError, ('INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s'
     68                         % (str(p1), str(p2), str(p3), str(p4)))
    5469
    5570#                     0s1    0e1    1s0    1e0   # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0
    … …  
    7287                   }
    7388
     89##
     90# @brief Finds intersection point of two line segments.
     91# @param line0 First line ((x1,y1), (x2,y2)).
     92# @param line1 Second line ((x1,y1), (x2,y2)).
     93# @param rtol Relative error for 'close'.
     94# @param atol Absolute error for 'close'.
     95# @return (status, value) where:
     96#             status = 0 - no intersection, value set to None
     97#                      1 - intersection found, value=(x,y)
     98#                      2 - lines collienar, overlap, value=overlap segment
     99#                      3 - lines collinear, no overlap, value is None
     100#                      4 - lines parallel, value is None
    74101def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8):
    75     """Returns intersecting point between two line segments or None
    76     (if parallel or no intersection is found).
    77 
    78     However, if parallel lines coincide partly (i.e. shara a common segment,
     102    """Returns intersecting point between two line segments.
     103
     104    However, if parallel lines coincide partly (i.e. share a common segment),
    79105    the line segment where lines coincide is returned
    80    
    81106
    82107    Inputs:
    … …  
    85110                      corresponding to a point.
    86111
    87 
    88112    Output:
    89         status, value
    90 
    91         where status and value is interpreted as follows
    92        
     113        status, value - where status and value is interpreted as follows:
    93114        status == 0: no intersection, value set to None.
    94115        status == 1: intersection point found and returned in value as [x,y].
    95         status == 2: Collinear overlapping lines found. Value takes the form [[x0,y0], [x1,y1]].
     116        status == 2: Collinear overlapping lines found.
     117                     Value takes the form [[x0,y0], [x1,y1]].
    96118        status == 3: Collinear non-overlapping lines. Value set to None.
    97         status == 4: Lines are parallel with a fixed distance apart. Value set to None.
    98    
     119        status == 4: Lines are parallel. Value set to None.
    99120    """
    100121
    … …  
    102123
    103124    line0 = ensure_numeric(line0, num.float)
    104     line1 = ensure_numeric(line1, num.float)   
     125    line1 = ensure_numeric(line1, num.float)
    105126
    106127    x0 = line0[0,0]; y0 = line0[0,1]
    … …  
    113134    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
    114135    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
    115        
     136
    116137    if num.allclose(denom, 0.0, rtol=rtol, atol=atol):
    117138        # Lines are parallel - check if they are collinear
    … …  
    123144                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
    124145
    125             return collinear_result[state_tuple]([x0,y0],[x1,y1],[x2,y2],[x3,y3])
     146            return collinear_result[state_tuple]([x0,y0], [x1,y1],
     147                                                 [x2,y2], [x3,y3])
    126148        else:
    127149            # Lines are parallel but aren't collinear
    128             return 4, None #FIXME (Ole): Add distance here instead of None 
     150            return 4, None #FIXME (Ole): Add distance here instead of None
    129151    else:
    130152        # Lines are not parallel, check if they intersect
    131153        u0 = u0/denom
    132         u1 = u1/denom       
     154        u1 = u1/denom
    133155
    134156        x = x0 + u0*(x1-x0)
    … …  
    137159        # Sanity check - can be removed to speed up if needed
    138160        assert num.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol)
    139         assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)       
     161        assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)
    140162
    141163        # Check if point found lies within given line segments
    142         if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 
     164        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
    143165            # We have intersection
    144166            return 1, num.array([x, y])
    … …  
    147169            return 0, None
    148170
    149 
     171##
     172# @brief Finds intersection point of two line segments.
     173# @param line0 First line ((x1,y1), (x2,y2)).
     174# @param line1 Second line ((x1,y1), (x2,y2)).
     175# @return (status, value) where:
     176#             status = 0 - no intersection, value set to None
     177#                      1 - intersection found, value=(x,y)
     178#                      2 - lines collienar, overlap, value=overlap segment
     179#                      3 - lines collinear, no overlap, value is None
     180#                      4 - lines parallel, value is None
     181# @note Wrapper for C function.
    150182def NEW_C_intersection(line0, line1):
    151     #FIXME(Ole): To write in C
    152     """Returns intersecting point between two line segments or None
    153     (if parallel or no intersection is found).
    154 
    155     However, if parallel lines coincide partly (i.e. shara a common segment,
     183    """Returns intersecting point between two line segments.
     184
     185    However, if parallel lines coincide partly (i.e. share a common segment),
    156186    the line segment where lines coincide is returned
    157    
    158187
    159188    Inputs:
    160189        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
    161                       A line can also be a 2x2 numeric array with each row
     190                      A line can also be a 2x2 numpy array with each row
    162191                      corresponding to a point.
    163192
    164 
    165193    Output:
    166         status, value
    167 
    168         where status is interpreted as follows
    169        
    170         status == 0: no intersection with value set to None
    171         status == 1: One intersection point found and returned in value as [x,y]
    172         status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
    173         status == 3: Lines would coincide but only if extended. Value set to None
    174         status == 4: Lines are parallel with a fixed distance apart. Value set to None.
    175    
    176     """
    177 
     194        status, value - where status and value is interpreted as follows:
     195        status == 0: no intersection, value set to None.
     196        status == 1: intersection point found and returned in value as [x,y].
     197        status == 2: Collinear overlapping lines found.
     198                     Value takes the form [[x0,y0], [x1,y1]].
     199        status == 3: Collinear non-overlapping lines. Value set to None.
     200        status == 4: Lines are parallel. Value set to None.
     201    """
    178202
    179203    line0 = ensure_numeric(line0, num.float)
    180     line1 = ensure_numeric(line1, num.float)   
     204    line1 = ensure_numeric(line1, num.float)
    181205
    182206    status, value = _intersection(line0[0,0], line0[0,1],
    … …  
    187211    return status, value
    188212
    189 
    190 
    191 
     213##
     214# @brief Determine if one point is inside a polygon.
     215# @param point The point of interest.
     216# @param polygon The polygon to test inclusion in.
     217# @param closed True if points on boundary are considered 'inside' polygon.
     218# @param verbose True if this function is to be verbose.
     219# @return True if point is inside the polygon.
     220# @note Uses inside_polygon() to do the work.
     221# @note Raises Exception if more than one point supplied.
    192222def is_inside_polygon(point, polygon, closed=True, verbose=False):
    193223    """Determine if one point is inside a polygon
    … …  
    195225    See inside_polygon for more details
    196226    """
    197 
    198 #    print 'polygon.py: 198, is_inside_polygon: point=%s' % str(point)
    199 #    print 'polygon.py: 198, is_inside_polygon: polygon=%s' % str(polygon)
    200227
    201228    indices = inside_polygon(point, polygon, closed, verbose)
    … …  
    208235        msg = 'is_inside_polygon must be invoked with one point only'
    209236        raise Exception, msg
    210    
    211 
     237
     238##
     239# @brief Determine which of a set of points are inside a polygon.
     240# @param points A set of points (tuple, list or array).
     241# @param polygon A set of points defining a polygon (tuple, list or array).
     242# @param closed True if points on boundary are considered 'inside' polygon.
     243# @param verbose True if this function is to be verbose.
     244# @return A list of indices of points inside the polygon.
    212245def inside_polygon(points, polygon, closed=True, verbose=False):
    213246    """Determine points inside a polygon
    214247
    215248       Functions inside_polygon and outside_polygon have been defined in
    216        terms af separate_by_polygon which will put all inside indices in
     249       terms of separate_by_polygon which will put all inside indices in
    217250       the first part of the indices array and outside indices in the last
    218251
    … …  
    222255       a list or a numeric array
    223256    """
    224 
    225     #if verbose: print 'Checking input to inside_polygon'
    226257
    227258    try:
    … …  
    232263        # If this fails it is going to be because the points can't be
    233264        # converted to a numeric array.
    234         msg = 'Points could not be converted to numeric array' 
     265        msg = 'Points could not be converted to numeric array'
    235266        raise Exception, msg
    236267
    … …  
    242273        # If this fails it is going to be because the points can't be
    243274        # converted to a numeric array.
    244         msg = 'Polygon %s could not be converted to numeric array' %(str(polygon))
     275        msg = ('Polygon %s could not be converted to numeric array'
     276               % (str(polygon)))
    245277        raise Exception, msg
    246278
    … …  
    253285                                                verbose=verbose)
    254286
    255 #    print 'polygon.py: 255, inside_polygon: indices=%s' % str(indices)
    256 #    print 'polygon.py: 255, inside_polygon: count=%s' % str(count)
    257287    # Return indices of points inside polygon
    258288    return indices[:count]
    259289
    260 
    261 
     290##
     291# @brief Determine if one point is outside a polygon.
     292# @param point The point of interest.
     293# @param polygon The polygon to test inclusion in.
     294# @param closed True if points on boundary are considered 'inside' polygon.
     295# @param verbose True if this function is to be verbose.
     296# @return True if point is outside the polygon.
     297# @note Uses inside_polygon() to do the work.
    262298def is_outside_polygon(point, polygon, closed=True, verbose=False,
    263299                       points_geo_ref=None, polygon_geo_ref=None):
    … …  
    268304
    269305    indices = outside_polygon(point, polygon, closed, verbose)
    270                               #points_geo_ref, polygon_geo_ref)
    271306
    272307    if indices.shape[0] == 1:
    … …  
    277312        msg = 'is_outside_polygon must be invoked with one point only'
    278313        raise Exception, msg
    279    
    280 
     314
     315##
     316# @brief Determine which of a set of points are outside a polygon.
     317# @param points A set of points (tuple, list or array).
     318# @param polygon A set of points defining a polygon (tuple, list or array).
     319# @param closed True if points on boundary are considered 'inside' polygon.
     320# @param verbose True if this function is to be verbose.
     321# @return A list of indices of points outside the polygon.
    281322def outside_polygon(points, polygon, closed = True, verbose = False):
    282323    """Determine points outside a polygon
    283324
    284325       Functions inside_polygon and outside_polygon have been defined in
    285        terms af separate_by_polygon which will put all inside indices in
     326       terms of separate_by_polygon which will put all inside indices in
    286327       the first part of the indices array and outside indices in the last
    287328
    … …  
    289330    """
    290331
    291     #if verbose: print 'Checking input to outside_polygon'
    292332    try:
    293333        points = ensure_numeric(points, num.float)
    … …  
    306346        raise Exception, msg
    307347
    308 
    309348    if len(points.shape) == 1:
    310349        # Only one point was passed in. Convert to array of points
    … …  
    321360    else:
    322361        return indices[count:][::-1]  #return reversed
    323        
    324 
    325 def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
     362
     363##
     364# @brief Separate a list of points into two sets inside+outside a polygon.
     365# @param points A set of points (tuple, list or array).
     366# @param polygon A set of points defining a polygon (tuple, list or array).
     367# @param closed True if points on boundary are considered 'inside' polygon.
     368# @param verbose True if this function is to be verbose.
     369# @return A tuple (in, out) of point indices for poinst inside amd outside.
     370def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
    326371    """Determine points inside and outside a polygon
    327372
    328373