working wse and U depth average with new code

Author: Rusty Holleman

stompy/model/openfoam/depth_average.py

@@ -926,14 +926,16 @@ def precalc_raster_weights_proc_by_faces(fld, xyz, face_nodes, face_cells, cell_
     #  compile into matrix
 
     mesh_state = (xyz,face_nodes,face_cells,cell_faces)
+    n_cells_orig = len(cell_faces)
 
     fld_x,fld_y = fld.xy() # I think these are pixel centers
 
+    cell_mapping=None
     for col in utils.progress(range(fld.shape[1])):
         if col==0: continue
         xmin=fld_x[col]-fld.dx/2
         xmax=fld_x[col]+fld.dx/2
-        mesh_state = mesh_ops.mesh_slice(np.r_[1,0,0], xmin, *mesh_state)
+        cell_mapping, mesh_state = mesh_ops.mesh_slice(np.r_[1,0,0], xmin, cell_mapping, *mesh_state)
 
     if 0: # debugging checks
         print("Checking cells")
@@ -946,7 +948,7 @@ def precalc_raster_weights_proc_by_faces(fld, xyz, face_nodes, face_cells, cell_
         if row==0: continue
         ymin=fld_y[row]-fld.dy/2
         ymax=fld_y[row]+fld.dy/2
-        mesh_state = mesh_ops.mesh_slice(np.r_[0,1,0], ymin, *mesh_state)
+        cell_mapping, mesh_state = mesh_ops.mesh_slice(np.r_[0,1,0], ymin, cell_mapping, *mesh_state)
 
     if 0: # debugging checks
         print("Checking cells after all slicing")
@@ -958,7 +960,7 @@ def precalc_raster_weights_proc_by_faces(fld, xyz, face_nodes, face_cells, cell_
     volumes,centers = mesh_ops.mesh_cell_volume_centers(*mesh_state)
 
     raster_weights = sparse.dok_matrix((fld.F.shape[0]*fld.F.shape[1],
-                                        centers.shape[0]),
+                                        n_cells_orig),
                                        np.float64)
 
     # check bounds, then compute clipped volume.
@@ -981,7 +983,7 @@ def precalc_raster_weights_proc_by_faces(fld, xyz, face_nodes, face_cells, cell_
 
             for cIdx in np.nonzero(sel)[0]:
                 # alpha-weighting is included per timestep
-                raster_weights[pix,cIdx] = volumes[cIdx]/Apixel
+                raster_weights[pix,cell_mapping[cIdx]] = volumes[cIdx]/Apixel
 
     return raster_weights
 

stompy/model/openfoam/mesh_ops.py

@@ -227,8 +227,11 @@ def mesh_check_cell(cIdx,verbose,mesh_state):
     return np.all(np.array(outwards)>0)
 
 
-def mesh_slice(slice_normal, slice_offset, xyz, face_nodes, face_cells, cell_faces):
+def mesh_slice(slice_normal, slice_offset, cell_mapping, xyz, face_nodes, face_cells, cell_faces):
     tol = 1e-10
+    if cell_mapping is None:
+        cell_mapping = np.arange(len(cell_faces))
+        
     # identify which faces to slice:
     if 1:
         offset_xyz = np.dot(xyz,slice_normal) - slice_offset
@@ -495,7 +498,10 @@ def mesh_slice(slice_normal, slice_offset, xyz, face_nodes, face_cells, cell_fac
                         # bad here.
                         if f_i==0 and side_xyz[f_nodes[f_i+1]]!=0:
                             # last-to-first edge is the one.
-                            assert f_nodes[-1]>=n_node_orig,"More trouble with slice nodes"
+                            # pretty sure this overly restrictive. face could have come in with a node
+                            # on the slice. Even if it's a new face, one of the nodes could have already
+                            # been on the slice.
+                            #assert f_nodes[-1]>=n_node_orig,"More trouble with slice nodes"
                             nodeA = f_nodes[-1]
                             nodeB = f_nodes[0]
                             assert nodeA!=nodeB
@@ -563,6 +569,7 @@ def mesh_slice(slice_normal, slice_offset, xyz, face_nodes, face_cells, cell_fac
         new_cIdx = len(cell_faces)
         cell_faces[cIdx,:] = cell_face_neg
         cell_faces = utils.array_append(cell_faces,cell_face_pos)
+        cell_mapping = utils.array_append(cell_mapping, cell_mapping[cIdx])
         face_cells = utils.array_append(face_cells,np.array([cIdx,new_cIdx]))
 
         #import pdb
@@ -602,7 +609,7 @@ def mesh_slice(slice_normal, slice_offset, xyz, face_nodes, face_cells, cell_fac
                 else:
                     assert False,"Failed to find face to flip"
 
-    return (xyz, face_nodes, face_cells, cell_faces)
+    return cell_mapping, (xyz, face_nodes, face_cells, cell_faces)
 
 
 def mesh_cell_bboxes(xyz,face_nodes,face_cells,cell_faces):