diff --git a/ntrfc/postprocessing/turbo/__init__.py b/ntrfc/postprocessing/turbo/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ntrfc/postprocessing/turbo/bladeloading.py b/ntrfc/postprocessing/turbo/bladeloading.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc18ccecadce72b12061699cd563d69877a8f9ef
--- /dev/null
+++ b/ntrfc/postprocessing/turbo/bladeloading.py
@@ -0,0 +1,73 @@
+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
diff --git a/ntrfc/utils/filehandling/mesh.py b/ntrfc/utils/filehandling/mesh.py
index 3d1bbfeae6b50651f35fa6dbb0c6b48554aa716d..2990442c71e54153780e474766c8954e0d068a05 100644
--- a/ntrfc/utils/filehandling/mesh.py
+++ b/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
diff --git a/ntrfc/utils/geometry/pointcloud_methods.py b/ntrfc/utils/geometry/pointcloud_methods.py
index bc1d6d93852d05a67faf49ccc779dd99c0c27226..feafbae2cb0a3b79b11643bb503ce70d619e90c2 100644
--- a/ntrfc/utils/geometry/pointcloud_methods.py
+++ b/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):
diff --git a/tests/test_ntrfc.py b/tests/test_ntrfc.py
index 283cb744c30a4d0b315e32a4f6e775789d52879b..99e198f10c9e6cc213aa8fa16c3dd478ba25f358 100644
--- a/tests/test_ntrfc.py
+++ b/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
diff --git a/tests/test_ntrfc_geometry.py b/tests/test_ntrfc_geometry.py
index 9065ca21df1b59e2cdec714b674089352fd72eba..23582abac60124a44f636d5d8959c0eb91f3cdcd 100644
--- a/tests/test_ntrfc_geometry.py
+++ b/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"