diff --git a/examples/profilepressure_igv.py b/examples/profilepressure_igv.py
new file mode 100644
index 0000000000000000000000000000000000000000..606b70dd78f930fe9c1b8ca7124f5f5b795dec06
--- /dev/null
+++ b/examples/profilepressure_igv.py
@@ -0,0 +1,30 @@
+from ntrfc.filehandling.mesh import load_mesh
+from ntrfc.geometry.plane import areaave_plane
+import matplotlib.pyplot as plt
+
+solution = load_mesh("igv_005_cgns.cgns")
+inlet =solution.slice(origin=(solution.bounds[0]+1e-12,0,0),normal=(1,0,0))
+inlet = inlet.ptc()
+pt1 = areaave_plane(inlet,"PressureStagnation")
+p1 = areaave_plane(inlet,"Pressure")
+
+outer = solution.extract_surface()
+
+outer_slice = outer.slice(origin=solution.center,normal=(0,0,1))
+regions = outer_slice.connectivity()
+blade = regions.extract_cells([i for i in range(regions.number_of_cells) if regions["RegionId"][i] == 1])
+
+camber_estimated = blade.bounds[1] - blade.bounds[0]
+
+xs = (blade.points[::,0]-blade.bounds[0])/camber_estimated
+cs = [(i-pt1)/(p1-pt1) for i in blade["Pressure"]]
+
+
+plt.figure(dpi=1200)
+plt.scatter(xs,cs,s=2)
+plt.xlabel("$c_{ax}$")
+plt.ylabel("$c_p$")
+plt.xlim([0, 1])
+plt.ylim([0.5, 2.5])
+plt.grid()
+plt.show()
diff --git a/ntrfc/geometry/alphashape.py b/ntrfc/geometry/alphashape.py
index 4c5df2adc856deeb7880cffd31f792b5085ff1fb..4e84bcf25170558e8fcf4dfaab9053a47e648f6f 100644
--- a/ntrfc/geometry/alphashape.py
+++ b/ntrfc/geometry/alphashape.py
@@ -1,110 +1,21 @@
import numpy as np
-from scipy.spatial import Delaunay
+def calc_concavehull(x, y, alpha=None):
-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]])
+ if alpha is not None:
+ import warnings
+ warnings.warn("The `alpha` parameter is deprecated and will be removed in a future version.", DeprecationWarning)
- points = np.asarray(points)
+ points = np.stack([x,y]).T
- 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.
- """
+ # Step 1: Find the center of the shape
+ center_x = points[:, 0].mean()
+ center_y = points[:, 1].mean()
- assert points.shape[0] > 3, "Need at least four points"
+ # Step 2: Calculate the angle of each point relative to the center
+ angles = np.arctan2(points[:, 1] - center_y, points[:, 0] - center_x)
- 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))
+ # Step 3: Sort the points based on their angles
+ angles, points = angles.argsort(), points[angles.argsort()]
- tri = Delaunay(points)
- edges = set()
- # Loop over triangles:
- # ia, ib, ic = indices of corner points of the triangle
- for ia, ib, ic in tri.simplices:
- 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
+ return points[::,0],points[::,1]