diff --git a/src/Meshing/Internal.py b/src/Meshing/Internal.py index 48a84a94..0e577533 100644 --- a/src/Meshing/Internal.py +++ b/src/Meshing/Internal.py @@ -171,8 +171,8 @@ class InternalMeshing(AMeshingTool): ltot1 = 0.0 ltot2 = 0.0 - alpha = [] - beta = [] + alpha = [0.0] + beta = [0.0] # remove names sect2.points[start2].name = '' sect2.points[end2].name = '' @@ -195,67 +195,38 @@ class InternalMeshing(AMeshingTool): else: # regular case alpha = list(map(lambda x: x/ltot1, alpha)) # target ratios beta = list(map(lambda x: x/ltot2, beta)) # current ratios - for i in range(len1 - len2): - beta.append(1.0) - beta.append(0.0) # beta[-1] - # loop on the points to insert - for k in range(len1-len2): - trouve = False - for j0 in range(len2): - if trouve: - break - undeplus = True - for j1 in range(j0, len2): - if undeplus: - if beta[j1] <= alpha[j1]: - undeplus = False - break - if undeplus: - if j0 == len2: - trouve = True - else: - j1 = j0+1 - if beta[j0]-alpha[j0] > alpha[j1]-beta[j0]: - trouve = True + match = self._match_indices(alpha, beta) - if trouve: - len2 += 1 - for i in reversed(range(j0, len2)): - beta[i+1] = beta[i] - if (beta[j0] == beta[j0-1]): - # rien - pass - else: - ratio = (alpha[j0] - beta[j0-1]) \ - / (beta[j0] - beta[j0-1]) - if ratio < 0.0: - logger.warning(f"ratio négatif {ratio}") - # on double le point a gauche - p = sect2.point(start2+j0-1).copy() - sect2.insert_point(start2+j0-1, p) - sect2.point(start2+j0-1).name = '' - beta[j0] = beta[j0-1] - else: - p = sect2.point(start2+j0).copy() - sect2.insert_point(start2+j0, p) - sect2.point(start2+j0+1).name = '' - sect2.point(start2+j0+1).x = \ - (1.0-ratio) * sect2.point(start2+j0).x + \ - ratio * sect2.point(start2+j0+2).x - sect2.point(start2+j0+1).y = \ - (1.0-ratio) * sect2.point(start2+j0).y + \ - ratio * sect2.point(start2+j0+2).y - sect2.point(start2+j0+1).z = \ - (1.0-ratio) * sect2.point(start2+j0).z + \ - ratio * sect2.point(start2+j0+2).z - beta[j0] = (1.-ratio)*beta[j0-1] \ - + ratio*beta[j0+1] + insertions_per_segment = [] # liste de (index, [points]) + for j in range(len2): + i_start = match[j] + i_end = match[j + 1] + q_start = sect2.point(start2 + j) + q_end = sect2.point(start2 + j + 1) + + new_points = [] + if i_end - i_start > 1: + u_start, u_end = alpha[i_start], alpha[i_end] + span = u_end - u_start + for k in range(i_start + 1, i_end): + ratio = (alpha[k] - u_start) / span # dans [0,1] + p = q_start.copy() + p.name = '' + p.x = (1.0-ratio) * q_start.x + ratio * q_end.x + p.y = (1.0-ratio) * q_start.y + ratio * q_end.y + p.z = (1.0-ratio) * q_start.z + ratio * q_end.z + new_points.append(p) + + # à insérer après l'indice j + insertions_per_segment.append((j, new_points)) + + # Insertion en ordre décroissant de j + # pour ne pas perturber les indices + for j, new_points in reversed(insertions_per_segment): + for offset, p in enumerate(new_points): + sect2.insert_point(start2 + j + 1 + offset, p) + len2 += 1 - if len1 > len2: - for i in range(len1 - len2): - p = sect2.point(start2+len2).copy() - sect2.insert_point(start2+len2, p) - sect2.point(start2+len2).name = '' if tag1 != "un": sect2.point(start2).name = tag1 if tag2 != "np": @@ -263,6 +234,71 @@ class InternalMeshing(AMeshingTool): sect2.modified() + def _match_indices(self, u_ref, v_target): + # appariement optimal (ordre croissant) par programmation dynamique + """ + Trouve, pour chaque v_target[j], un indice i_j dans u_ref tel que + i_0 < i_1 < ... < i_{m-1} et + Sum |u_ref[i_j] - v_target[j]| soit minimal. + + Retourne la liste des indices i_j (longueur m), + avec i_0 = 0 et i_{m-1} = n-1 imposés + (les extrémités se correspondent forcément). + + """ + n = len(u_ref) + m = len(v_target) + + # Les extrémités sont fixées. + fixed_first, fixed_last = 0, n - 1 + + # DP sur les indices disponibles entre fixed_first+1 et fixed_last-1 + # dp[j][i] = coût minimal pour apparier v_target[1..j] aux indices + # se terminant à l'indice i (i croissant strict), j de 1 à m-2. + if m == 2: + return [fixed_first, fixed_last] + + # indices internes possibles + candidates = list(range(fixed_first + 1, fixed_last)) + inner_targets = v_target[1:-1] # cibles internes (j = 1..m-2) + + INF = float("inf") + # dp[j] : dict indice_candidat -> (coût_min, provenance) + dp = [{} for _ in range(len(inner_targets))] + + for j, vt in enumerate(inner_targets): + for idx_pos, i in enumerate(candidates): + cost_here = abs(u_ref[i] - vt) + if j == 0: + # doit être > fixed_first + dp[j][i] = (cost_here, fixed_first) + else: + best = (INF, None) + for i_prev in candidates[:idx_pos]: + if i_prev in dp[j - 1]: + c = dp[j - 1][i_prev][0] + cost_here + if c < best[0]: + best = (c, i_prev) + if best[1] is not None: + dp[j][i] = best + + # Reconstruction : on part du dernier niveau, on choisit le + # meilleur indice tel que la suite jusqu'à fixed_last + # reste croissante et valide. + last_level = dp[-1] + valid = {i: c for i, c in last_level.items() if i < fixed_last} + if not valid: + raise RuntimeError("Alignement impossible avec ces données.") + best_i = min(valid, key=lambda i: valid[i][0]) + + chain = [best_i] + for j in range(len(inner_targets) - 1, 0, -1): + prev_i = dp[j][chain[-1]][1] + chain.append(prev_i) + chain.reverse() + + return [fixed_first] + chain + [fixed_last] + def update_rk(self, reach, begin_rk, end_rk, step: float = 50, limites=[-1, -1],