debug mesher

src/Meshing/Internal.py

@@ -57,6 +57,13 @@ class InternalMeshing(AMeshingTool):
 
         # we never modify original profiles, we work on the copies
         m_profiles = self.st_to_m(profiles, guide_list)
+
+        # Adapt the spacing to the kilometre-point direction. The value from
+        # the dialog is a distance, so only its magnitude is significant.
+        step = self.signed_step(
+            step, m_profiles[0].rk, m_profiles[-1].rk
+        )
+
         new_profiles = self.interpolate_transversal_step(m_profiles, step)
 
         for new_profiles2 in new_profiles:
@@ -65,6 +72,12 @@ class InternalMeshing(AMeshingTool):
 
         return new_profiles
 
+    @staticmethod
+    def signed_step(step, begin_rk, end_rk):
+        if end_rk < begin_rk:
+            return -abs(step)
+        return abs(step)
+
     def st_to_m(self, profiles, guide_list):
         guide_list = ["un"] + guide_list + ["np"]
         max_values = [0] * (len(guide_list) - 1)
@@ -158,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 = ''
@@ -182,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:
+            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
-                        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]
 
-        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":
@@ -250,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],

src/Model/Geometry/PointXYZ.py

@@ -95,10 +95,9 @@ class PointXYZ(Point):
                     "ADD COLUMN deleted BOOLEAN NOT NULL DEFAULT FALSE"
                 )
 
-        if major == "0" and minor == "2":
-            if int(release) < 5:
+        # Try to create index at each update
         execute(
-                    "CREATE INDEX idx_point " +
+            "CREATE INDEX IF NOT EXISTS idx_point " +
             "ON geometry_pointXYZ(profile, scenario);"
         )
 

src/Solver/CommandLine.py

@@ -158,11 +158,15 @@ class CommandLineSolver(AbstractSolver):
             txt = add_file.text
             txt = txt.replace("@version", get_version())
             txt = txt.replace("@date", datetime.now().isoformat(sep=" "))
+            txt = self._prepare_additional_file_text(add_file.path, txt)
 
             f.write(txt)
 
         return files
 
+    def _prepare_additional_file_text(self, path, text):
+        return text
+
     def export_study_description(self, study, repertory, qlog, name="0"):
         files = []
 

src/Solver/Mage.py

@@ -679,6 +679,23 @@ class Mage(CommandLineSolver):
 
         return res
 
+    def _prepare_additional_file_text(self, path, text):
+        if os.path.splitext(path)[1].lower() != ".var":
+            return text
+
+        lines = []
+        for line in text.splitlines(keepends=True):
+            stripped = line.lstrip()
+            fields = stripped.split()
+            if stripped.startswith("$"):
+                # MAGE only recognizes rule markers in the first column.
+                line = stripped
+            elif (not stripped.startswith("*") and len(fields) >= 2
+                  and fields[1] == "W2Z"):
+                line = "$ " + stripped
+            lines.append(line)
+        return "".join(lines)
+
     def _export_VAR(self, study, repertory, qlog, name="0"):
         files = []
 
@@ -729,9 +746,11 @@ class Mage(CommandLineSolver):
                         name = bhs.name
                         if name == "":
                             name = f"HS_{bhs.id:>3}".replace(" ", "0")
+                        else:
+                            name = name.replace(" ", "_")
 
                         f.write(
-                            f"${name} clapet"
+                            f"${name} clapet\n"
                         )
         return files
 

src/View/BoundaryCondition/Edit/Table.py

@@ -143,21 +143,43 @@ class TableModel(PamhyrTableModel):
         self.update()
 
     def delete(self, rows, parent=QModelIndex()):
-        self.beginRemoveRows(parent, rows[0], rows[-1])
+        if not rows:
+            return
 
-        rows = list(map(
-            lambda r: self.get_true_data_row(r),
-            rows))
+        rows = sorted(set(rows))
+        rows = [
+            self.get_true_data_row(row)
+            for row in rows
+        ]
 
+        # A multiple selection may contain non-contiguous rows. A model reset
+        # accurately represents this atomic undo command, unlike one
+        # beginRemoveRows() call spanning rows that are not all removed.
+        self.beginResetModel()
+        try:
             self._undo.push(
                 DelCommand(
                     self._data, rows
                 )
             )
-
-        self.endRemoveRows()
+        finally:
+            self.endResetModel()
         self.update()
 
+    def undo(self):
+        self.beginResetModel()
+        try:
+            self._undo.undo()
+        finally:
+            self.endResetModel()
+
+    def redo(self):
+        self.beginResetModel()
+        try:
+            self._undo.redo()
+        finally:
+            self.endResetModel()
+
     def sort(self, _reverse, parent=QModelIndex()):
         self.layoutAboutToBeChanged.emit()
 

src/View/BoundaryCondition/Edit/UndoCommand.py

@@ -97,12 +97,11 @@ class DelCommand(QUndoCommand):
         QUndoCommand.__init__(self)
 
         self._data = data
-        self._rows = rows
-
-        self._bc = []
-        for row in rows:
-            self._bc.append(self._data.get_i(row))
-        self._bc.sort()
+        self._rows = sorted(set(rows))
+        self._bc = [
+            self._data.get_i(row)
+            for row in self._rows
+        ]
 
     def undo(self):
         for el in self._bc:
@@ -120,13 +119,17 @@ class SortCommand(QUndoCommand):
         self._data = data
         self._reverse = _reverse
 
-        self._old = self._data.data
-        self._indexes = None
+        # Keep the order of the complete underlying list, including soft
+        # deleted items. ``data`` only exposes visible items and therefore
+        # cannot reliably restore a sort after a delete/undelete operation.
+        self._old_positions = {
+            id(item): index
+            for index, item in enumerate(self._data._data)
+        }
 
     def undo(self):
-        ll = self._data.data
         self._data.sort(
-            key=lambda x: self._indexes[ll.index(x)]
+            key=lambda item: self._old_positions[id(item)]
         )
 
     def redo(self):
@@ -134,14 +137,6 @@ class SortCommand(QUndoCommand):
             _reverse=self._reverse,
             key=lambda x: x[0]
         )
-        if self._indexes is None:
-            self._indexes = list(
-                map(
-                    lambda p: self._old.index(p),
-                    self._data.data
-                )
-            )
-            self._old = None
 
 
 class MoveCommand(QUndoCommand):

src/View/BoundaryCondition/UndoCommand.py

@@ -124,7 +124,7 @@ class DelCommand(QUndoCommand):
             self._bcs.insert(self._tab, row, el)
 
     def redo(self):
-        self._bcs.delete_i(self._tab, self._rows)
+        self._bcs.delete(self._tab, [el for _, el in self._bc])
 
 
 class SortCommand(QUndoCommand):

src/View/Geometry/Table.py

@@ -259,6 +259,9 @@ class GeometryReachTableModel(PamhyrTableModel):
                 **data
             )
 
+        self.apply_meshing(new_profiles, data, tableView)
+
+    def apply_meshing(self, new_profiles, data, tableView):
         if new_profiles is None:
             return
 

src/View/Geometry/Window.py

@@ -350,14 +350,22 @@ class GeometryWindow(PamhyrWindow):
             tableview = self.tableView
 
             def fn():
-                self._table.meshing(mesher, data, tableview)
-            # self._table.meshing(mesher, data, tableview)
+                return mesher.meshing(self._reach, **data)
+
             dlg2 = WaitingDialog(
                 payload_fn=fn,
                 title="waiting_mesh",
                 parent=self
             )
             dlg2.exec_()
+
+            if dlg2.worker_exception is not None:
+                raise dlg2.worker_exception
+
+            # Models, undo stacks and selections belong to the GUI thread.
+            self._table.apply_meshing(
+                dlg2.worker_result, data, tableview
+            )
             self._table.update()
 
         except Exception as e:

src/View/Results/CustomPlot/Plot.py

@@ -95,7 +95,7 @@ class CustomPlot(PamhyrPlot):
         table = results.get("table")["zfd"]
         ts = results.get_timestamp_id(self._current_timestamp)
 
-        rk = reach.geometry.get_rk()
+        # rk = reach.geometry.get_rk()
         z = list(
             map(
                 lambda p: table[ts, p.global_index],

src/View/Results/Window.py

@@ -1150,7 +1150,6 @@ class ResultsWindow(PamhyrWindow):
         z = table["Z"][:, profile.global_index]
         q = table["Q"][:, profile.global_index]
         v = table["V"][:, profile.global_index]
-        zfd = table["zfd"][:, profile.global_index]
 
         if self._current_results == 2:
             reach1 = self._results[0].river.reach(self._reach)
@@ -1169,8 +1168,8 @@ class ResultsWindow(PamhyrWindow):
             v2 = table["V"][:, profile2.global_index]
 
         if "bed_elevation" in y:
-            if reach.has_bedload():
-                z_min = zfd
+            if reach.has_bedload() and "zfd" in table:
+                z_min = table["zfd"][:, profile.global_index]
             else:
                 z_min = [profile.geometry.z_min()] * len(self._timestamps)
 

src/View/WaitingDialog.py

@@ -47,9 +47,15 @@ class Worker(QObject):
         super(self.__class__, self).__init__(parent)
 
         self._payload_fn = payload_fn
+        self.result = None
+        self.exception = None
 
     def process(self):
-        self._payload_fn()
+        try:
+            self.result = self._payload_fn()
+        except Exception as exception:
+            self.exception = exception
+        finally:
             self.signalStatus.emit('end')
 
 
@@ -123,6 +129,14 @@ class WaitingDialog(PamhyrDialog):
         self._timer = QTimer(self)
         self._timer.timeout.connect(self.update_spinner)
 
+    @property
+    def worker_result(self):
+        return self._worker.result
+
+    @property
+    def worker_exception(self):
+        return self._worker.exception
+
     def update_spinner(self):
         self._spinner_step += 1
         self._spinner_step %= len(self._spin)