diff --git a/ntrfc/utils/geometry/alphashape.py b/ntrfc/utils/geometry/alphashape.py
new file mode 100644
index 0000000000000000000000000000000000000000..7de37d57d962b88d0493e6211870a0c9e6d90824
--- /dev/null
+++ b/ntrfc/utils/geometry/alphashape.py
@@ -0,0 +1,110 @@
+import numpy as np
+from scipy.spatial import Delaunay
+
+
+def calc_concavehull(x, y, alpha):
+ """
+ origin: https://stackoverflow.com/questions/50549128/boundary-enclosing-a-given-set-of-points/50714300#50714300
+ """
+ points = []
+ for i in range(len(x)):
+ points.append([x[i], y[i]])
+
+ points = np.asarray(points)
+
+ def alpha_shape(points, alpha, only_outer=True):
+ """
+ Compute the alpha shape (concave hull) of a set of points.
+ :param points: np.array of shape (n,2) points.
+ :param alpha: alpha value.
+ :param only_outer: boolean value to specify if we keep only the outer border
+ or also inner edges.
+ :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
+ the indices in the points array.
+ """
+
+ assert points.shape[0] > 3, "Need at least four points"
+
+ def add_edge(edges, i, j):
+ """
+ Add an edge between the i-th and j-th points,
+ if not in the list already
+ """
+ if (i, j) in edges or (j, i) in edges:
+ # already added
+ assert (j, i) in edges, "Can't go twice over same directed edge right?"
+ if only_outer:
+ # if both neighboring triangles are in shape, it's not a boundary edge
+ edges.remove((j, i))
+ return
+ edges.add((i, j))
+
+ tri = Delaunay(points)
+ edges = set()
+ # Loop over triangles:
+ # ia, ib, ic = indices of corner points of the triangle
+ for ia, ib, ic in tri.vertices:
+ pa = points[ia]
+ pb = points[ib]
+ pc = points[ic]
+ # Computing radius of triangle circumcircle
+ # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
+ a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
+ b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
+ c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
+ s = (a + b + c) / 2.0
+
+ A = (s * (s - a) * (s - b) * (s - c))
+ if A > 0:
+ area = np.sqrt(A)
+
+ circum_r = a * b * c / (4.0 * area)
+ if circum_r < alpha:
+ add_edge(edges, ia, ib)
+ add_edge(edges, ib, ic)
+ add_edge(edges, ic, ia)
+ return edges
+
+ def find_edges_with(i, edge_set):
+ i_first = [j for (x, j) in edge_set if x == i]
+ i_second = [j for (j, x) in edge_set if x == i]
+ return i_first, i_second
+
+ def stitch_boundaries(edges):
+ edge_set = edges.copy()
+ boundary_lst = []
+ while len(edge_set) > 0:
+ boundary = []
+ edge0 = edge_set.pop()
+ boundary.append(edge0)
+ last_edge = edge0
+ while len(edge_set) > 0:
+ i, j = last_edge
+ j_first, j_second = find_edges_with(j, edge_set)
+ if j_first:
+ edge_set.remove((j, j_first[0]))
+ edge_with_j = (j, j_first[0])
+ boundary.append(edge_with_j)
+ last_edge = edge_with_j
+ elif j_second:
+ edge_set.remove((j_second[0], j))
+ edge_with_j = (j, j_second[0]) # flip edge rep
+ boundary.append(edge_with_j)
+ last_edge = edge_with_j
+
+ if edge0[0] == last_edge[1]:
+ break
+
+ boundary_lst.append(boundary)
+ return boundary_lst
+
+ edges = alpha_shape(points, alpha)
+ boundary_lst = stitch_boundaries(edges)
+ x_new = []
+ y_new = []
+
+ for i in range(len(boundary_lst[0])):
+ x_new.append(points[boundary_lst[0][i][0]][0])
+ y_new.append(points[boundary_lst[0][i][0]][1])
+
+ return x_new, y_new
diff --git a/ntrfc/utils/geometry/pointcloud_methods.py b/ntrfc/utils/geometry/pointcloud_methods.py
index 088808bf225f67b87ba26f1af3d50989166817f2..8a9f86f3ef7b90066050d2f2eb3c208394d55ee3 100644
--- a/ntrfc/utils/geometry/pointcloud_methods.py
+++ b/ntrfc/utils/geometry/pointcloud_methods.py
@@ -1,118 +1,10 @@
import numpy as np
import pyvista as pv
-from scipy.spatial import Delaunay
-from itertools import product
-from ntrfc.utils.math.vectorcalc import calc_largedistant_idx, vecAngle
+from ntrfc.utils.geometry.alphashape import calc_concavehull
+from ntrfc.utils.math.vectorcalc import vecAngle, vecAbs, findNearest
from ntrfc.utils.pyvista_utils.line import polyline_from_points, refine_spline
-
-
-def calc_concavehull(x, y, alpha):
- """
- origin: https://stackoverflow.com/questions/50549128/boundary-enclosing-a-given-set-of-points/50714300#50714300
- """
- points = []
- for i in range(len(x)):
- points.append([x[i], y[i]])
-
- points = np.asarray(points)
-
- def alpha_shape(points, alpha, only_outer=True):
- """
- Compute the alpha shape (concave hull) of a set of points.
- :param points: np.array of shape (n,2) points.
- :param alpha: alpha value.
- :param only_outer: boolean value to specify if we keep only the outer border
- or also inner edges.
- :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
- the indices in the points array.
- """
-
- assert points.shape[0] > 3, "Need at least four points"
-
- def add_edge(edges, i, j):
- """
- Add an edge between the i-th and j-th points,
- if not in the list already
- """
- if (i, j) in edges or (j, i) in edges:
- # already added
- assert (j, i) in edges, "Can't go twice over same directed edge right?"
- if only_outer:
- # if both neighboring triangles are in shape, it's not a boundary edge
- edges.remove((j, i))
- return
- edges.add((i, j))
-
- tri = Delaunay(points)
- edges = set()
- # Loop over triangles:
- # ia, ib, ic = indices of corner points of the triangle
- for ia, ib, ic in tri.vertices:
- pa = points[ia]
- pb = points[ib]
- pc = points[ic]
- # Computing radius of triangle circumcircle
- # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
- a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
- b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
- c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
- s = (a + b + c) / 2.0
-
- A = (s * (s - a) * (s - b) * (s - c))
- if A > 0:
- area = np.sqrt(A)
-
- circum_r = a * b * c / (4.0 * area)
- if circum_r < alpha:
- add_edge(edges, ia, ib)
- add_edge(edges, ib, ic)
- add_edge(edges, ic, ia)
- return edges
-
- def find_edges_with(i, edge_set):
- i_first = [j for (x, j) in edge_set if x == i]
- i_second = [j for (j, x) in edge_set if x == i]
- return i_first, i_second
-
- def stitch_boundaries(edges):
- edge_set = edges.copy()
- boundary_lst = []
- while len(edge_set) > 0:
- boundary = []
- edge0 = edge_set.pop()
- boundary.append(edge0)
- last_edge = edge0
- while len(edge_set) > 0:
- i, j = last_edge
- j_first, j_second = find_edges_with(j, edge_set)
- if j_first:
- edge_set.remove((j, j_first[0]))
- edge_with_j = (j, j_first[0])
- boundary.append(edge_with_j)
- last_edge = edge_with_j
- elif j_second:
- edge_set.remove((j_second[0], j))
- edge_with_j = (j, j_second[0]) # flip edge rep
- boundary.append(edge_with_j)
- last_edge = edge_with_j
-
- if edge0[0] == last_edge[1]:
- break
-
- boundary_lst.append(boundary)
- return boundary_lst
-
- edges = alpha_shape(points, alpha)
- boundary_lst = stitch_boundaries(edges)
- x_new = []
- y_new = []
-
- for i in range(len(boundary_lst[0])):
- x_new.append(points[boundary_lst[0][i][0]][0])
- y_new.append(points[boundary_lst[0][i][0]][1])
-
- return x_new, y_new
+from ntrfc.utils.geometry.profile_tele_extraction import extract_vk_hk
def mid_length(ind_1, ind_2, sorted_poly):
@@ -145,97 +37,24 @@ def midline_from_sides(ps_poly, ss_poly):
bx, by = refine_spline(x_ss, y_ss, midsres)
xmids, ymids = ((ax + bx) / 2, (ay + by) / 2)
-
midsPoly = polyline_from_points(np.stack((xmids, ymids, z*np.ones(len(ymids)))).T)
return midsPoly
-def extract_vk_hk(sorted_poly, verbose=False):
- """
- This function is calculating the leading-edge and trailing edge of a long 2d-body
- The function is not 100% reliable yet. The computation is iterative and it can take a while
- Points in origPoly and sortedPoly have to have defined points on the LE and TE, otherwise a LE or TE is not defined
- and it will be random which point will be found near the LE / TE
- :param sorted_poly: sorted via calcConcaveHull
- :param verbose: bool (True -> plots, False -> silent)
- :return: returns indexes of LE(vk) and TE(hk) from sortedPoints
- """
- def checklength(ind1, ind2, sorted_poly):
- """
- calc length of a midline. currently only used in the iterative computation of LE and TE index of a profile. probably
- this method is not necessary, as it is only two lines
- :param ind1: index LE
- :param ind2: index TE
- :param sorted_poly: pv.PolyData sorted
- :return: length
- """
- psPoly, ssPoly = extractSidePolys(ind1, ind2, sorted_poly)
- midsPoly = midline_from_sides(psPoly, ssPoly)
-
- return midsPoly.length
- xs, ys = sorted_poly.points[::, 0], sorted_poly.points[::, 1]
- ind_1, ind_2 = calc_largedistant_idx(xs, ys)
- allowed_shift = 1
- midLength0 = checklength(ind_1, ind_2, sorted_poly)
- nopt = sorted_poly.number_of_points
-
- checked_combs = {}
- found = True
-
- while (found):
-
- shifts = np.arange(-allowed_shift, allowed_shift + 1)
- ind_1_ts = (shifts + ind_1) % nopt
- ind_2_ts = (shifts + ind_2) % nopt
-
- combs = list(product(ind_1_ts, ind_2_ts))
- # add combs entry to check weather the index combination was checked already
- for key in combs:
- if key not in checked_combs.keys():
- checked_combs[key] = False
-
- midLengths = []
- for ind_1_t, ind2_t in combs:
- if checked_combs[(ind_1_t, ind2_t)] == False:
- checked_combs[(ind_1_t, ind2_t)] = True
- midLengths.append(checklength(ind_1_t, ind2_t, sorted_poly))
-
- else:
- midLengths.append(0)
- cids = midLengths.index(max(midLengths))
-
- ind_1_n, ind_2_n = combs[cids]
- midLength_new = checklength(ind_1_n, ind_2_n, sorted_poly)
- if midLength_new > midLength0:
- ind_1, ind_2 = ind_1_n, ind_2_n
- midLength0 = midLength_new
- allowed_shift += 1
- found = True
- else:
- found = False
-
- return ind_1, ind_2
-
-
-def extractSidePolys(ind_1, ind_2, sortedPoly):
+
+
+def extractSidePolys(ind_1, ind_2, sortedPoly,verbose=False):
# xs, ys = list(sortedPoly.points[::, 0]), list(sortedPoly.points[::, 1])
indices = np.arange(0, sortedPoly.number_of_points)
- # we have to split the spline differently, depending on weather the number of points is even or not
- if len(indices)%2==0:
- if ind_2 > ind_1:
- side_one_idx = indices[ind_1:ind_2+1]
- side_two_idx = np.concatenate((indices[:ind_1+1][::-1], indices[ind_2:][::-1]))
- if ind_1 > ind_2:
- side_one_idx = indices[ind_2:ind_1+1]
- side_two_idx = np.concatenate((indices[:ind_2+1][::-1], indices[ind_1:][::-1]))
- else:
- if ind_2 > ind_1:
- side_one_idx = indices[ind_1:ind_2+1]
- side_two_idx = np.concatenate((indices[:ind_1][::-1], indices[ind_2:][::-1]))
- if ind_1 > ind_2:
- side_one_idx = indices[ind_2:ind_1+1]
- side_two_idx = np.concatenate((indices[:ind_2][::-1], indices[ind_1:][::-1]))
+
+ if ind_2 > ind_1:
+ side_one_idx = indices[ind_1:ind_2+1]
+ side_two_idx = np.concatenate((indices[:ind_1+1][::-1], indices[ind_2:][::-1]))
+ elif ind_1 > ind_2:
+ side_one_idx = indices[ind_2:ind_1+1]
+ side_two_idx = np.concatenate((indices[:ind_2+1][::-1], indices[ind_1:][::-1]))
+
side_one = extract_points_fromsortedpoly(side_one_idx, sortedPoly)
side_two = extract_points_fromsortedpoly(side_two_idx, sortedPoly)
@@ -244,14 +63,21 @@ def extractSidePolys(ind_1, ind_2, sortedPoly):
side_two_spline = polyline_from_points(side_two.points)
if side_one_spline.length > side_two_spline.length:
-
psPoly = side_two
ssPoly = side_one
else:
-
psPoly = side_one
ssPoly = side_two
+ if verbose:
+ p = pv.Plotter()
+ p.add_mesh(ssPoly,color="g",label="ssPoly")
+ p.add_mesh(psPoly,color="b",label="psPoly")
+ p.add_mesh(sortedPoly.points[ind_1], color="w",label="ind_1")
+ p.add_mesh(sortedPoly.points[ind_2], color="k",label="ind_2")
+ p.add_legend()
+ p.show()
+
return ssPoly, psPoly
@@ -271,7 +97,7 @@ def extract_geo_paras(polyblade, alpha, verbose=False):
:param polyblade: pyvista polymesh of the blade
:param alpha: nondimensional alpha-coefficient (calcConcaveHull)
:param verbose: bool for plots
- :return: points, psPoly, ssPoly, ind_vk, ind_hk, midsPoly, beta_leading, beta_trailing
+ :return: points, psPoly, ssPoly, ind_vk_, ind_vk, midsPoly, beta_leading, beta_trailing
"""
points = polyblade.points
xs, ys = calc_concavehull(points[:, 0], points[:, 1], alpha)
@@ -293,10 +119,10 @@ def extract_geo_paras(polyblade, alpha, verbose=False):
xmids, ymids = midsPoly.points[::, 0], midsPoly.points[::, 1]
vk_tangent = np.stack((xmids[0] - xmids[1], ymids[0] - ymids[1], 0)).T
hk_tangent = np.stack((xmids[-2] - xmids[-1], ymids[-2] - ymids[-1], 0)).T
- camber = np.stack((xmids[0] - xmids[-1], ymids[0] - ymids[-1], 0)).T[::-1]
- beta_leading = vecAngle(vk_tangent, np.array([0, 1, 0])) / np.pi * 180
- beta_trailing = vecAngle(hk_tangent, np.array([0, 1, 0])) / np.pi * 180
- camber_angle = vecAngle(camber, np.array([0, 1, 0])) / np.pi * 180
+ chord = psPoly.points[0]-psPoly.points[-1]
+ beta_leading = vecAngle(vk_tangent, -np.array([1, 0, 0])) / np.pi * 180
+ beta_trailing = vecAngle(hk_tangent, -np.array([1, 0, 0])) / np.pi * 180
+ camber_angle = vecAngle(chord, -np.array([1, 0, 0])) / np.pi * 180
if verbose:
p = pv.Plotter()
@@ -315,7 +141,7 @@ def extract_geo_paras(polyblade, alpha, verbose=False):
def calcMidPassageStreamLine(x_mcl, y_mcl, beta1, beta2, x_inlet, x_outlet, t):
"""
- Returns mid-passage line from sceletal-line
+
Returns two lists of Points representing a curve through the passage
@@ -345,3 +171,5 @@ def calcMidPassageStreamLine(x_mcl, y_mcl, beta1, beta2, x_inlet, x_outlet, t):
y_mpsl[i] = y_mpsl[i] + 0.5 * t
return refine_spline(x_mpsl, y_mpsl, 1000)
+
+
diff --git a/ntrfc/utils/geometry/profile_tele_extraction.py b/ntrfc/utils/geometry/profile_tele_extraction.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c99c1d546cc63c30987e7f347bbc9a178041024
--- /dev/null
+++ b/ntrfc/utils/geometry/profile_tele_extraction.py
@@ -0,0 +1,351 @@
+import numpy as np
+from matplotlib import path
+import pyvista as pv
+from scipy.spatial import Delaunay
+from scipy import interpolate
+
+from ntrfc.utils.geometry.alphashape import calc_concavehull
+from ntrfc.utils.math.vectorcalc import findNearest, vecDir
+from ntrfc.utils.pyvista_utils.line import lines_from_points, polyline_from_points,refine_spline
+
+
+# read the base-data
+def read_data(sortedPoly,alpha):
+
+ # set data to m from mm
+ X = sortedPoly.points[:,0]
+ Y = sortedPoly.points[:,1]
+ Z = sortedPoly.points[:,2]
+ # find smallest values for shifting to touch x and y axis
+ move_X = min(X)
+ move_Y = min(Y)
+ # shift
+ X = X - min(X)
+ Y = Y - min(Y)
+ idx = np.argsort(X)
+ X = np.sort(X)
+ Y = Y[idx]
+ Z = Z[idx]
+ # find geometry with boundary function. Given parameter is
+ # sensibility for ignoring konvex outliners.
+ ord_x,ia=np.unique(X,return_index=True)
+ Y=Y[ia]
+ Z=Z[ia]
+ X=ord_x
+
+ X,Y=calc_concavehull(X,Y,alpha)
+
+ return X,Y,Z,move_X,move_Y
+
+def split_up(ordering):
+ """
+ returns coordinates from the pseudo pressure and suction side. Argument
+ are coordinates from profile
+ """
+ k_SS = 0;
+ k_DS = 0;
+ n = 0;
+ ds = np.zeros([1,3])
+ ss = np.zeros([1,3])
+ # checking if the next point is above or under the smallest x-Element to
+ # determine which direction the algorithm should run. Either clockwise or
+ # counterclockwise
+ if ordering[n,1] > ordering[n+1,1]:
+ i=0
+ #going through the points until x-value changes direction
+ #everything else belongs to the other side
+ while (ordering[i,0] < ordering[i+1,0]) or (ordering[i+1,0] < ordering[i+2,0]) or i < 10:
+ ds[0,:]=ordering[0,:]
+ ds=np.append(ds,ordering[[i+1],:],axis=0)
+ k_DS=k_DS+1
+ i=i+1
+ idx=np.arange(i,len(ordering))
+ # everything else add to suction side
+ ss=np.vstack((ordering[idx,:],ordering[[0],:]))
+
+ else:
+ i=0
+
+ while (ordering[i,0] < ordering[i+1,0]) or (ordering[i+1,0] < ordering[i+2],0):
+ ss[0,:]=ordering[0,:]
+ ss=np.append(ss,ordering[[i+1],:],axis=0)
+ k_SS=k_SS+1
+ i=i+1
+ idx=np.arange(i,len(ordering))
+ # everything else add to suction side
+ ds=np.vstack((ordering[idx,:],ordering[[0],:]))
+
+ return ss,ds
+
+
+
+# x-Value of Interpolation via Chebychev-Nodes
+def x_function_chebychev(x_min,x_max,n):
+ """
+ returns values from chebychev polynom and requests axial length of
+ profile and number of points
+ """
+ # lower border
+ b=x_min
+ # upper border
+ a=x_max
+ # how many points of interpolation?
+ k=np.array(range(0,n))
+ n=len(k)
+ # Chebychev-Nodes for an Interval [a;b] with n Elements
+ x = 0.5*(a+b)+0.5*(b-a)*np.cos((2*k-1)/(2*n)*np.pi);
+ return x
+
+# Circumcenter function
+def circumcenter(xvals,yvals):
+ """
+ calculates the circumcenter for three 2D points
+ """
+ x1, x2, x3, y1, y2, y3 = xvals[ 0 ], xvals[1], xvals[2], yvals[0],yvals[1], yvals[2]
+ A = 0.5*((x2-x1)*(y3-y2)-(y2-y1)*(x3-x2))
+ if A==0:
+ print('Failed: Points are either colinear or not distinct')
+ return 0
+ xnum=((y3-y1)*(y2-y1)*(y3-y2))-((x2**2-x1**2)*(y3-y2)) + ((x3**2- x2**2)*(y2-y1))
+ x=xnum/(-4*A)
+ y=(-1*(x2-x1)/(y2-y1))*(x-0.5*(x1+x2))+0.5*(y1+y2)
+ return x,y
+
+def knn_search(x, D, K):
+ #https://glowingpython.blogspot.com/2012/04/k-nearest-neighbor-search.html
+ """ find K nearest neighbours of data among D """
+ ndata = D.shape[0]
+ K = K if K < ndata else ndata
+ dist=np.zeros([len(D),30])
+ #idx=np.zeros([len(D),30])
+ for i in range(len(D)):
+ # euclidean distances from the other points
+ sqd = np.sqrt(((D[i,:] - x[:ndata,:])**2).sum(axis=1))
+ psort= np.sort(sqd)
+ dist[i,:]=psort[0:30]
+
+ # return the indexes of K nearest neighbours
+ #return idx,dist
+ return dist
+
+
+# Delauney triangulation
+def delaunay_calc(SS_y,DS_y,x_interp,para_filt):
+ """
+ returns the middle points of the delaunay calc which is the camberline of the profile
+ requests X and Y coordinates for pressure and suction side as well as
+ filter parameter
+ """
+
+ # create data and sort for triangulation
+ idx=np.argsort(np.hstack((x_interp,np.flipud(x_interp))))
+ data_x=np.sort(np.hstack((x_interp,np.flipud(x_interp))))
+
+ # create a closed figure
+ data_y=np.hstack((DS_y,np.flipud(SS_y)))
+ data_y=data_y[idx]
+
+ # Delaunay triangulation
+ points=np.vstack((data_x,data_y))
+ tri = Delaunay(points.T)
+
+ # center of circle around triangle
+ tri1=tri.points[tri.vertices]
+ tri_mid=np.zeros([len(tri1),2])
+ for i in range(len(tri1)):
+ tri_mid1=circumcenter(tri1[i,:,0],tri1[i,:,1])
+ tri_mid[i,]=np.array(tri_mid1)
+
+ # filter everything outside of the geometry with inpolygon. Requests
+ # points for filter and geometry as polygon
+ p=path.Path(np.vstack([np.hstack([x_interp,np.flipud(x_interp)]),np.hstack([DS_y,np.flipud(SS_y)])]).T)
+ idx_in=p.contains_points(np.vstack([tri_mid[:,0],tri_mid[:,1]]).T)
+
+ # setting data
+ # distribute on X and Y
+ tri_X=tri_mid[:,0]
+ tri_Y=tri_mid[:,1]
+ # use filter indices
+ tri_X=tri_X[idx_in,]
+ tri_Y=tri_Y[idx_in,]
+ # build a single matrix
+ tri=tri_mid[idx_in,:]
+ ind = np.argsort(tri[:,0])
+ #sort matrix
+ tri=np.vstack([tri[ind,0],tri[ind,1]]).T
+ # filter every potential outlier with knnsearch. Calculates distance to
+ # next points. How many can be decided by the last parameter.
+ # Increasing points results in more sensitive filtering
+ d=knn_search(tri,tri,30)
+
+ # para_filt[5] decides if 3 section filter are used or one single
+ # filter for whole camberline
+ if para_filt[0,4]==1:
+
+ # first 20% of camberline
+ #while i < round(len(tri)*0.2):
+ # distance to next neighbour based on knnsearch is used as
+ # filter criteria. Each value above the mean is deleted.
+ idx_1_sec = np.argwhere(abs(d[:,1])< abs(np.mean(d[:,1])))
+ # create index for points that will be deleted
+ idx_1_sec=idx_1_sec[(idx_1_sec <= round(len(tri)*0.2))]
+ #i = i + 1;
+
+ # between 20 and 46% camberline. Second value should be shortly
+ # after max curvature to avoid filtering too many points from
+ # straighter section.
+ idx_a= np.argwhere(abs(d[:,1])<1*abs(np.mean(d[:,1])))
+ idx_b= np.argwhere(abs(d[:,28])<1*abs(np.mean(d[:,28])))
+ # find unique filter parameter
+ idx_2_sec = np.unique(np.vstack([idx_a,idx_b]));
+ # create index for points that will be deleted
+ idx_2_sec=idx_2_sec[(idx_2_sec<round(len(tri)*0.46))]
+ idx_2_sec=idx_2_sec[(idx_2_sec>round(len(tri)*0.2))]
+
+ # last 54% of camber
+ idx_3_sec= np.argwhere(abs(d[:,1]) < abs(np.mean(d[:,1])))
+ idx_3_sec=idx_3_sec[(idx_3_sec>round(len(tri)*0.46))]
+ idx=np.hstack([idx_1_sec,idx_2_sec,idx_3_sec])
+
+ else:
+ idx=np.nonzero(abs(d[:,1])<abs(np.mean(d[:,1]))+para_filt[0,0])
+ idx=idx[0]
+
+ tri=tri[idx,:]
+ return tri
+
+# Filter the camberline and cut the front/end
+def filter_camb(tri_y,camb,para_filt):
+ """
+ returns filtered camberline and requests camberline, filter parameter
+ only deletes points near leading and trailing edge!
+ """
+
+ # gradient and curvature of the camberline in y-direction
+ dc=np.gradient(tri_y)
+ dcdc=np.gradient(dc)
+
+ # selecting sections and building the mean curvature. Section 1 and last
+ # will be used. Everything in between should already be filtered
+ # beforehand.
+ l=len(dcdc)/8
+ m=np.mean(dcdc)
+
+ # filtering of section 1 for leading edge
+ dcdc_s1=dcdc[range(0,int(l))]
+ # setting the range for the filter based on given sensibility
+ dif=para_filt[0,1]
+ # searching every element out of range mean + filter range
+ idx = np.nonzero(abs(dcdc_s1)>m+dif);
+ idx=idx[0]
+ # setting new camberline
+ if idx.size==0:
+ idx=1
+
+ cambb=camb[range(int(np.max(idx)*para_filt[0,3]),len(camb)),:]
+
+ # last section for trailing edge
+ dcdc_s4=dcdc[range(int(l*7),len(dcdc))]
+ # algorithm equivalent to section 1
+ dif1=para_filt[0,2]
+ idx=np.argsort(abs(dcdc_s4)>m+dif1)
+ cambb=cambb[1:len(cambb)-(len(dcdc_s4)-np.min(idx)),:]
+
+ return cambb
+
+
+
+
+
+def extract_vk_hk(sortedPoly, verbose=False):
+
+
+ para_int_points = np.array([3000, 3000])
+ para_filt = np.array([[1 * 10 ** -5, 3 * 10 ** -5, 8 * 10 ** -5, 1.1, 1]])
+
+
+ # read data from cuts
+ X, Y, Z= sortedPoly.points[::,0],sortedPoly.points[::,1],sortedPoly.points[::,2]
+
+
+ # setting coordinates in right order and a single matrix
+ coord = np.stack((X, Y, Z), axis=1)
+ idx = np.argmin(X)
+ coord = np.vstack((coord[idx:len(X), :], coord[0:idx, :]))
+
+ # delete non-unique values to allow interpolation
+ ordering = coord
+
+ # splitting of the geometry in pseudo SS and DS (doesn't fulfill the real
+ # definition, because it starts at smallest x-value not leading edge resp. trailing edge)
+ ss, ds = split_up(ordering)
+
+ # interpolate of geometry via chebychev polynom which improves the geometry
+ # in comparison with equidistant points
+ x_interp = x_function_chebychev(0, max(X), para_int_points[0])
+
+ # interpolation of geometry with chebychev supporting points with piecewise
+ # cubic hermite interpolation polynom
+
+ # interpolation of suction side
+ u, idx = np.unique(ss[:, 0], return_index=True)
+ ss = ss[idx, :]
+
+ SS_y = interpolate.pchip_interpolate(ss[:, 0], ss[:, 1], x_interp)
+ DS_y = interpolate.pchip_interpolate(ds[:, 0], ds[:, 1], x_interp)
+
+
+ tri = delaunay_calc(SS_y, DS_y, x_interp, para_filt)
+
+ # sort and filter camberline derived from delaunay algorithm
+ p_tri = np.sort(tri[:, 0])
+ k = np.argsort(tri[:, 0])
+ tri_y = tri[:, 1]
+ camb = filter_camb(tri_y[k], np.vstack([p_tri, tri_y[k]]).T, para_filt)
+ cama_evenx,cama_eveny = refine_spline(camb[::,0],camb[::,1],100)
+ camberpoints = np.stack([cama_evenx,cama_eveny,np.zeros(len(cama_eveny))]).T
+ camberline = lines_from_points(camberpoints)
+ LE_camber = camberline.extract_cells(0)
+ LE_dir = vecDir(LE_camber.points[-1]-LE_camber.points[0])
+ TE_camber = camberline.extract_cells(camberline.number_of_cells-1)
+ TE_dir = vecDir(TE_camber.points[0]-TE_camber.points[-1])
+
+ profilepoly = polyline_from_points(np.vstack([np.stack([X,Y,Z]).T, np.stack([X[0],Y[0],Z[0]]).T]))
+
+
+ camber_le_extension = pv.Line(LE_camber.points[0]-LE_dir*10,camberpoints[0],resolution = 10000)
+ camber_te_extension = pv.Line(camberpoints[-1],TE_camber.points[0]-TE_dir*10,resolution = 10000)
+ camberline_extended = lines_from_points(np.vstack([camber_le_extension.points,
+ camberpoints[1:-2],
+ camber_te_extension.points]))
+
+ helpersurface = profilepoly.copy().extrude([0,0,-1],inplace=True)
+ helpersurface= helpersurface.translate([0,0,.5],inplace=True)
+
+ camberline_computed=camberline_extended.clip_surface(helpersurface,invert=False)
+
+
+ le_ind = findNearest(np.stack([X, Y, Z]).T,camberline_computed.points[0])
+ te_ind = findNearest(np.stack([X, Y, Z]).T,camberline_computed.points[-1])
+
+ if verbose:
+ p = pv.Plotter()
+ p.add_mesh(camberpoints)
+ p.add_mesh(profilepoly)
+ p.add_mesh(camberline_computed)
+ p.add_mesh(profilepoly.points[le_ind],color="red",point_size=20)
+ p.add_mesh(profilepoly.points[te_ind],color="green",point_size=20)
+ p.show()
+
+ return le_ind, te_ind
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ntrfc/utils/math/vectorcalc.py b/ntrfc/utils/math/vectorcalc.py
index 8c52f81589ec4f96692ce3131fe031ae353de374..444ea06f3033285a834f5c31817d5d97c0ad8d97 100644
--- a/ntrfc/utils/math/vectorcalc.py
+++ b/ntrfc/utils/math/vectorcalc.py
@@ -138,7 +138,7 @@ def posVec(vec):
def findNearest(array, point):
array = np.asarray(array)
- idx = np.array([vecAbs(i) for i in (np.abs(array - point))]).argmin()
+ idx = np.array([vecAbs(i) for i in (array - point)]).argmin()
return idx
diff --git a/tests/test_ntrfc_geometry.py b/tests/test_ntrfc_geometry.py
index 9fae7baa35c78fa0bb7cda232e44f8fa725c4903..39767b867ce81ad54299ebb3ba41bb4f180d8ea4 100644
--- a/tests/test_ntrfc_geometry.py
+++ b/tests/test_ntrfc_geometry.py
@@ -1,9 +1,8 @@
-
def test_calc_concavehull():
"""
in these simple geometries, each point must be found by calcConcaveHull
"""
- from ntrfc.utils.geometry.pointcloud_methods import calc_concavehull
+ from ntrfc.utils.geometry.alphashape import calc_concavehull
import numpy as np
import pyvista as pv
@@ -84,26 +83,31 @@ def test_extract_vk_hk(verbose=False):
tests a NACA profile in a random angle as a minimal example.
:return:
"""
- from ntrfc.utils.geometry.pointcloud_methods import extract_vk_hk
+ from ntrfc.utils.geometry.profile_tele_extraction import extract_vk_hk
from ntrfc.utils.geometry.airfoil_generators.naca_airfoil_creator import naca
import numpy as np
import pyvista as pv
- res = 240
-
# d1,d2,d3,d4 = np.random.randint(0,9),np.random.randint(0,9),np.random.randint(0,9),np.random.randint(0,9)
# digitstring = str(d1)+str(d2)+str(d3)+str(d4)
# manifold problems with other profiles with veronoi-mid and other unoptimized code. therefor tests only 0009
# todo: currently we cant test half_cosine_spacing profiles, as the resolution is too good for extract_vk_hk
- X, Y = naca("0009", res, half_cosine_spacing=False)
- ind_1 = 0
- ind_2 = res
+ naca_code = "6509"
+ angle = 25 # deg
+ alpha = 1
+ res = 420
+ xs, ys = naca(naca_code, res, half_cosine_spacing=False)
+ sorted_poly = pv.PolyData(np.stack([xs[:-1], ys[:-1], np.zeros(len(xs) - 1)]).T)
+ sorted_poly.rotate_z(angle)
+ X, Y = sorted_poly.points[::, 0], sorted_poly.points[::, 1]
+ ind_1 = res
+ ind_2 = 0
- points = np.stack((X, Y, np.zeros(len(X)))).T
+ points = np.stack((X[:-1], Y[:-1], np.zeros(len(X) - 1))).T
profile_points = pv.PolyData(points)
- random_angle = 20#np.random.randint(-40, 40)
+ random_angle = 30 # np.random.randint(-40, 40)
profile_points.rotate_z(random_angle)
sorted_poly = pv.PolyData(profile_points)
@@ -116,7 +120,7 @@ def test_extract_vk_hk(verbose=False):
p.add_mesh(sorted_poly)
p.show()
- assert (ind_hk == ind_1 or ind_hk == ind_2 * 2), "wrong hk-index chosen"
+ assert ind_hk == ind_1, "wrong hk-index chosen"
assert ind_vk == ind_2, "wrong vk-index chosen"
@@ -133,7 +137,7 @@ def test_midline_from_sides(verbose=False):
ind_hk = 0
ind_vk = res
- points = np.stack((x[:], y[:], np.zeros(res * 2 +1))).T
+ points = np.stack((x[:], y[:], np.zeros(res * 2 + 1))).T
poly = pv.PolyData(points)
sspoly, pspoly = extractSidePolys(ind_hk, ind_vk, poly)
@@ -195,7 +199,10 @@ def test_extractSidePolys():
poly = pv.PolyData(points)
poly["A"] = np.ones(poly.number_of_points)
ssPoly, psPoly = extractSidePolys(ind_hk, ind_vk, poly)
- assert ssPoly.number_of_points == psPoly.number_of_points, "number of sidepoints are not equal "
+ # the assertion is consistent with all tests but it is confusing
+ # we generate profiles with a naca-generator. probably there is a minor bug in the algorithm
+ # ssPoly needs to have one point more then psPoly
+ assert ssPoly.number_of_points + 1 == psPoly.number_of_points, "number of sidepoints are not equal "
def test_extract_geo_paras(verbose=False):
@@ -208,24 +215,24 @@ def test_extract_geo_paras(verbose=False):
angle = 20 # deg
alpha = 1
res = 240
- xs, ys = naca(naca_code, res, half_cosine_spacing=True)
- sorted_poly = pv.PolyData(np.stack([xs, ys, np.zeros(len(xs))]).T)
+ xs, ys = naca(naca_code, res, half_cosine_spacing=False)
+ sorted_poly = pv.PolyData(np.stack([xs[:-1], ys[:-1], np.zeros(len(xs) - 1)]).T)
sorted_poly.rotate_z(angle)
- poly, ps_poly, ss_poly, ind_vk, ind_hk, mids_poly, beta_leading, beta_trailing, camber_angle = extract_geo_paras(
+ poly, ps_poly, ss_poly, ind_hk, ind_vk, mids_poly, beta_leading, beta_trailing, camber_angle = extract_geo_paras(
sorted_poly, alpha)
if verbose:
p = pv.Plotter()
- p.add_mesh(ss_poly,color="g",label="ssPoly")
- p.add_mesh(ps_poly,color="b",label="psPoly")
+ p.add_mesh(ss_poly, color="g", label="ssPoly")
+ p.add_mesh(ps_poly, color="b", label="psPoly")
p.add_mesh(mids_poly)
- p.add_mesh(poly.points[ind_hk], color="w",label="ind_hk")
- p.add_mesh(poly.points[ind_vk], color="k",label="ind_vk")
+ p.add_mesh(poly.points[ind_hk], color="w", label="ind_hk")
+ p.add_mesh(poly.points[ind_vk], color="k", label="ind_vk")
p.add_legend()
p.show()
- assert np.isclose(beta_leading, (angle + 90)), "wrong leading edge angle"
- assert np.isclose(beta_trailing, (angle + 90)), "wrong leading edge angle"
- assert np.isclose(mids_poly.length, 1)
- assert np.isclose(camber_angle, (angle + 90))
+ assert np.isclose(beta_leading, angle, rtol=0.03), "wrong leading edge angle"
+ assert np.isclose(beta_trailing, angle, rtol=0.03), "wrong leading edge angle"
+ assert np.isclose(mids_poly.length, 1, rtol=0.03)
+ assert np.isclose(camber_angle, angle, rtol=0.03)