Changeset 8398

Timestamp:

Apr 12, 2012, 9:00:50 PM (13 years ago)

Author:

davies

Message:

New util

File:

• trunk/anuga_work/development/gareth/tests/util.py (modified) (2 diffs)

trunk/anuga_work/development/gareth/tests/util.py

@@ -16 +16 @@
                                them according to their 'time'. This means that
                                they can be stuck together using 'combine_outputs' correctly
- 
+
+    util.triangle_areas -- compute the areas of every triangle in a get_outputs object [ must be vertex-based]
+
+    util.water_volume -- compute the water volume at every time step in an sww file (needs both vertex and centroid value input). 
 
     Here is an example ipython session which uses some of these functions:
@@ -372 +375 @@
 
     return near_points, local_coord
+
+########################
+# TRIANGLE AREAS, WATER VOLUME
+def triangle_areas(p, subset='null'):
+    # Compute areas of triangles in p -- assumes p contains vertex information
+    # subset = vector of centroid indices to include in the computation. 
+
+    if(subset=='null'):
+        subset=range(len(p.vols[:,0]))
+    
+    x0=p.x[p.vols[subset,0]]
+    x1=p.x[p.vols[subset,1]]
+    x2=p.x[p.vols[subset,2]]
+    
+    y0=p.y[p.vols[subset,0]]
+    y1=p.y[p.vols[subset,1]]
+    y2=p.y[p.vols[subset,2]]
+    
+    # Vectors for cross-product
+    v1_x=x0-x1
+    v1_y=y0-y1
+    #
+    v2_x=x2-x1
+    v2_y=y2-y1
+    # Area
+    area=(v1_x*v2_y-v1_y*v2_x)*0.5
+    area=abs(area)
+    return area
+
+###
+
+def water_volume(p,p2, per_unit_area=False, subset='null'):
+    # Compute the water volume from p(vertex values) and p2(centroid values)
+
+    if(subset=='null'):
+        subset=range(len(p2.x))
+
+    l=len(p2.time)
+    area=triangle_areas(p, subset=subset)
+    
+    total_area=area.sum()
+    volume=p2.time*0.
+   
+    # This accounts for how volume is measured in ANUGA 
+    # Compute in 2 steps to reduce precision error (important sometimes)
+    for i in range(l):
+        #volume[i]=((p2.stage[i,subset]-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+        volume[i]=((p2.stage[i,subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+        volume[i]=volume[i]((-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
+    
+    if(per_unit_area):
+        volume=volume/total_area 
+    
+    return volume
+
+
+