From 396f0573518c9cd066c31c1b59505489dd0529d4 Mon Sep 17 00:00:00 2001
From: MatthieuJacobs <matthieu.jacobs@student.kuleuven.be>
Date: Mon, 7 Mar 2022 16:39:56 +0100
Subject: [PATCH] Files from Philip for grid quality assessment
---
GridQuality/grid_analyze.py | 203 ++++++++++++++++++++++++++++++++++++
GridQuality/quick_g.py | 14 +++
2 files changed, 217 insertions(+)
create mode 100644 GridQuality/grid_analyze.py
create mode 100644 GridQuality/quick_g.py
diff --git a/GridQuality/grid_analyze.py b/GridQuality/grid_analyze.py
new file mode 100644
index 0000000..31fc5fa
--- /dev/null
+++ b/GridQuality/grid_analyze.py
@@ -0,0 +1,203 @@
+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
diff --git a/GridQuality/quick_g.py b/GridQuality/quick_g.py
new file mode 100644
index 0000000..5824cbd
--- /dev/null
+++ b/GridQuality/quick_g.py
@@ -0,0 +1,14 @@
+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()
--
GitLab