Merge branch 'new_concavehull' into 'master'

new concaveHull function

alpha value is now deprecated

See merge request !18

examples/profilepressure_igv.py

+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()

ntrfc/geometry/alphashape.py

 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]