Skip to content
Snippets Groups Projects
Commit a696c4c4 authored by Malte Nyhuis's avatar Malte Nyhuis
Browse files

New geo extraction

parent 88e4150b
No related branches found
No related tags found
2 merge requests!10Merge recent changes in master to we_ahrens,!6New geo extraction
Hide whitespace changes
Inline Side-by-side
ntrfc/utils/geometry/alphashape.py 0 → 100644
+ 110
− 0
View file @ a696c4c4
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
ntrfc/utils/geometry/pointcloud_methods.py
+ 31
− 203
View file @ a696c4c4
import numpy as np import numpy as np
import pyvista as pv 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 from ntrfc.utils.pyvista_utils.line import polyline_from_points, refine_spline
from ntrfc.utils.geometry.profile_tele_extraction import extract_vk_hk
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
def mid_length(ind_1, ind_2, sorted_poly): def mid_length(ind_1, ind_2, sorted_poly):
...@@ -145,97 +37,24 @@ def midline_from_sides(ps_poly, ss_poly): ...@@ -145,97 +37,24 @@ def midline_from_sides(ps_poly, ss_poly):
bx, by = refine_spline(x_ss, y_ss, midsres) bx, by = refine_spline(x_ss, y_ss, midsres)
xmids, ymids = ((ax + bx) / 2, (ay + by) / 2) xmids, ymids = ((ax + bx) / 2, (ay + by) / 2)
midsPoly = polyline_from_points(np.stack((xmids, ymids, z*np.ones(len(ymids)))).T) midsPoly = polyline_from_points(np.stack((xmids, ymids, z*np.ones(len(ymids)))).T)
return midsPoly 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 def extractSidePolys(ind_1, ind_2, sortedPoly,verbose=False):
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):
# xs, ys = list(sortedPoly.points[::, 0]), list(sortedPoly.points[::, 1]) # xs, ys = list(sortedPoly.points[::, 0]), list(sortedPoly.points[::, 1])
indices = np.arange(0, sortedPoly.number_of_points) 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:
if ind_2 > ind_1: side_one_idx = indices[ind_1:ind_2+1]
side_one_idx = indices[ind_1:ind_2+1] side_two_idx = np.concatenate((indices[:ind_1+1][::-1], indices[ind_2:][::-1]))
side_two_idx = np.concatenate((indices[:ind_1+1][::-1], indices[ind_2:][::-1])) elif ind_1 > ind_2:
if ind_1 > ind_2: side_one_idx = indices[ind_2:ind_1+1]
side_one_idx = indices[ind_2:ind_1+1] side_two_idx = np.concatenate((indices[:ind_2+1][::-1], indices[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]))
side_one = extract_points_fromsortedpoly(side_one_idx, sortedPoly) side_one = extract_points_fromsortedpoly(side_one_idx, sortedPoly)
side_two = extract_points_fromsortedpoly(side_two_idx, sortedPoly) side_two = extract_points_fromsortedpoly(side_two_idx, sortedPoly)
...@@ -244,14 +63,21 @@ def extractSidePolys(ind_1, ind_2, sortedPoly): ...@@ -244,14 +63,21 @@ def extractSidePolys(ind_1, ind_2, sortedPoly):
side_two_spline = polyline_from_points(side_two.points) side_two_spline = polyline_from_points(side_two.points)
if side_one_spline.length > side_two_spline.length: if side_one_spline.length > side_two_spline.length:
psPoly = side_two psPoly = side_two
ssPoly = side_one ssPoly = side_one
else: else:
psPoly = side_one psPoly = side_one
ssPoly = side_two 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 return ssPoly, psPoly
...@@ -271,7 +97,7 @@ def extract_geo_paras(polyblade, alpha, verbose=False): ...@@ -271,7 +97,7 @@ def extract_geo_paras(polyblade, alpha, verbose=False):
:param polyblade: pyvista polymesh of the blade :param polyblade: pyvista polymesh of the blade
:param alpha: nondimensional alpha-coefficient (calcConcaveHull) :param alpha: nondimensional alpha-coefficient (calcConcaveHull)
:param verbose: bool for plots :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 points = polyblade.points
xs, ys = calc_concavehull(points[:, 0], points[:, 1], alpha) xs, ys = calc_concavehull(points[:, 0], points[:, 1], alpha)
...@@ -293,10 +119,10 @@ def extract_geo_paras(polyblade, alpha, verbose=False): ...@@ -293,10 +119,10 @@ def extract_geo_paras(polyblade, alpha, verbose=False):
xmids, ymids = midsPoly.points[::, 0], midsPoly.points[::, 1] xmids, ymids = midsPoly.points[::, 0], midsPoly.points[::, 1]
vk_tangent = np.stack((xmids[0] - xmids[1], ymids[0] - ymids[1], 0)).T 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 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] chord = psPoly.points[0]-psPoly.points[-1]
beta_leading = vecAngle(vk_tangent, np.array([0, 1, 0])) / np.pi * 180 beta_leading = vecAngle(vk_tangent, -np.array([1, 0, 0])) / np.pi * 180
beta_trailing = vecAngle(hk_tangent, np.array([0, 1, 0])) / np.pi * 180 beta_trailing = vecAngle(hk_tangent, -np.array([1, 0, 0])) / np.pi * 180
camber_angle = vecAngle(camber, np.array([0, 1, 0])) / np.pi * 180 camber_angle = vecAngle(chord, -np.array([1, 0, 0])) / np.pi * 180
if verbose: if verbose:
p = pv.Plotter() p = pv.Plotter()
...@@ -315,7 +141,7 @@ def extract_geo_paras(polyblade, alpha, verbose=False): ...@@ -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): 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 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): ...@@ -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 y_mpsl[i] = y_mpsl[i] + 0.5 * t
return refine_spline(x_mpsl, y_mpsl, 1000) return refine_spline(x_mpsl, y_mpsl, 1000)
ntrfc/utils/geometry/profile_tele_extraction.py 0 → 100644
+ 351
− 0
View file @ a696c4c4
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
ntrfc/utils/math/vectorcalc.py
+ 1
− 1
View file @ a696c4c4
...@@ -138,7 +138,7 @@ def posVec(vec): ...@@ -138,7 +138,7 @@ def posVec(vec):
def findNearest(array, point): def findNearest(array, point):
array = np.asarray(array) 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 return idx
......
tests/test_ntrfc_geometry.py
+ 31
− 24
View file @ a696c4c4
def test_calc_concavehull(): def test_calc_concavehull():
""" """
in these simple geometries, each point must be found by calcConcaveHull 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 numpy as np
import pyvista as pv import pyvista as pv
...@@ -84,26 +83,31 @@ def test_extract_vk_hk(verbose=False): ...@@ -84,26 +83,31 @@ def test_extract_vk_hk(verbose=False):
tests a NACA profile in a random angle as a minimal example. tests a NACA profile in a random angle as a minimal example.
:return: :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 from ntrfc.utils.geometry.airfoil_generators.naca_airfoil_creator import naca
import numpy as np import numpy as np
import pyvista as pv 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) # 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) # 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 # 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 # 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) naca_code = "6509"
ind_1 = 0 angle = 25 # deg
ind_2 = res 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) 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) profile_points.rotate_z(random_angle)
sorted_poly = pv.PolyData(profile_points) sorted_poly = pv.PolyData(profile_points)
...@@ -116,7 +120,7 @@ def test_extract_vk_hk(verbose=False): ...@@ -116,7 +120,7 @@ def test_extract_vk_hk(verbose=False):
p.add_mesh(sorted_poly) p.add_mesh(sorted_poly)
p.show() 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" assert ind_vk == ind_2, "wrong vk-index chosen"
...@@ -133,7 +137,7 @@ def test_midline_from_sides(verbose=False): ...@@ -133,7 +137,7 @@ def test_midline_from_sides(verbose=False):
ind_hk = 0 ind_hk = 0
ind_vk = res 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) poly = pv.PolyData(points)
sspoly, pspoly = extractSidePolys(ind_hk, ind_vk, poly) sspoly, pspoly = extractSidePolys(ind_hk, ind_vk, poly)
...@@ -195,7 +199,10 @@ def test_extractSidePolys(): ...@@ -195,7 +199,10 @@ def test_extractSidePolys():
poly = pv.PolyData(points) poly = pv.PolyData(points)
poly["A"] = np.ones(poly.number_of_points) poly["A"] = np.ones(poly.number_of_points)
ssPoly, psPoly = extractSidePolys(ind_hk, ind_vk, poly) 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): def test_extract_geo_paras(verbose=False):
...@@ -208,24 +215,24 @@ def test_extract_geo_paras(verbose=False): ...@@ -208,24 +215,24 @@ def test_extract_geo_paras(verbose=False):
angle = 20 # deg angle = 20 # deg
alpha = 1 alpha = 1
res = 240 res = 240
xs, ys = naca(naca_code, res, half_cosine_spacing=True) xs, ys = naca(naca_code, res, half_cosine_spacing=False)
sorted_poly = pv.PolyData(np.stack([xs, ys, np.zeros(len(xs))]).T) sorted_poly = pv.PolyData(np.stack([xs[:-1], ys[:-1], np.zeros(len(xs) - 1)]).T)
sorted_poly.rotate_z(angle) 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) sorted_poly, alpha)
if verbose: if verbose:
p = pv.Plotter() p = pv.Plotter()
p.add_mesh(ss_poly,color="g",label="ssPoly") p.add_mesh(ss_poly, color="g", label="ssPoly")
p.add_mesh(ps_poly,color="b",label="psPoly") p.add_mesh(ps_poly, color="b", label="psPoly")
p.add_mesh(mids_poly) p.add_mesh(mids_poly)
p.add_mesh(poly.points[ind_hk], color="w",label="ind_hk") 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_vk], color="k", label="ind_vk")
p.add_legend() p.add_legend()
p.show() p.show()
assert np.isclose(beta_leading, (angle + 90)), "wrong leading edge angle" assert np.isclose(beta_leading, angle, rtol=0.03), "wrong leading edge angle"
assert np.isclose(beta_trailing, (angle + 90)), "wrong leading edge angle" assert np.isclose(beta_trailing, angle, rtol=0.03), "wrong leading edge angle"
assert np.isclose(mids_poly.length, 1) assert np.isclose(mids_poly.length, 1, rtol=0.03)
assert np.isclose(camber_angle, (angle + 90)) assert np.isclose(camber_angle, angle, rtol=0.03)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment