Files from Philip for grid quality assessment

GridQuality/grid_analyze.py

+import numpy as np
+from matplotlib import pyplot as plt
+from combine.quick_g import quick_g_plot
+
+# not whole loop, only tenth -> phi doesnt loop    
+def yield_grid_cells(g, r_lahead=1, tht_lahead=1, phi_lahead=0):
+    g = g[:,:-1] # remove double tht=0,2pi values
+    # implement general looping edges g[r,tht,phi,c]; tht loops, r, phi not
+    ors, othts, ophis, _ = g.shape  # original xxx size
+    g = np.pad(g, [(0,0), (0,tht_lahead), (0,0), (0,0)], mode='wrap')
+    
+    # create smaller output array
+    h = np.empty((ors-r_lahead, othts, ophis-phi_lahead,
+                  1+r_lahead, tht_lahead+1, phi_lahead+1, 3))
+    
+    # iterate over all parameter coordinates
+    for cr in range(ors-r_lahead):
+        for ctht in range(othts): #tht loops
+            for cphi in range(ophis-phi_lahead):
+                h[cr, ctht, cphi] = g[cr:cr+r_lahead+1,
+                                      ctht:ctht+tht_lahead+1,
+                                      cphi:cphi+phi_lahead+1]    
+    return h
+     
+           
+def volume(g):
+    wv = yield_grid_cells(g, 1, 1, 1)  # window view
+    # calculate coord differences 
+    r = wv[...,0,0,0,0]
+    dr = np.diff(wv[...,:,0,0,:], axis=-2)[...,0,:]
+    dz = np.diff(wv[...,0,:,0,:], axis=-2)[...,0,:]
+    drphi = np.diff(wv[...,0,0,:,:], axis=-2)[...,0,:]
+    # change phi angle difference to space difference by multiplying with r
+    drphi *= r.reshape(r.shape + (1,)) # dphi * r
+    # calculate volume from column vector determinant
+    vectors = np.empty(dr.shape + (3,))
+    vectors[...,0], vectors[...,1], vectors[...,2] = dr, dz, drphi
+    vols = np.linalg.det(vectors)
+    return vols
+        
+def non_orthogonality(g):
+    wv = yield_grid_cells(g, 2, 2, 0)  # window view
+    # the two directions of cell edges, along r and along tht
+    rs, thts, phis = wv.shape[:3]
+    # angles on r and tht faces
+    a_r, a_tht = np.empty((rs, thts, phis)), np.empty((rs, thts, phis))
+    
+    # using numpy vector functions is fast, but takes ages to code. will loop
+    for r in range(rs):
+        for tht in range(thts):
+            for phi in range(phis):
+                x = wv[r,tht,phi]
+                # cell centers
+                mid = np.mean(x[0:2,0:2], axis=(0,1))[0]
+                mid_ar = np.mean(x[1:3,0:2], axis=(0,1))[0]
+                mid_atht = np.mean(x[0:2,1:3], axis=(0,1))[0]
+                # direction of cell midpoint connections 
+                v_ar = (mid_ar - mid)/np.linalg.norm(mid_ar - mid)
+                v_atht = (mid_atht - mid)/np.linalg.norm(mid_atht - mid)
+                # direction of edges 
+                t_ar = (x[1,1]-x[1,0])/np.linalg.norm(x[1,1]-x[1,0])
+                t_atht = (x[1,1]-x[0,1])/np.linalg.norm(x[1,1]-x[0,1])
+                # angle between edge direction and midpoint direction
+                a_r[r,tht,phi] = np.abs(np.arcsin(np.dot(t_ar, v_ar)))
+                a_tht[r,tht,phi] = np.abs(np.arcsin(np.dot(t_atht, v_atht)))
+                
+    a_r *= 180/np.pi; a_tht *= 180/np.pi
+    return a_r, a_tht
+    
+def unevenness(g):
+    wv = yield_grid_cells(g, 2, 2, 0)  # window view
+    # the two directions of cell edges, along r and along tht
+    rs, thts, phis = wv.shape[:3]
+    # unevenness on r and tht faces, also named a
+    a_r, a_tht = np.empty((rs, thts, phis)), np.empty((rs, thts, phis))
+    
+    # using numpy vector functions is fast, but takes ages to code. will loop
+    for r in range(rs):
+        for tht in range(thts):
+            for phi in range(phis):
+                x = wv[r,tht,phi,...,:2]
+                # cell centers
+                mid = np.mean(x[0:2,0:2], axis=(0,1))[0]
+                mid_ar = np.mean(x[1:3,0:2], axis=(0,1))[0]
+                mid_atht = np.mean(x[0:2,1:3], axis=(0,1))[0]
+                # direction of cell midpoint connections 
+                v_ar = (mid_ar - mid)/np.linalg.norm(mid_ar - mid)
+                v_atht = (mid_atht - mid)/np.linalg.norm(mid_atht - mid)
+                # cell midpoint line midpoint m_f
+                linemid_ar = np.mean([mid, mid_ar], axis=0)
+                linemid_atht = np.mean([mid, mid_atht], axis=0)
+                # cell corner connection line (edge) midpoint c_f
+                facemid_ar = np.mean(x[1,0:2], axis=0)[0]
+                facemid_atht = np.mean(x[0:2,1], axis=0)[0]
+                # line distance to point
+                ld_ar = np.linalg.norm(np.cross(v_ar, facemid_ar-mid_ar))
+                ld_atht = np.linalg.norm(np.cross(v_atht, facemid_atht-mid_atht))
+                # rotated cell mitpoint dir by 90 deg
+                n_ar = np.array([-v_ar[1], v_ar[0]])
+                n_atht = np.array([-v_atht[1], v_atht[0]])
+                # check if needed to add or remove
+                f_ar = np.sign(np.dot(n_ar, facemid_ar-mid_ar))
+                f_atht = np.sign(np.dot(n_atht, facemid_atht-mid_atht))
+                # move corner connection line (edge) midpoint on line c_f'
+                proj_facemid_ar = facemid_ar - f_ar*ld_ar*n_ar
+                proj_facemid_atht = facemid_atht - f_atht*ld_atht*n_atht
+                # distance of c_f' and m_f
+                num_dist_ar = np.linalg.norm(proj_facemid_ar-linemid_ar)
+                num_dist_atht = np.linalg.norm(proj_facemid_atht-linemid_atht)
+                # distance between two cell midpoints 
+                den_dist_ar = np.linalg.norm(mid_ar - mid)
+                den_dist_atht = np.linalg.norm(mid_atht - mid)
+                # measure of unevenness
+                a_r[r,tht,phi] = num_dist_ar/den_dist_ar
+                a_tht[r,tht,phi] = num_dist_atht/den_dist_atht
+                
+    return a_r, a_tht
+
+def skewness(g):
+    wv = yield_grid_cells(g, 2, 2, 0)  # window view
+    # the two directions of cell edges, along r and along tht
+    rs, thts, phis = wv.shape[:3]
+    # skewness on r and tht faces, also named a
+    a_r, a_tht = np.empty((rs, thts, phis)), np.empty((rs, thts, phis))
+    
+    # using numpy vector functions is fast, but takes ages to code. will loop
+    for r in range(rs):
+        for tht in range(thts):
+            for phi in range(phis):
+                x = wv[r,tht,phi,...,:2]
+                # cell centers
+                mid = np.mean(x[0:2,0:2], axis=(0,1))[0]
+                mid_ar = np.mean(x[1:3,0:2], axis=(0,1))[0]
+                mid_atht = np.mean(x[0:2,1:3], axis=(0,1))[0]
+                # direction of cell midpoint connections 
+                v_ar = (mid_ar - mid)/np.linalg.norm(mid_ar - mid)
+                v_atht = (mid_atht - mid)/np.linalg.norm(mid_atht - mid)
+                # cell corner connection line (edge) midpoint c_f
+                facemid_ar = np.mean(x[1,0:2], axis=0)[0]
+                facemid_atht = np.mean(x[0:2,1], axis=0)[0]
+                # line distance to point == ||c_f' - c_f||
+                ld_ar = np.abs(np.cross(v_ar, facemid_ar-mid_ar))
+                ld_atht = np.abs(np.cross(v_atht, facemid_atht-mid_atht))
+                # distance between two cell midpoints 
+                den_dist_ar = np.linalg.norm(mid_ar - mid)
+                den_dist_atht = np.linalg.norm(mid_atht - mid) 
+                # measure of skewness
+                a_r[r,tht,phi] = ld_ar/den_dist_ar
+                a_tht[r,tht,phi] = ld_atht/den_dist_atht
+                
+    return a_r, a_tht
+
+def nonconvex(g):
+    wv = yield_grid_cells(g, 1, 1, 0)  # window view
+    # the two directions of cell edges, along r and along tht
+    rs, thts, phis = wv.shape[:3]
+    # skewness on r and tht faces, also named a
+    is_convex = np.empty((rs, thts, phis))
+    
+    # using numpy vector functions is fast, but takes ages to code. will loop
+    for r in range(rs):
+        for tht in range(thts):
+            for phi in range(phis):
+                x = wv[r,tht,phi,...,0,:2]
+                # from point, to point, *other point indices to check against
+                # its just looping through the subarray, feel free to improve
+                indices = [[(0,0), (1,0), (1,1), (0,1)],
+                           [(1,0), (1,1), (0,1), (0,0)],
+                           [(1,1), (0,1), (0,0), (1,0)],
+                           [(0,1), (0,0), (1,0), (1,1)]]
+                for (of, to, A, B) in indices:
+                    di = x[to] - x[of] # direction of line
+                    # check if A and B on same side of line. If so, then 
+                    # the crossproduct with the direction has the same sign
+                    vA = np.sign(np.cross(di, x[A]-x[of]))
+                    vB = np.sign(np.cross(di, x[B]-x[of]))
+                    if vA == vB:
+                        is_convex[r,tht,phi] = vA
+                    else:
+                        is_convex[r,tht,phi] = 0
+    
+    return is_convex
+    
+
+if __name__ == "__main__":
+    try:
+        g = np.load("g_possible.npy")
+        # g = np.load("g_kaputt.npy")
+    except:
+        print("Need to generate grid array 'g.npy' with another program")
+        raise
+    
+    
+    #nonc = nonconvex(g) # TODO  1 == good, 0 == nonconvex, -1 == weird
+    
+    # volumes = volume(g)
+    # angles_r, angles_tht = non_orthogonality(g)
+    # uneven_r, uneven_tht = unevenness(g)
+    # skewness_r, skewness_tht = skewness(g)
+    # print(f"{np.sum(volumes):.3E} <- Whole volume in m^3")
+    # print(f"{np.min(volumes):.3E}, {np.max(volumes):.3E} <- Min, Max cell volume in m^3")
+    # print(f"{np.mean(volumes):.3E}, {np.std(volumes):.3E} <- Mean, Std of cell volume in m^3")
+    
\ No newline at end of file

GridQuality/quick_g.py

+import numpy as np
+from matplotlib import pyplot as plt
+
+def quick_g_plot(g, phi=-1, title=None):
+    R, z, _ = np.transpose(g[:,:,phi,:], (2,0,1))
+    fig, ax = plt.subplots(figsize=(8,8))
+    ax.set_aspect("equal")
+    ax.plot(R,z, c="blue", lw=0.3)
+    ax.plot(R.T,z.T, c="red", lw=0.3)
+    if title:
+        ax.set_title(title)
+    elif title is None:
+        ax.set_title(f"Crossection at the {phi}. slice")
+    fig.tight_layout()