Changeset 8945


Ignore:
Timestamp:
Jul 17, 2013, 11:11:19 PM (11 years ago)
Author:
steve
Message:

Put in a region argument to get_integral.

Location:
trunk/anuga_core/source/anuga
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/__init__.py

    r8757 r8945  
    3434
    3535from anuga.abstract_2d_finite_volumes.quantity import Quantity
     36
     37from anuga.operators.region import Region
    3638
    3739from anuga.abstract_2d_finite_volumes.util import file_function, \

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r8944 r8945  
    10651065
    10661066        Usage:
    1067             i = get_extreme_index()
     1067            i = Q.get_extreme_index()
    10681068
    10691069        Notes:
    … …  
    15201520        self.y_gradient[:] = 0.0
    15211521
    1522     def get_integral(self, full_only=True):
    1523 
    1524         """Compute the integral of quantity across entire domain."""
     1522    def get_integral(self, full_only=True, region=None, indices=None):
     1523
     1524        """Compute the integral of quantity across entire domain,
     1525        or over a region. Eg
     1526        my_region = anuga.Region(polygon = user_polygon)
     1527
     1528        Q.get_integral(region = my_region)"""
     1529
     1530        assert region is None or indices is None
    15251531
    15261532       
    15271533        areas = self.domain.get_areas()
    15281534
    1529         if full_only:
    1530             indices = num.where(self.domain.tri_full_flag ==1)[0]
     1535        if region is None and indices is None:
     1536            if full_only:
     1537                indices = num.where(self.domain.tri_full_flag ==1)[0]
     1538                return num.sum(areas[indices]*self.centroid_values[indices])
     1539            else:
     1540                return num.sum(areas*self.centroid_values)
     1541
     1542        if indices is None:
     1543            indices = region.get_indices(full_only)
     1544
     1545        if indices == []:
     1546            return 0.0
     1547        elif indices is None:
     1548            return num.sum(areas*self.centroid_values)
     1549        else:
    15311550            return num.sum(areas[indices]*self.centroid_values[indices])
    1532         else:
    1533             return num.sum(areas*self.centroid_values)
    15341551
    15351552

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r8930 r8945  
    425425        assert num.allclose (quantity.get_integral(), ref_integral)
    426426
    427 
     427    def test_integral_with_region(self):
     428        quantity = Quantity(self.mesh4)
     429
     430        # Try constants first
     431        const = 5
     432        quantity.set_values(const, location = 'vertices')
     433        #print 'Q', quantity.get_integral()
     434
     435        assert num.allclose(quantity.get_integral(), self.mesh4.get_area() * const)
     436
     437        # Try with a linear function
     438        def f(x, y):
     439            return x+y
     440
     441        quantity.set_values(f, location = 'vertices')
     442
     443
     444        from anuga import Region
     445
     446        reg1 = Region(self.mesh4, indices=[2])
     447        ref_integral = (10.0/3) * 2
     448
     449        assert num.allclose (quantity.get_integral(region=reg1), ref_integral)
     450
     451        reg2 = Region(self.mesh4, indices=[2,3])
     452        ref_integral = (10.0/3 + 10.0/3) * 2
     453
     454        assert num.allclose (quantity.get_integral(region=reg2), ref_integral)
     455
     456        id=[2,3]
     457        ref_integral = (10.0/3 + 10.0/3) * 2
     458
     459        assert num.allclose (quantity.get_integral(indices=id), ref_integral)
    428460
    429461    def test_set_vertex_values(self):
    … …  
    25402572if __name__ == "__main__":
    25412573    suite = unittest.makeSuite(Test_Quantity, 'test')
    2542     runner = unittest.TextTestRunner(verbosity=2)
     2574    runner = unittest.TextTestRunner(verbosity=1)
    25432575    runner.run(suite)

  • trunk/anuga_core/source/anuga/operators/erosion_operators.py

    r8944 r8945  
    1111
    1212import numpy as num
    13 from anuga.geometry.polygon import inside_polygon
    14 from anuga import indent
     13
    1514
    1615from anuga import Domain

  • trunk/anuga_core/source/anuga/operators/rate_operators.py

    r8944 r8945  
    1313from anuga import indent
    1414import numpy as num
    15 from warnings import warn
    1615import anuga.utilities.log as log
    17 from anuga.geometry.polygon import inside_polygon
    18 
    19 
    20 from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    21                                               Modeltime_too_late
    2216from anuga.utilities.function_utils import evaluate_temporal_function
    2317
    24 from anuga import Domain
    25 from anuga import Quantity
     18
    2619from anuga.operators.base_operator import Operator
    2720from anuga.operators.region import Region
    … …  
    10295
    10396        if self.rate_spatial:
    104             if self.indices is None:
     97            if indices is None:
    10598                x = self.coord_c[:,0]
    10699                y = self.coord_c[:,1]
    107100            else:
    108                 x = self.coord_c[self.indices,0]
    109                 y = self.coord_c[self.indices,1]
     101                x = self.coord_c[indices,0]
     102                y = self.coord_c[indices,1]
     103
    110104            rate = self.get_spatial_rate(x,y,t)
    111105        else:
    … …  
    117111
    118112        if num.all(rate >= 0.0):
    119             if self.indices is None:
     113            if indices is None:
    120114                self.stage_c[:] = self.stage_c[:]  \
    121115                       + factor*rate*timestep
    … …  
    124118                       + factor*rate*timestep
    125119        else: # Be more careful if rate < 0
    126             if self.indices is None:
     120            if indices is None:
    127121                self.stage_c[:] = num.maximun(self.stage_c  \
    128122                       + factor*rate*timestep, self.elev_c )

  • trunk/anuga_core/source/anuga/operators/region.py

    r8938 r8945  
    1010from anuga import Quantity
    1111import numpy as num
    12 import anuga.utilities.log as log
     12
    1313
    1414from anuga.geometry.polygon import inside_polygon
    … …  
    1616from anuga.utilities.function_utils import determine_function_type
    1717
    18 from anuga import indent
     18#from anuga import indent
    1919
    2020
    … …  
    5656        if self.indices is not None:
    5757            # This overrides polygon, center and radius
    58             pass
     58            self.indices = num.asarray(self.indices)
     59
     60            if self.indices.size == 0:
     61                self.indices = []
     62
    5963        elif (self.center is not None) and (self.radius is not None):
    6064
    … …  
    7175        else:
    7276            assert self.indices is None or self.indices is []
     77
     78
     79
     80        if self.indices == []:
     81            self.full_indices = []
     82        elif self.indices is None:
     83            self.full_indices = num.where(self.domain.tri_full_flag ==1)[0]
     84        else:
     85            self.full_indices = self.indices[self.domain.tri_full_flag[self.indices]==1]
    7386
    7487
    … …  
    93106                indices.append(k)
    94107
    95         self.indices = indices
     108        if indices is []:
     109            self.indices = indices
     110        else:
     111            self.indices = num.asarray(indices)
    96112
    97113        msg = 'No centroids intersect circle center'+str(c)+' radius '+str(r)
    … …  
    104120        points = self.domain.get_centroid_coordinates(absolute=True)
    105121
     122        indices = num.asarray(inside_polygon(points, self.polygon))
    106123
    107         self.indices = inside_polygon(points, self.polygon)
     124        if indices is []:
     125            self.indices = indices
     126        else:
     127            self.indices = num.asarray(indices)
    108128
    109129
     130    def get_indices(self, full_only=True):
     131
     132        if full_only:
     133            return self.full_indices
     134        else:
     135            return self.indices
     136
     137       
    110138
    111139
    112 

  • trunk/anuga_core/source/anuga/operators/run_collect_max_stage_operator.py

    r8944 r8945  
    6262Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
    6363Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
    64 Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
     64Bd = anuga.Dirichlet_boundary([1,0.,0.])    # Constant boundary values
    6565
    6666
    6767# Associate boundary tags with boundary objects
    6868domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    69 
    70 
    7169
    7270#-----------------------------------------------------------------------------
    … …  
    8583
    8684
    87 
    8885# save the max_stage centroid data to a text file
    8986max_operator.save_centroid_data_to_csv()
    9087
     88# Let's have a look at the max_stage
    9189max_operator.plot_quantity()
    9290

  • trunk/anuga_core/source/anuga/operators/run_polygon_erosion.py

    r8944 r8945  
    2727    z = -x/100
    2828
    29     N = len(x)
    30     for i in range(N):
    31         # Step
    32         if 2 < x[i] < 4:
    33             z[i] += 0.4 - 0.05*y[i]
     29    # Step
     30    id = (2 < x) & (x < 4)
     31    z[id] += 0.4 - 0.05*y[id]
    3432
    35         # Permanent pole
    36         if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
    37             z[i] += 1
     33    # Permanent pole
     34    id = (x - 8)**2 + (y - 2)**2 < 0.4**2
     35    z[id] += 1
    3836
    39         # Dam
    40         if 12 < x[i] < 13:
    41             z[i] += 0.4
    42 
     37    # Dam
     38    id = (12 < x) & (x < 13)
     39    z[id] += 0.4
    4340
    4441    return z
    … …  
    7471        t = self.get_time()
    7572        dt = self.get_timestep()
    76 
    7773
    7874        updated = True

  • trunk/anuga_core/source/anuga/operators/run_set_depth_friction.py

    r8480 r8945  
    1212from anuga import Dirichlet_boundary
    1313from anuga import Time_boundary
     14from anuga import Region
     15from anuga import indent
    1416
    1517#------------------------------------------------------------------------------
    … …  
    2325                                               len1=length, len2=width)
    2426domain = Domain(points, vertices, boundary)
    25 domain.set_name('set_depth_friction') # Output name
     27domain.set_name() # Output name
    2628print domain.statistics()
    2729
    … …  
    3537    z = -x/100
    3638
    37     N = len(x)
    38     for i in range(N):
    39         # Step
    40         if 2 < x[i] < 4:
    41             z[i] += 0.4 - 0.05*y[i]
     39    # Step
     40    id = ( 2 < x ) & (x < 4)
     41    z[id] +=  0.4 - 0.05*y[id]
    4242
    43         # Permanent pole
    44         if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
    45             z[i] += 1
    46            
    47 #        # Pole 2
    48 #        if (x[i] - 14)**2 + (y[i] - 3.5)**2 < 0.4**2:
    49 #            z[i] += 1.0
     43    # Permanent pole
     44    id = (x - 8)**2 + (y - 2)**2 < 0.4**2
     45    z[id] += 1
     46
     47    # Pole 2
     48    #id =  (x - 14)**2 + (y - 3.5)**2 < 0.4**2
     49    #z[id] += 1.0
     50
    5051
    5152    return z
    52 
    5353
    5454
    … …  
    6868
    6969#------------------------------------------------------------------------------
     70# Setup operators which are applied each inner step
     71#------------------------------------------------------------------------------
     72from anuga.operators.set_friction_operators import Depth_friction_operator
     73
     74op1 = Depth_friction_operator(domain)
     75
     76p1 = [ [12.0, 2.5], [13.5, 2.5], [13.5, 4.0], [12.0, 4.0] ]
     77op2 = Depth_friction_operator(domain,
     78                                  friction_max = 10,
     79                                  friction_min = 0.0,
     80                                  polygon=p1)
     81
     82# Setup region for integrating quantities
     83p2 = [ [8.0, 2.5], [9.5, 2.5], [9.5, 4.0], [8.0, 4.0] ]
     84reg = Region(domain, polygon=p2)
     85
     86# Some useful aliases
     87stage = domain.quantities['stage']
     88elev = domain.quantities['elevation']
     89
     90#------------------------------------------------------------------------------
    7091# Evolve system through time
    7192#------------------------------------------------------------------------------
    72 
    73 from anuga.operators.set_friction_operators import Set_depth_friction_operator
    74 from anuga.operators.set_friction_operators import Polygonal_depth_friction_operator
    75 
    76 #op1 = Set_depth_friction_operator(domain)
    77 
    78 polygon1 = [ [12.0, 2.5], [13.5, 2.5], [13.5, 4.0], [12.0, 4.0] ]
    79 op2 = Polygonal_depth_friction_operator(domain, friction_max = 10000, friction_min = 0.0, polygon=polygon1)
    80 
    81 
    82 
    83 for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
     93for t in domain.evolve(yieldstep=0.1, finaltime=10.0):
    8494    domain.print_timestepping_statistics()
    8595    domain.print_operator_timestepping_statistics()
     96
     97    # let calculate the integral of height over a region
     98    height = stage-elev
     99    print indent+'Int_p2(h) = '+str(height.get_integral(region=reg))
    86100
    87101

  • trunk/anuga_core/source/anuga/operators/set_elevation_operator.py

    r8929 r8945  
    99
    1010
    11 from anuga import Domain
    12 from anuga import Quantity
    13 import numpy as num
    14 import anuga.utilities.log as log
    15 
    16 from anuga.geometry.polygon import inside_polygon
    1711
    1812from anuga.operators.base_operator import Operator

  • trunk/anuga_core/source/anuga/operators/set_friction_operators.py

    r8480 r8945  
    1010
    1111
    12 from anuga import Domain
    13 from anuga import Quantity
    14 import numpy as num
    15 import anuga.utilities.log as log
    1612
    17 from anuga.geometry.polygon import inside_polygon
    1813
    1914from anuga.operators.base_operator import Operator
    20 from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    21                                               Modeltime_too_late
     15from anuga.operators.region import Region
     16
     17
    2218from anuga import indent
    2319
    … …  
    2622default_friction_max = 0.035
    2723
    28 class Set_depth_friction_operator(Operator):
     24class Depth_friction_operator(Operator, Region):
    2925    """
    3026    Set the friction in a region
    31 
    32     indices: None == all triangles, Empty list [] no triangles
    33 
    34     rate can be a function of time.
    35 
    3627    """
    3728
    … …  
    4031                 friction_min=default_friction_min,
    4132                 friction_max=default_friction_max,
    42                  indices = None,
     33                 indices=None,
     34                 polygon=None,
     35                 center=None,
     36                 radius=None,
    4337                 description = None,
    4438                 label = None,
    … …  
    4842
    4943        Operator.__init__(self, domain, description, label, logging, verbose)
     44
     45        Region.__init__(self, domain,
     46                indices=indices,
     47                polygon=polygon,
     48                center=center,
     49                radius=radius,
     50                verbose=verbose)
    5051
    5152        #------------------------------------------
    … …  
    5960        self.friction_max = friction_max
    6061        self.friction_c = self.domain.get_quantity('friction').centroid_values
    61         self.indices = indices
     62
    6263
    6364
    … …  
    106107
    107108#===============================================================================
    108 # Specific Bed Operators for circular region.
     109# Specific Operator for circular region.
    109110#===============================================================================
    110 class Circular_set_depth_friction_operator(Set_depth_friction_operator):
     111class Circular_depth_friction_operator(Depth_friction_operator):
    111112    """
    112     Set elevation over a circular region
     113    Set friction over a circular region
    113114
    114115    """
    … …  
    121122                 verbose=False):
    122123
    123         assert center is not None
    124         assert radius is not None
    125 
    126 
    127         # Determine indices in update region
    128         N = domain.get_number_of_triangles()
    129         points = domain.get_centroid_coordinates(absolute=True)
    130 
    131 
    132         indices = []
    133 
    134         c = center
    135         r = radius
    136 
    137         self.center = center
    138         self.radius = radius
    139 
    140         intersect = False
    141         for k in range(N):
    142             x, y = points[k,:]    # Centroid
    143 
    144             if ((x-c[0])**2+(y-c[1])**2) < r**2:
    145                 intersect = True
    146                 indices.append(k)
    147 
    148 
    149         msg = 'No centroids intersect circle center'+str(center)+' radius '+str(radius)
    150         assert intersect, msg
    151124
    152125
    153126
    154 
    155         # It is possible that circle doesn't intersect with mesh (as can happen
    156         # for parallel runs
    157 
    158 
    159         Set_depth_friction_operator.__init__(self,
     127        Dpth_friction_operator.__init__(self,
    160128                                    domain,
    161129                                    friction_min=friction_min,
    162130                                    friction_max=friction_max,
    163                                     indices=indices,
     131                                    center=center,
     132                                    radius=radius,
    164133                                    verbose=verbose)
    165134
    166135
    167 
    168 
    169 
    170136#===============================================================================
    171 # Specific Bed Operators for polygonal region.
     137# Specific Operator for polygonal region.
    172138#===============================================================================
    173 class Polygonal_depth_friction_operator(Set_depth_friction_operator):
     139class Polygonal_depth_friction_operator(Depth_friction_operator):
    174140    """
    175     Add water at certain rate (ms^{-1} = vol/Area/sec) over a
    176     polygonal region
    177 
    178     rate can be a function of time.
     141    Set Friction over a polygon
    179142
    180143    """
    … …  
    187150
    188151
    189         # Determine indices in update region
    190         N = domain.get_number_of_triangles()
    191         points = domain.get_centroid_coordinates(absolute=True)
    192152
    193 
    194         indices = inside_polygon(points, polygon)
    195         self.polygon = polygon
    196 
    197         # It is possible that circle doesn't intersect with mesh (as can happen
    198         # for parallel runs
    199 
    200 
    201         Set_depth_friction_operator.__init__(self,
     153        Depth_friction_operator.__init__(self,
    202154                               domain,
    203155                               friction_min=friction_min,
    204156                               friction_max=friction_max,
    205                                indices=indices,
     157                               polygon=polygon,
    206158                               verbose=verbose)
    207159

  • trunk/anuga_core/source/anuga/operators/set_quantity_operator.py

    r8931 r8945  
    99
    1010
    11 from anuga import Domain
    12 from anuga import Quantity
    13 import numpy as num
    14 import anuga.utilities.log as log
    1511
    16 from anuga.geometry.polygon import inside_polygon
    1712
    1813from anuga.operators.base_operator import Operator

  • trunk/anuga_core/source/anuga/operators/set_stage_operator.py

    r8930 r8945  
    2626class Set_stage_operator(Set_quantity_operator):
    2727    """
    28     Set the stage in a region
    29 
    30     indices: None == all triangles, Empty list [] no triangles
    31 
    32     rate can be a function of time.
    33 
     28    Set the stage over a region
    3429    """
    3530
    … …  
    9287    """
    9388    Set stage over a circular region
    94 
    9589    """
    9690
    … …  
    116110class Polygonal_set_stage_operator(Set_stage_operator):
    117111    """
    118     Add water at certain rate (ms^{-1} = vol/Area/sec) over a
    119     polygonal region
    120 
    121     rate can be a function of time.
    122 
     112    Set stage over a polygonal region
    123113    """
    124114

  • trunk/anuga_core/source/anuga/operators/test_rate_operators.py

    r8854 r8945  
    1414from anuga.file_conversion.file_conversion import timefile2netcdf
    1515from anuga.config import time_format
     16
     17from anuga.fit_interpolate.interpolate import Modeltime_too_early
     18from anuga.fit_interpolate.interpolate import Modeltime_too_late
    1619
    1720from rate_operators import *
Note: See TracChangeset for help on using the changeset viewer.