blade-loading implementation

ntrfc/postprocessing/turbo/bladeloading.py

+from ntrfc.utils.geometry.pointcloud_methods import  extract_profilepoints
+from ntrfc.utils.pyvista_utils.plane import areaave_plane
+import matplotlib.pyplot as plt
+import numpy as np
+import pyvista as pv
+
+def calc_loading_volmesh(volmesh, bladesurface,alpha, verbose=False):
+    #todo: test method is missing
+    bounds = volmesh.bounds
+
+    zm = (bounds[5]-bounds[4])/2
+    midslice_z = bladesurface.slice(origin=(0,0,zm),normal=(0,0,1))
+
+    sortedPoly, psVals, ssVals, ind_vk, ind_hk, midsPoly, beta_leading, beta_trailing, camber_angle= extract_profilepoints(midslice_z,alpha)
+    sortedPoints = sortedPoly.points
+
+    x1 = bounds[0] + 1e-5 * bounds[1]
+    x2 = bounds[1] + 1e-5 * bounds[0]
+
+    inlet = volmesh.slice(normal="x", origin=(x1, 0, 0)).compute_cell_sizes().point_data_to_cell_data()
+    outlet = volmesh.slice(normal="x", origin=(x2, 0, 0)).compute_cell_sizes().point_data_to_cell_data()
+    p1 = areaave_plane(inlet,"p")
+    p2 = areaave_plane(outlet,"p")
+
+    camber = pv.Line((0, 0, 0), -(sortedPoints[ind_vk] - sortedPoints[ind_hk]))
+
+    shift = sortedPoints[ind_vk]
+    shift -= psVals.points[0][-1]
+
+    ssVals.points -= shift
+    psVals.points -= shift
+
+    ssVals.rotate_z(-camber_angle)
+    psVals.rotate_z(-camber_angle)
+
+    psVals = psVals.cell_data_to_point_data()
+    ssVals = ssVals.cell_data_to_point_data()
+
+    ps_xc = np.zeros(psVals.number_of_points)
+    ps_cp = np.zeros(psVals.number_of_points)
+
+    for idx, pt in enumerate(psVals.points):
+        ps_xc[idx] = pt[0] / camber.length
+        ps_cp[idx] = calc_inflow_cp(psVals.point_data["p"][idx], p2, p1)
+
+    ss_xc = np.zeros(ssVals.number_of_points)
+    ss_cp = np.zeros(ssVals.number_of_points)
+
+    for idx, pt in enumerate(ssVals.points):
+        ss_xc[idx] = pt[0] / camber.length
+        ss_cp[idx] = calc_inflow_cp(ssVals.point_data["p"][idx], p2, p1)
+
+    ssVals["xc"] = ss_xc
+    ssVals["cp"] = ss_cp
+    psVals["xc"] = ps_xc
+    psVals["cp"] = ps_cp
+
+    if verbose:
+        plt.figure()
+        plt.plot(ss_xc, ss_cp)
+        plt.plot(ps_xc, ps_cp)
+        plt.show()
+    return psVals, ssVals
+
+def calc_inflow_cp(px, pt1, p1):
+    """
+    :param px: pressure at position
+    :param pt1: total pressure inlet
+    :param p1: pressure inlet
+    :return: lift coefficient
+    """
+    cp = (px - pt1) / (p1 - pt1)
+    return cp

ntrfc/utils/filehandling/mesh.py

@@ -6,9 +6,8 @@ import pyvista as pv
 def load_mesh(path_to_mesh):
     assert os.path.isfile(path_to_mesh), path_to_mesh + " is not a valid file"
     extension = os.path.splitext(path_to_mesh)[1]
-    if extension == ".vtk":
+    if extension == ".vtk" or extension == ".vtp" or extension == ".vtu":
         mesh = pv.read(path_to_mesh)
-
     elif extension == ".cgns":
         cgns = cgnsReader(path_to_mesh)
         if str(type(cgns)).find('vtkStructuredGrid') != -1:
@@ -16,28 +15,8 @@ def load_mesh(path_to_mesh):
         elif str(type(cgns)).find('vtkUnstructuredGrid') != -1:
             mesh = pv.UnstructuredGrid(cgns)
         elif str(type(cgns)).find('vtkMultiBlockDataSet') != -1:
-            #mesh = pv.UnstructuredGrid()
             multiBlockMesh = pv.MultiBlock(cgns)
-            """for domainId in range(multiBlockMesh.GetNumberOfBlocks()):
-                domain = multiBlockMesh.GetBlock(domainId)
-                for blockId in range(domain.GetNumberOfBlocks()):
-                    block = domain.GetBlock(blockId)
-                    mesh = mesh.merge(block)
-            """
             mesh = multiBlockMesh.combine()
-
-    array_names = mesh.array_names
-    for name in array_names:
-        if name == "Velocity":
-            mesh.rename_array("Velocity", "U")
-        if name == "Pressure":
-            mesh.rename_array("Pressure", "p")
-        if name == "Density":
-            mesh.rename_array("Density", "rho")
-        if name == "Temperature":
-            mesh.rename_array("Temperature", "T")
-    mesh = mesh.cell_data_to_point_data()
-
     return mesh
 
 

ntrfc/utils/geometry/pointcloud_methods.py

@@ -243,22 +243,30 @@ def extractSidePolys(ind_hk, ind_vk, sortedPoly):
         y_ps = ys[ind_vk:] + ys[:ind_hk + 1]
         x_ps = xs[ind_vk:] + xs[:ind_hk + 1]
 
-    psl_helper = polyline_from_points(np.stack((x_ps, y_ps, np.zeros(len(x_ps)))).T)
-    ssl_helper = polyline_from_points(np.stack((x_ss, y_ss, np.zeros(len(x_ss)))).T)
-
-    if psl_helper.length > ssl_helper.length:
-
-        psPoly = pv.PolyData(ssl_helper.points)
-        ssPoly = pv.PolyData(psl_helper.points)
-    else:
-
-        psPoly = pv.PolyData(psl_helper.points)
-        ssPoly = pv.PolyData(ssl_helper.points)
+    ss_ids = []
+    ps_ids = []
+
+    for ix,iy in zip(x_ss,y_ss):
+        if ix in xs and iy in ys:
+            xid = np.where(xs == ix)[0][0]
+            yid = np.where(ys == iy)[0][0]
+            if xid==yid:
+                ss_ids.append(xid)
+
+    for ix, iy in zip(x_ps, y_ps):
+        if ix in xs and iy in ys:
+            xid = np.where(xs == ix)[0][0]
+            yid = np.where(ys == iy)[0][0]
+            if xid == yid:
+                ps_ids.append(xid)
+
+    psPoly = sortedPoly.extract_points(ps_ids)
+    ssPoly = sortedPoly.extract_points(ss_ids)
 
     return ssPoly, psPoly
 
 
-def extract_geo_paras(points, alpha, verbose=False):
+def extract_profilepoints(poly, alpha, verbose=False):
     """
     This function is extracting profile-data as stagger-angle, midline, psPoly, ssPoly and more from a set of points
     Be careful, you need a suitable alpha-parameter in order to get the right geometry
@@ -269,10 +277,19 @@ def extract_geo_paras(points, alpha, verbose=False):
     :param verbose: bool for plots
     :return: points, psPoly, ssPoly, ind_vk, ind_hk, midsPoly, beta_leading, beta_trailing
     """
-
+    points = poly.points
+    oxs,oys = points[:,0],points[:,1]
     xs, ys = calc_concavehull(points[:, 0], points[:, 1], alpha)
+
+    bladeids = []
+    for ix,iy in zip(xs,ys):
+        if ix in oxs and iy in oys:
+            xid = np.where(oxs == ix)[0][0]
+            yid = np.where(oys == iy)[0][0]
+            if xid==yid:
+                bladeids.append(xid)
     points = np.stack((xs, ys, np.zeros(len(xs)))).T
-    sortedPoly = pv.PolyData(points)
+    sortedPoly = poly.extract_points(bladeids)
 
     ind_hk, ind_vk = extract_vk_hk(sortedPoly)
     psPoly, ssPoly = extractSidePolys(ind_hk, ind_vk, sortedPoly)
@@ -297,7 +314,7 @@ def extract_geo_paras(points, alpha, verbose=False):
         p.add_legend()
         p.show()
 
-    return points, psPoly, ssPoly, ind_vk, ind_hk, midsPoly, beta_leading, beta_trailing, camber_angle
+    return sortedPoly,psPoly, ssPoly, ind_vk, ind_hk, midsPoly, beta_leading, beta_trailing, camber_angle
 
 
 def calcMidPassageStreamLine(x_mcl, y_mcl, beta1, beta2, x_inlet, x_outlet, t):

tests/test_ntrfc.py

@@ -68,7 +68,7 @@ def test_loadmesh_vtk(tmpdir):
 
     test_file = tmpdir / "tmp.vtk"
     mesh = pv.Box()
-    mesh["Density"] = np.ones(mesh.number_of_cells)
+    mesh["rho"] = np.ones(mesh.number_of_cells)
     mesh.save(test_file)
     mesh_load = load_mesh(test_file)
     assert "rho" in mesh_load.array_names

tests/test_ntrfc_geometry.py

@@ -192,7 +192,7 @@ def test_extractSidePolys():
 
 
 def test_extract_geo_paras():
-    from ntrfc.utils.geometry.pointcloud_methods import extract_geo_paras
+    from ntrfc.utils.geometry.pointcloud_methods import extract_profilepoints
     from ntrfc.utils.geometry.airfoil_generators.naca_airfoil_creator import naca
     import numpy as np
     import pyvista as pv
@@ -205,8 +205,8 @@ def test_extract_geo_paras():
     sorted_poly = pv.PolyData(np.stack([xs, ys, np.zeros(len(xs))]).T)
     sorted_poly.rotate_z(angle)
 
-    points, ps_poly, ss_poly, ind_vk, ind_hk, mids_poly, beta_leading, beta_trailing, camber_angle = extract_geo_paras(
-        sorted_poly.points, alpha)
+    points, ps_poly, ss_poly, ind_vk, ind_hk, mids_poly, beta_leading, beta_trailing, camber_angle = extract_profilepoints(
+        sorted_poly, alpha)
 
     assert np.isclose(beta_leading, (angle + 90)), "wrong leading edge angle"
     assert np.isclose(beta_trailing, (angle + 90)), "wrong leading edge angle"