Adding calculate volumes as a task to geomodelr.

geomodelr/__main__.py

@@ -24,6 +24,7 @@
 import shared
 import unittest
 import cpp
+import utils
 from model import GeologicalModel
 
 class ParametersException(Exception):
@@ -105,6 +106,28 @@ def intersect_plane( geojson, verbose=False ):
             raise
         line = sys.stdin.readline()
 
+def calculate_volumes( geojson, params, verbose=False ):
+    model = prep_model(geojson)
+    t = {}
+
+    def query_func(p):
+        q = model.closest_topo(p)[0]
+        if q in t:
+            return t[q]
+        else:
+            l = len(t)
+            t[q] = l
+            return l
+
+    vols, elems = utils.octtree_volume_calculation(query_func, model.geojson['bbox'], params[0], params[1])
+
+    r = { v: k for k, v in t.iteritems() }
+    print "VOLUMES"
+    vols = sorted( [ ( v, k ) for k, v in vols.iteritems() ] )
+    for v, i in vols:
+        if r[i] != "AIR":
+            print "%s: %s" % ( r[i], v )
+
 def get_information(geojson, verbose):
     # Show map, cross sections, polygons, etc.
     model = GeologicalModel(json.loads(geojson.read()))
@@ -146,13 +169,19 @@ def main(args=None):
                        help = """ intersects a plane with the faults of the model and 
                                   returns the lines using the coordinate system dictated by the plane. """)
 
+    group.add_argument("-vol", "--volume", nargs=2, metavar=("INIT_REFS", "OCTREE_REFS"),
+                        default=None, type=int,
+                        help=""" approximates the raw volumes of every matherial in the model
+                                 using an octree. Parameters are the initial number of grid
+                                 refinements, and the number of octree subdivisions after that. """)
+
     parser.add_argument("-v", "--verbose", action="store_true",
                         help = """shows more information to the user.""")
 
     parser.add_argument("-p", "--profile", action="store_true",
                         help = """profiles geomodelr.""")
 
-    parser.add_argument("model", nargs="?", type=argparse.FileType('r'))
+    parser.add_argument("model", type=argparse.FileType('r'))
     args = parser.parse_args()
 
     if args.verbose:
@@ -173,6 +202,8 @@ def main(args=None):
         args.info(args.model, args.verbose)
     elif args.intersect_plane:
         args.intersect_plane(args.model, args.verbose)
+    elif args.volume:
+        calculate_volumes(args.model, args.volume)
     else:
         parser.print_help()
 

geomodelr/test.py

@@ -452,10 +452,7 @@ def query_func(p):
         self.assertGreaterEqual(ref_grid['grid'].shape[2], ref_grid['units'].shape[2])
 
         grid = utils.generate_octtree_grid(query_func, bbox, 3, 3, 3)
-        vols, elems = utils.octtree_volume_calculation(query_func, bbox, 100, 3)
-        tot = sum( vols.values() )
-        exp = (bbox[3]-bbox[0])*(bbox[4]-bbox[1])*(bbox[5]-bbox[2])
-
+        vols, elems = utils.octtree_volume_calculation(query_func, bbox, 10, 3)
 
 def main(args=None):
     unittest.main()

geomodelr/utils.py

@@ -350,7 +350,8 @@ def octtree_volume_calculation(query_func, bbox, grid_divisions, oct_refine):
     def percent_aggregated():
         tot = sum( total_volumes.values() )
         exp = (bbox[3]-bbox[0])*(bbox[4]-bbox[1])*(bbox[5]-bbox[2])
-        print tot, exp, "%s%%" % (100*tot/exp)
+
+        print "VOLUME AGGREGATED", tot, "VOLUME EXPECTED", exp, "PERCENTAGE %s%%" % (100*tot/exp)
 
     # Generate a simple grid.
     simple = generate_simple_grid(query_func, bbox, grid_divisions)
@@ -362,7 +363,7 @@ def percent_aggregated():
     points = np.zeros((pdims[0]*pdims[1]*pdims[2], 3))
     mdims = ts(pdims, (-1,-1,-1))
     elems = np.zeros((mdims[0]*mdims[1]*mdims[2], 8), dtype=int)
-    units = np.zeros((pdims[0]*pdims[1]*pdims[2],), dtype=object)
+    units = np.zeros((pdims[0]*pdims[1]*pdims[2],), dtype=int)
     for i in xrange(pdims[0]):
         for j in xrange(pdims[1]):
             for k in xrange(pdims[2]):
@@ -401,15 +402,26 @@ def should_refine_cell(pts, uns):
             return True
         return False
 
-    def add_cell_volume(pts, uns):
+    def get_volume( pts ):
         diff = pts[7]-pts[0]
-        sm = diff[0]*diff[1]*diff[2]
-        unit = uns[0]
+        return diff[0]*diff[1]*diff[2]
+
+    def add_volume( unit, sm ):
         if unit in total_volumes:
             total_volumes[unit] += sm    
         else:
             total_volumes[unit] = sm
 
+    def add_cell_volume(pts, uns):
+        sm = get_volume( pts )
+        unit = uns[0]
+        add_volume(unit, sm)
+
+    def cell_size( ):
+        pts = points[elems[0,:]]
+        dif = pts[7]-pts[0]
+        print "CELL SIZE", dif[0], dif[1], dif[2]
+
     # Find which elements should remain, wich should be removed, and reduce points and units accordingly.
     rem_idx = []
     for i in xrange(elems.shape[0]):
@@ -449,6 +461,7 @@ def reduce_points( ):
 
     print "FILTERED", elems.shape[0]
     percent_aggregated()
+    cell_size()
     for r in xrange(oct_refine):
         print "ITERATION", r
 
@@ -530,5 +543,12 @@ def reduce_points( ):
         gc.collect()
         print "ELEMENTS", elems.shape[0]
         percent_aggregated()
+        cell_size()
 
+    for i in xrange(elems.shape[0]):
+        sm = get_volume(points[elems[i,:]])
+        for ie in elems[i,:].tolist():
+            add_volume( units[ie], sm / 8.0 )
+    print "APPROXIMATED"
+    percent_aggregated()
     return ( total_volumes, { 'points': points, 'units': units, 'elems': elems } )