Ignore:
Timestamp:
Jun 18, 2010, 5:49:57 PM (12 years ago)
Author:
steve
Message:

Continuing to numpy the for loops

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py

    r7852 r7860  
    1010
    1111
    12 def gravity(domain):
     12def gravity_for_loops(domain):
    1313    """Apply gravitational pull in the presence of bed slope
    1414    """
     
    1717
    1818    xmom  = domain.quantities['xmomentum'].explicit_update
    19     stage = domain.quantities['stage'].explicit_update
    20 #    ymom = domain.quantities['ymomentum'].explicit_update
    2119
    2220    Stage = domain.quantities['stage']
    2321    Elevation = domain.quantities['elevation']
    24     #h = Stage.edge_values - Elevation.edge_values
     22
    2523    h = Stage.vertex_values - Elevation.vertex_values
    2624    b = Elevation.vertex_values
     
    3129
    3230    for k in range(domain.number_of_elements):
    33 #        avg_h = sum( h[k,:] )/3
    3431        avg_h = sum( h[k,:] )/2
    3532
     
    4037        b0, b1 = b[k,:]
    4138
    42         w0, w1 = w[k,:]
    43         wx = gradient(x0, x1, w0, w1)
    4439
    4540        #zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
     
    5146        #stage[k] = 0.0
    5247
     48
     49def gravity(domain):
     50    """Apply gravitational pull in the presence of bed slope
     51    """
     52
     53    xmom  = domain.quantities['xmomentum'].explicit_update
     54
     55    Elevation = domain.quantities['elevation']
     56    Height    = domain.quantities['height']
     57
     58    hc = Height.centroid_values
     59    b  = Elevation.vertex_values
     60
     61    x = domain.vertices
     62    g = domain.g
     63
     64    x0 = x[:,0]
     65    x1 = x[:,1]
     66
     67    b0 = b[:,0]
     68    b1 = b[:,1]
     69
     70    bx = (b1-b0)/(x1-x0)
     71
     72    xmom[:] = xmom - g*bx*hc
    5373
    5474def manning_friction(domain):
Note: See TracChangeset for help on using the changeset viewer.