calculate cross-sectional area of elements and panels.

1 files changed, 21 insertions, 14 deletions

diff --git a/makeBoundary.py b/makeBoundary.py
index 6eece9d..d81a86c 100755
--- a/makeBoundary.py
+++ b/makeBoundary.py
@@ -243,22 +243,29 @@ print(r_h[panel[ind_p]][ind_q-1], r_h[panel[ind_p]][ind_q])
 velocity_panel    = panel_target_flow/A_extra
 print(velocity_channel, velocity_panel)
 
-csa = np.zeros(len(xin))             # cross-sectional area
-csa_west = 0
-csa_chan = 0
-csa_east = 0
-for index, xitem in enumerate(xin):
-    if xitem < xin[280]:
-        csa[index] = max(0, (h_extra - zin[index])*dX)
-        csa_west += csa[index]
-    elif xitem < xin[301]:
-        csa[index] = max(0, (zmax - zin[index])*dX)
-        csa_chan += csa[index]
+csa = np.zeros(len(xin))             # cross-sectional area of element
+csa_p = np.zeros(num_panels)         # cross-sectional area of panel
+for i, p in enumerate(panel):
+    if i < ind_p:
+        # panels are filled
+        area_sum = 0
+        for m in range(marker_ind[p], marker_ind[p+1]):
+            csa[m] = max(0, (zmax - zin[m])*dX)
+            area_sum += csa[m]
+        csa_p[p] = area_sum
+    elif i == ind_p:
+        # panel is part-filled
+        area_sum = 0
+        for m in range(marker_ind[p], marker_ind[p+1]):
+            csa[m] = max(0, (h_extra - zin[m])*dX)
+            area_sum += csa[m]
+        csa_p[p] = area_sum
     else:
-        csa[index] = max(0, (zmax_east - zin[index])*dX)
-        csa_east += csa[index]
+        # panel is empty
+        csa_p[p] = 0
 
-#print('csa_west = {} csa_chan = {} csa_east = {}'.format(csa_west, csa_chan, csa_east))
+
+#print('csa_p[0] = {} csa_p[1] = {}'.format(csa_p[0], csa_p[1]))
 #print('A_i_west = {} A_i = {} A_i_east = {}'.format(A_extra, A_i[-1], A_i_east[-1]))
 
 def save_bc():