intermediate

examples/stuff.py

+import numpy as np
+from scipy import ndimage
+from skimage.morphology import medial_axis
+import matplotlib.pyplot as plt
+
+
+def microstructure(l=256):
+    """
+    Synthetic binary data: binary microstructure with blobs.
+
+    Parameters
+    ----------
+
+    l: int, optional
+        linear size of the returned image
+
+    """
+    n = 5
+    x, y = np.ogrid[0:l, 0:l]
+    mask_outer = (x - l/2)**2 + (y - l/2)**2 < (l/2)**2
+    mask = np.zeros((l, l))
+    generator = np.random.RandomState(1)
+    points = l * generator.rand(2, n**2)
+    mask[(points[0]).astype(int), (points[1]).astype(int)] = 1
+    mask = ndimage.gaussian_filter(mask, sigma=l/(4.*n))
+    return mask > mask.mean()
+
+data = microstructure(l=64)
+
+# Compute the medial axis (skeleton) and the distance transform
+skel, distance = medial_axis(data, return_distance=True)
+
+# Distance to the background for pixels of the skeleton
+dist_on_skel = distance * skel
+
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
+ax1.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
+ax1.axis('off')
+ax2.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
+ax2.contour(data, [0.5], colors='w')
+ax2.axis('off')
+
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+plt.show()

ntrfc/geometry/plane.py

@@ -6,9 +6,9 @@ from ntrfc.math.vectorcalc import vecAbs, vecProjection
 
 def massflow_plane(mesh, rhoname="rho", velocityname="U"):
     normals = mesh.compute_normals(point_normals=False, cell_normals=True)["Normals"]
-    rhos = mesh[rhoname]
+    rhos = mesh.cell_data[rhoname]
     areas = mesh.compute_cell_sizes()["Area"]
-    velocities = mesh[velocityname]
+    velocities = mesh.cell_data[velocityname]
 
     massflow = np.array(
         [vecAbs(vecProjection(velocities[i], normals[i])) for i in range(mesh.number_of_cells)]) ** 2 * rhos * areas
@@ -19,7 +19,7 @@ def massflow_plane(mesh, rhoname="rho", velocityname="U"):
 def massflowave_plane(mesh, valname, rhoname="rho", velocityname="U"):
     massflow = massflow_plane(mesh, rhoname=rhoname, velocityname=velocityname)
 
-    mass_ave = np.sum(mesh[valname] * massflow) / np.sum(massflow)
+    mass_ave = np.nansum(mesh[valname] * massflow) / np.nansum(massflow)
 
     return mass_ave
 

ntrfc/turbo/cascade_case/post.py

@@ -21,8 +21,6 @@ def blade_loading_cp(
     if not figpath:
         figpath = tempfile.mkdtemp() + "/blade_loading_cp.png"
 
-    camber_length = case_instance.blade.camber_length
-
     inlet = case_instance.mesh_dict["inlet"]
     inlet["u"] = inlet[velvar][::, 0]
     inlet["v"] = inlet[velvar][::, 1]
@@ -33,6 +31,11 @@ def blade_loading_cp(
     v = massflowave_plane(inlet, valname="v", rhoname=densityvar, velocityname=velvar)
     w = massflowave_plane(inlet, valname="w", rhoname=densityvar, velocityname=velvar)
     U = vecAbs([u, v, w])
+
+    px_min = np.min(case_instance.blade.ps_pv.points[:, 0])
+    px_max = np.max(case_instance.blade.ps_pv.points[:, 0])
+    cax_len = px_max - px_min
+
     pt1 = p1 + 1 / 2 * rho * U ** 2
 
     ssmeshpoints = case_instance.blade.ss_pv
@@ -42,14 +45,16 @@ def blade_loading_cp(
     ps_cp = np.zeros(psmeshpoints.number_of_points)
 
     for idx, pts1 in enumerate(psmeshpoints.points):
-        ps_xc[idx] = pts1[0] / camber_length
+        #minmax norm ps_xc
+
+        ps_xc[idx] = (pts1[0] - px_min) / cax_len
         ps_cp[idx] = calc_inflow_cp(psmeshpoints.point_data[pressurevar][idx], pt1, p1)
 
     ss_xc = np.zeros(ssmeshpoints.number_of_points)
     ss_cp = np.zeros(ssmeshpoints.number_of_points)
 
     for idx, pts1 in enumerate(ssmeshpoints.points):
-        ss_xc[idx] = pts1[0] / camber_length
+        ss_xc[idx] = (pts1[0] - px_min) / cax_len
         ss_cp[idx] = calc_inflow_cp(ssmeshpoints.point_data[pressurevar][idx], pt1, p1)
 
     plt.figure()

ntrfc/turbo/cascade_case/solution.py

@@ -71,5 +71,6 @@ class GenericCascadeCase():
         bladepoly = pv.PolyData(self.active_blade_slice.points)
         for arr in self.active_blade_slice.point_data.keys():
             bladepoly[arr] = self.active_blade_slice.point_data[arr]
+
         self.blade = Blade2D(bladepoly)
         self.blade.compute_all_frompoints(alpha)

ntrfc/turbo/cascade_case/utils/domain_utils.py

@@ -97,9 +97,9 @@ class Blade2D:
         self.compute_skeleton()
         self.compute_stats()
 
-    def plot(self, figurepath=tempfile.mkdtemp() + "/plot.png", window_size=[2400, 2400]):
+    def plot(self, figurepath=tempfile.mkdtemp() + "/plot.png", window_size=[2400, 2400],off_screen=True):
         sfactor = window_size[0] / 200
-        p = pv.Plotter(off_screen=True)
+        p = pv.Plotter(off_screen=off_screen)
         p.add_mesh(self.ps_pv, color="red", label="pressure side", line_width=1)
         p.add_mesh(self.ss_pv, color="blue", label="suction side", line_width=1)
         p.add_mesh(self.skeletonline_pv, color="black", label="skeleton line", line_width=1)

ntrfc/turbo/profile_tele_extraction.py

@@ -86,13 +86,12 @@ def clean_sites(voronoi_sites, skeletonize_sites):
 def compute_midline_error(ind_1_start, ind_2_start, sortedPoly, shapelypoly):
     ps_pv, ss_pv = extractSidePolys(ind_1_start, ind_2_start, sortedPoly)
 
-    mid = lines_from_points(midline_from_sides(ps_pv, ss_pv, res=100).points)
+    mid = lines_from_points(midline_from_sides(ps_pv, ss_pv, res=200).points)
 
     mid_normals = mid.extrude([0, 0, 1]).compute_normals().slice("z").compute_cell_sizes()
 
 
     normals = mid_normals["Normals"]
-    lengths = mid_normals["Length"]
     centers = mid_normals.cell_centers().points
 
     raylength=np.sum(mid_normals["Length"])
@@ -113,8 +112,19 @@ def compute_midline_error(ind_1_start, ind_2_start, sortedPoly, shapelypoly):
     differences = np.abs(lengths_positive - lengths_negative) / ((lengths_positive + lengths_negative) / 2)
 
     rms_error = np.sqrt(np.mean(differences ** 2))
+    max_error = np.max(differences)
+    end_error = (sum(differences[-5:])+sum(differences[:5]))/10
 
-    return rms_error
+    p= pv.Plotter()
+    midpoly = pv.PolyData(centers+(ps_pv.points[0]-centers[0]))
+    midpoly["lengths"] = differences
+    p.add_mesh(midpoly,scalars="lengths")
+    p.add_mesh(sortedPoly)
+    p.add_mesh(ps_pv, color="red", label="pressure side", line_width=1)
+    p.add_mesh(ss_pv, color="blue", label="suction side", line_width=1)
+    p.add_text(f"RMS Error: {rms_error:.2e}\nMax Error: {max_error:.2e}\nEnd Error: {end_error:.2e}", font_size=10)
+    p.show()
+    return end_error
 
 
 def extract_vk_hk(sortedPoly: pv.PolyData) -> (int, int):
@@ -128,14 +138,13 @@ def extract_vk_hk(sortedPoly: pv.PolyData) -> (int, int):
 
     #create sheply polygon from refined pointset
     shapelypoly = Polygon(points_2d_closed_refined_voronoi)
-    raylength = 1000
 
     voronoi_sites = voronoi_skeleton_sites(points_2d_closed_refined_voronoi)
     skeletonize_sites = skeletonize_skeleton_sites(points_2d_closed_refined_skeletonize)
     valid_midline = clean_sites(voronoi_sites, skeletonize_sites)
     sort_indices = np.argsort(valid_midline[:, 0])
     valid_midline_sorted = valid_midline[sort_indices]
-    smoothing = 0.00001
+    smoothing = 0.000001
 
     u_new = np.arange(0, 1, 1 / 1024)
     (tck, u), fp, ier, msg = splprep((valid_midline_sorted[::, 0], valid_midline_sorted[::, 1]), u=None, per=0, k=3,
@@ -147,7 +156,6 @@ def extract_vk_hk(sortedPoly: pv.PolyData) -> (int, int):
     le_ind, te_ind, skeletonline_complete = skeletonline_completion(2, points_orig, points_2d_closed_refined_voronoi,
                                                                     np.stack([x_new, y_new]).T)
 
-    rms_error_ref = compute_midline_error(le_ind, te_ind,sortedPoly,shapelypoly)
 
     # Define the objective function
     def objective_function(indices, sortedPoly, shapelypoly):
@@ -155,95 +163,32 @@ def extract_vk_hk(sortedPoly: pv.PolyData) -> (int, int):
         rms_error_le_up = compute_midline_error(le_ind, te_ind, sortedPoly, shapelypoly)
         return rms_error_le_up
 
-    @dataclass
-    class SkeletonResult:
-        le_ind: int
-        te_ind: int
-        rms_error: float
-
-
-    # optimize te
-    search_te_pos = 3
-    te_pos_ind = te_ind
-    search_te_neg = 3
-    te_neg_ind = te_ind
-    search_le_pos = 3
-    le_pos_ind = le_ind
-    search_le_neg = 3
-    le_neg_ind = le_ind
-    search_le_neg_te_pos = 3
-    le_neg_te_pos_le_ind = le_ind
-    le_neg_te_pos_te_ind = te_ind
-    search_le_pos_te_neg = 3
-    le_pos_te_neg_le_ind = le_ind
-    le_pos_te_neg_te_ind = te_ind
-    search_le_neg_te_neg = 3
-    le_neg_te_neg_le_ind = le_ind
-    le_neg_te_neg_te_ind = te_ind
-    search_le_pos_te_pos = 3
-    le_pos_te_pos_le_ind = le_ind
-    le_pos_te_pos_te_ind = te_ind
-
-    skeletonresults = []
-    skeletonresults.append(SkeletonResult(le_ind, te_ind, rms_error_ref))
-
-
-    while search_te_pos > 0:
-        te_pos_ind = (te_pos_ind + 1) % npts
-        rms_error_te_pos = objective_function([le_ind, te_pos_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_ind, te_pos_ind, rms_error_te_pos))
-        if rms_error_te_pos > rms_error_ref:
-            search_te_pos -= 1
-    while search_te_neg > 0:
-        te_neg_ind = (te_neg_ind - 1) % npts
-        rms_error_te_neg = objective_function([le_ind, te_neg_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_ind, te_neg_ind, rms_error_te_neg))
-        if rms_error_te_neg > rms_error_ref:
-            search_te_neg -= 1
-    while search_le_pos > 0:
-        le_pos_ind = (le_pos_ind + 1) % npts
-        rms_error_le_pos = objective_function([le_pos_ind, te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_pos_ind, te_ind, rms_error_le_pos))
-        if rms_error_le_pos > rms_error_ref:
-            search_le_pos -= 1
-    while search_le_neg > 0:
-        le_neg_ind = (le_neg_ind - 1) % npts
-        rms_error_le_neg = objective_function([le_neg_ind, te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_neg_ind, te_ind, rms_error_le_neg))
-        if rms_error_le_neg > rms_error_ref:
-            search_le_neg -= 1
-    while search_le_neg_te_pos > 0:
-        le_neg_te_pos_le_ind = (le_neg_te_pos_le_ind - 1) % npts
-        le_neg_te_pos_te_ind = (le_neg_te_pos_te_ind + 1) % npts
-        rms_error_le_neg_te_pos = objective_function([le_neg_te_pos_le_ind, le_neg_te_pos_te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_neg_te_pos_le_ind, le_neg_te_pos_te_ind, rms_error_le_neg_te_pos))
-        if rms_error_le_neg_te_pos > rms_error_ref:
-            search_le_neg_te_pos -= 1
-    while search_le_pos_te_neg > 0:
-        le_pos_te_neg_le_ind = (le_pos_te_neg_le_ind + 1) % npts
-        le_pos_te_neg_te_ind = (le_pos_te_neg_te_ind - 1) % npts
-        rms_error_le_pos_te_neg = objective_function([le_pos_te_neg_le_ind, le_pos_te_neg_te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_pos_te_neg_le_ind, le_pos_te_neg_te_ind, rms_error_le_pos_te_neg))
-        if rms_error_le_pos_te_neg > rms_error_ref:
-            search_le_pos_te_neg -= 1
-    while search_le_neg_te_neg > 0:
-        le_neg_te_neg_le_ind = (le_neg_te_neg_le_ind - 1) % npts
-        le_neg_te_neg_te_ind = (le_neg_te_neg_te_ind - 1) % npts
-        rms_error_le_neg_te_neg = objective_function([le_neg_te_neg_le_ind, le_neg_te_neg_te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_neg_te_neg_le_ind, le_neg_te_neg_te_ind, rms_error_le_neg_te_neg))
-        if rms_error_le_neg_te_neg > rms_error_ref:
-            search_le_neg_te_neg -= 1
-    while search_le_pos_te_pos > 0:
-        le_pos_te_pos_le_ind = (le_pos_te_pos_le_ind + 1) % npts
-        le_pos_te_pos_te_ind = (le_pos_te_pos_te_ind + 1) % npts
-        rms_error_le_pos_te_pos = objective_function([le_pos_te_pos_le_ind, le_pos_te_pos_te_ind], sortedPoly, shapelypoly)
-        skeletonresults.append(SkeletonResult(le_pos_te_pos_le_ind, le_pos_te_pos_te_ind, rms_error_le_pos_te_pos))
-        if rms_error_le_pos_te_pos > rms_error_ref:
-            search_le_pos_te_pos -= 1
-
-    best_combination = min(skeletonresults, key=lambda x: x.rms_error)
-
-    le_ind_best, te_ind_best = best_combination.le_ind, best_combination.te_ind
+    le_search_interval = 12
+    te_search_interval = 12
+
+    le_min = (le_ind - le_search_interval) % npts
+    te_min = (te_ind - te_search_interval) % npts
+
+    #    le_max = (le_ind + 20) % npts
+    #    te_max = (te_ind + 20) % npts
+
+    allowed_le = [(i + le_min) % npts for i in range(2 * le_search_interval)]
+    allowed_te = [(i + te_min) % npts for i in range(2 * te_search_interval)]
+
+    specific_combinations = np.array(np.meshgrid(allowed_le, allowed_te)).T.reshape(-1, 2)
+
+    # Initialize the best values
+    best_combination = None
+    best_rms_error = float('inf')
+
+    # Evaluate the objective function for each combination
+    for indices in specific_combinations:
+        rms_error = objective_function(indices, sortedPoly, shapelypoly)
+        if rms_error < best_rms_error:
+            best_rms_error = rms_error
+            best_combination = indices
+
+    le_ind_best, te_ind_best = best_combination[0], best_combination[1]
     return le_ind_best, te_ind_best
 
 
@@ -290,6 +235,7 @@ def voronoi_skeleton_sites(points_2d_closed_refined_voronoi):
 
 
 def skeletonize_skeleton_sites(points_2d_closed_refined_skeletonize):
+
     res = len(points_2d_closed_refined_skeletonize)
     dx = np.min(points_2d_closed_refined_skeletonize[:, 0])
     dy = np.min(points_2d_closed_refined_skeletonize[:, 1])

requirements.txt

 pip>=23
+numpy==1.26.4
 matplotlib==3.9.0
 pytest==8.2.0
 trame==3.6.0

tests/geometry/test_ntrfc_geometry.py

@@ -21,7 +21,6 @@ def test_extract_vk_hk(verbose=False):
     xs, ys = naca(naca_code, res, half_cosine_spacing=False, finite_te=False)
     sorted_poly = pv.PolyData(np.stack([xs, ys, np.zeros(len(xs))]).T)
     sorted_poly.rotate_z(angle, inplace=True)
-    X, Y = sorted_poly.points[::, 0], sorted_poly.points[::, 1]
     ind_1 = res
     ind_2 = 0
 
@@ -33,6 +32,7 @@ def test_extract_vk_hk(verbose=False):
         p.add_mesh(sorted_poly.points[ind_hk], color="yellow", point_size=20)
         p.add_mesh(sorted_poly.points[ind_vk], color="red", point_size=20)
         p.add_mesh(sorted_poly)
+        p.view_xy()
         p.show()
 
     assert all(sorted_poly.points[ind_2] == sorted_poly.points[ind_hk])

tests/turbo/cascadecase/domain/test_ntrfc_blade2d.py

@@ -102,7 +102,7 @@ def test_t106():
 
     blade = Blade2D(points)
     blade.compute_all_frompoints()
-    blade.plot()
+
     bladepts = blade.sortedpointsrolled_pv.points
 
     xs, ys = bladepts[:, 0], bladepts[:, 1]