diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 43bd3e9b..715364bd 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -60,10 +60,6 @@ dl-mage8-linux: artifacts: paths: - mage8-linux/mage - - mage8-linux/mage_as7 - - mage8-linux/mage_extraire - - mage8-linux/mailleurTT - - mage8-linux/libbief.so dl-mage8-windows: stage: downloads @@ -79,10 +75,6 @@ dl-mage8-windows: artifacts: paths: - mage8-windows/mage.exe - - mage8-windows/mage_as7.exe - - mage8-windows/mage_extraire.exe - - mage8-windows/mailleurTT.exe - - mage8-windows/libbief.dll dl-adists-linux: stage: downloads diff --git a/src/Meshing/Internal.py b/src/Meshing/Internal.py new file mode 100644 index 00000000..ed9a5c48 --- /dev/null +++ b/src/Meshing/Internal.py @@ -0,0 +1,331 @@ +# Internal.py -- Pamhyr +# Copyright (C) 2023-2025 INRAE +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# -*- coding: utf-8 -*- + +import logging +from numpy import sign +from copy import deepcopy + +from tools import logger_color_red, logger_color_reset, logger_exception +from Meshing.AMeshingTool import AMeshingTool + +from Model.Geometry.ProfileXYZ import ProfileXYZ + +logger = logging.getLogger() + + +class InternalMeshing(AMeshingTool): + def __init__(self): + super(InternalMeshing, self).__init__() + + ########### + # Meshing # + ########### + + def meshing(self, reach, + step: float = 50, + limites=[-1, -1], + linear: bool = True): + if reach is None or len(reach.profiles) == 0: + return reach + + if limites[0] > limites[1]: + lim = limites[1] + limites[1] = limites[0] + limites[0] = lim + elif limites[0] == limites[1]: + return reach + + self.st_to_m(reach, limites) + self.interpolate_transversal_step(reach, + limites, + step) + self.purge_aligned_points(reach) + + return reach + + def st_to_m(self, reach, limites): + + gl = reach.compute_guidelines() + profiles = reach.profiles[limites[0]:limites[1]+1] + + print(profiles) + + # we make sure that the lines are in the left-to-right order + guide_list = [ + x.name + for x in profiles[0].named_points() + if x.name in gl[0] + ] + + gl1 = [""] + guide_list + gl2 = guide_list + [""] + max_values = [0] * (len(guide_list) + 1) + max_values_index = [0] * (len(guide_list) + 1) + # we work between each couple of guidelines + for j, p in enumerate(profiles): + index1 = 0 + index2 = p.named_point_index(guide_list[0]) + if index2 - index1 > max_values[0]: + max_values[0] = index2 - index1 + max_values_index[0] = j + for i in range(len(guide_list) - 1): + index1 = p.named_point_index(guide_list[i]) + index2 = p.named_point_index(guide_list[i+1]) + if index2 - index1 > max_values[i + 1]: + max_values[i + 1] = index2 - index1 + max_values_index[i + 1] = j + index1 = p.named_point_index(guide_list[-1]) + index2 = len(p) - 1 + if index2 - index1 > max_values[-1]: + max_values[-1] = index2 - index1 + max_values_index[-1] = j + + for i in range(len(max_values)): + for isect in range(max_values_index[i], len(profiles)-1): + self.compl_sect(profiles[isect], + profiles[isect+1], + gl1[i], gl2[i]) + for isect in reversed(range(0, max_values_index[i])): + self.compl_sect(profiles[isect+1], + profiles[isect], + gl1[i], gl2[i]) + + def interpolate_transversal_step(self, + reach, + limites, + step): + # calcul number of intermediate profiles + np = [] + for i in range(limites[0], limites[1]): + np.append(int((reach.profiles[i+1].rk - + reach.profiles[i].rk) / step) - 1) + if np[-1] < 0: + np[-1] = 0 + + d = [] # ratios + p = [] # new profiles + ptr = int(limites[0]) + for i in range(limites[1] - limites[0]): + d = 1.0/float(np[i]+1) + ptr0 = ptr + for j in range(np[i]): + p = reach.profiles[ptr0].copy() + # RATIO entre les deux sections initiales + dj = float(j+1)*d + ptr += 1 # next profile, original + for k in range(len(reach.profiles[ptr0].points)): + p.points[k].x = reach.profiles[ptr0].points[k].x + \ + dj*(reach.profiles[ptr].points[k].x - + reach.profiles[ptr0].points[k].x) + p.points[k].y = reach.profiles[ptr0].points[k].y + \ + dj*(reach.profiles[ptr].points[k].y - + reach.profiles[ptr0].points[k].y) + p.points[k].z = reach.profiles[ptr0].points[k].z + \ + dj*(reach.profiles[ptr].points[k].z - + reach.profiles[ptr0].points[k].z) + p.rk = reach.profiles[ptr0].rk + \ + dj*(reach.profiles[ptr].rk - + reach.profiles[ptr0].rk) + p.name = f'interpol{p.rk}' + reach.insert_profile(ptr, p) + ptr += 1 # next profile, original + + def compl_sect(self, sect1, sect2, tag1, tag2): + # sect1: reference section + # sect2: section to complete + # start1 and 2: indices of the first point of the bed to complete + # end = start + len + # len1 and 2: number of intervals + # len1: target len + + if tag1 == '': # left end point + start1 = 0 + start2 = 0 + else: + start1 = sect1.named_point_index(tag1) + start2 = sect2.named_point_index(tag1) + + if tag2 == '': # right end point + end1 = sect1.nb_points-1 + end2 = sect2.nb_points-1 + else: + end1 = sect1.named_point_index(tag2) + end2 = sect2.named_point_index(tag2) + + len1 = end1 - start1 + len2 = end2 - start2 + + if len1 < len2: + self.compl_sect(sect2, sect1, start2, start1, + len2, len1, tag1, tag2) + return + elif len1 == len2: + return + + ltot1 = 0.0 + ltot2 = 0.0 + alpha = [] + beta = [] + # remove names + sect2.points[start2].name = '' + sect2.points[end2].name = '' + for i in range(len1): + ltot1 += sect1.point(start1+i).dist(sect1.point(start1+i+1)) + alpha.append(ltot1) + alpha = list(map(lambda x: x/ltot1, alpha)) # target ratios + for i in range(len2): + ltot2 += sect2.point(start2+i).dist(sect2.point(start2+i+1)) + beta.append(ltot2) + 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] + if len2 < 1 or ltot2 < 0.0001: + # unique point (copy len1 times) + # we are at an edge of the profile + for i in range(len1-len2): + p = sect2.point(start2).copy() + sect2.insert_point(start2, p) + if tag1 == '': # left end point + sect2.point(start2).name = '' + if tag2 == '': # left end point + sect2.point(start2+1).name = '' + elif ltot1 < 0.0001: + sect2.add_npoints(len1-len2) + else: # regular case + # 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 + + 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: + print(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 = '' + sect2.point(start2).name = tag1 + sect2.point(start2+len2).name = tag2 + sect2.modified() + + def update_rk(self, reach, + step: float = 50, + limites=[-1, -1], + origin=0, + directrices=['un', 'np'], + lplan: bool = False, + lm: int = 3, + linear: bool = False, + origin_value=0.0, + orientation=0): + if reach is None or len(reach.profiles) == 0: + return reach + + nprof = reach.number_profiles + if origin >= nprof or origin < 0: + logger.warning( + f"Origin outside reach" + ) + return reach + + # orientation is the orintation of the bief: + # 1: up -> downstream + # 2: down -> upstream + # else: keep current orientation + if orientation == 1: + sgn = 1.0 + elif orientation == 2: + sgn = -1.0 + else: + sgn = sign(reach.profiles[-1].rk - reach.profiles[0].rk) + + reach.profiles[origin].rk = origin_value + + if origin < nprof - 1: + for i in range(origin+1, nprof): + # 2D + d1 = reach.profiles[i-1].named_point( + directrices[0] + ).dist_2d(reach.profiles[i].named_point(directrices[0])) + d2 = reach.profiles[i-1].named_point( + directrices[1] + ).dist_2d(reach.profiles[i].named_point(directrices[1])) + reach.profiles[i].rk = reach.profiles[i-1].rk \ + + (sgn * (d1 + d2) / 2) + if origin > 0: + for i in reversed(range(0, origin)): + # 2D + d1 = reach.profiles[i+1].named_point( + directrices[0] + ).dist_2d(reach.profiles[i].named_point(directrices[0])) + d2 = reach.profiles[i+1].named_point( + directrices[1] + ).dist_2d(reach.profiles[i].named_point(directrices[1])) + reach.profiles[i].rk = reach.profiles[i+1].rk \ + - (sgn * (d1 + d2) / 2) + + def purge_aligned_points(self, reach): + for p in reach.profiles: + p.purge_aligned_points() diff --git a/src/Meshing/Mage.py b/src/Meshing/Mage.py index d02122bd..94ef5385 100644 --- a/src/Meshing/Mage.py +++ b/src/Meshing/Mage.py @@ -19,6 +19,7 @@ import os import logging import tempfile +from numpy import sign from ctypes import ( cdll, diff --git a/src/Model/Geometry/PointXYZ.py b/src/Model/Geometry/PointXYZ.py index e1bea68c..b82e51e9 100644 --- a/src/Model/Geometry/PointXYZ.py +++ b/src/Model/Geometry/PointXYZ.py @@ -394,3 +394,25 @@ class PointXYZ(Point): ) ** 0.5 return res + + def dist_from_seg(self, p1, p2): + a2 = pow(PointXYZ.distance(self, p1), 2) + b2 = pow(PointXYZ.distance(self, p2), 2) + c2 = pow(PointXYZ.distance(p1, p2), 2) + if c2 < 0.00000001: + # si les deux points sont confondus + d = 0.0 + else: + d = np.sqrt(abs(a2 - pow((c2 - b2 + a2) / (2*np.sqrt(c2)), 2))) + + return d + + def copy(self): + p = PointXYZ(self.x, + self.y, + self.z, + self.name, + self._profile, + self._status) + p.sl = self.sl + return p diff --git a/src/Model/Geometry/Profile.py b/src/Model/Geometry/Profile.py index af1e8f5e..04310d8c 100644 --- a/src/Model/Geometry/Profile.py +++ b/src/Model/Geometry/Profile.py @@ -75,6 +75,14 @@ class Profile(object): def point(self, index): return self.points[index] + def named_point(self, name): + return next((p for p in self.points if p.name == name), None) + + def named_point_index(self, name): + return next( + (p for p in enumerate(self.points) if p[1].name == name), None + )[0] + @property def reach(self): return self._reach @@ -193,7 +201,7 @@ class Profile(object): """Insert point at index. Args: - index: The index of new profile. + index: The index of new point. point: The point. Returns: diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py index 2cab60b4..6f5f05fd 100644 --- a/src/Model/Geometry/ProfileXYZ.py +++ b/src/Model/Geometry/ProfileXYZ.py @@ -1087,3 +1087,50 @@ class ProfileXYZ(Profile, SQLSubModel): p.z = p.z + z self.modified() + + def modified(self): + self.tab_up_to_date = False + self.station_up_to_date = False + + def add_npoints(self, npoints): + # add npoints in a profile + for k in range(npoints): + disti = 0.0 + i = 1 + for j in range(len(self.points)-1): + distj = self.point(j).dist(self.point(j+1)) + if distj > disti: + disti = distj + i = j + self.insert_point(i, self.point(i).copy()) + self.point(i+1).name = '' + self.point(i+1).x = 0.5 * self.point(i).x + 0.5 * self.point(i+2).x + self.point(i+1).y = 0.5 * self.point(i).y + 0.5 * self.point(i+2).y + self.point(i+1).z = 0.5 * self.point(i).z + 0.5 * self.point(i+2).z + + def copy(self): + p = ProfileXYZ(self.id, + self.name, + self.rk, + self.reach, + self.num, + 0, + 0, 0, + self._status) + for i, k in enumerate(self.points): + p.insert_point(i, k.copy()) + + return p + + def purge_aligned_points(self): + + align = True + while (align): + align = False + for i in range(1, self.number_points - 1): + d = self.point(i).dist_from_seg(self.point(i-1), + self.point(i+1)) + if d < 0.00001 and self.point(i).name == "": + align = True + self.delete_i([i]) + break diff --git a/src/Model/Results/River/River.py b/src/Model/Results/River/River.py index 20a14613..b94b06d6 100644 --- a/src/Model/Results/River/River.py +++ b/src/Model/Results/River/River.py @@ -280,6 +280,12 @@ class Reach(SQLSubModel): ) ) + self._profile_mask = list( + map( + lambda p: p.name[0:8] != 'interpol', self._profiles + ) + ) + def __len__(self): return len(self._profiles) @@ -295,6 +301,10 @@ class Reach(SQLSubModel): def profiles(self): return self._profiles.copy() + @property + def profile_mask(self): + return self._profile_mask + def profile(self, id): return self._profiles[id] diff --git a/src/View/Geometry/MeshingDialog.py b/src/View/Geometry/MeshingDialog.py index e2bf0310..4ecb1862 100644 --- a/src/View/Geometry/MeshingDialog.py +++ b/src/View/Geometry/MeshingDialog.py @@ -23,12 +23,12 @@ from PyQt5.QtGui import ( ) from PyQt5.QtCore import ( - Qt, QVariant, QAbstractTableModel, + Qt, QVariant, QAbstractTableModel, QItemSelectionModel, ) from PyQt5.QtWidgets import ( QDialogButtonBox, QComboBox, QUndoStack, QShortcut, - QDoubleSpinBox, + QDoubleSpinBox, QRadioButton, ) @@ -55,7 +55,7 @@ class MeshingDialog(PamhyrDialog): self._space_step = 50.0 self._lplan = False self._lm = "3" - self._linear = False + self._linear = True self._begin_cs = -1 self._end_cs = -1 self._begin_dir = "un" @@ -71,7 +71,7 @@ class MeshingDialog(PamhyrDialog): self._origin = self._reach.profile(0) self._init_default_values_profiles() - self._init_default_values_guidelines() + # self._init_default_values_guidelines() self.set_double_spin_box( "doubleSpinBox_space_step", @@ -83,14 +83,24 @@ class MeshingDialog(PamhyrDialog): else: self.set_radio_button("radioButton_spline", True) + self.find(QRadioButton, "radioButton_spline").setEnabled(False) + self.find(QRadioButton, "radioButton_linear").setEnabled(False) + def _init_default_values_profiles(self): profiles = self.profiles self.combobox_add_items("comboBox_begin_rk", profiles) self.combobox_add_items("comboBox_end_rk", profiles) - self.set_combobox_text("comboBox_begin_rk", profiles[0]) - self.set_combobox_text("comboBox_end_rk", profiles[-1]) + r = self.parent.tableView\ + .selectionModel()\ + .selectedRows() + if len(r) <= 1: + self.set_combobox_text("comboBox_begin_rk", profiles[0]) + self.set_combobox_text("comboBox_end_rk", profiles[-1]) + else: + self.set_combobox_text("comboBox_begin_rk", profiles[r[0].row()]) + self.set_combobox_text("comboBox_end_rk", profiles[r[-1].row()]) @property def profiles(self): @@ -167,8 +177,8 @@ class MeshingDialog(PamhyrDialog): self._begin_cs = self.profiles.index(p1) self._end_cs = self.profiles.index(p2) - self._begin_dir = self.get_combobox_text("comboBox_begin_gl") - self._end_dir = self.get_combobox_text("comboBox_end_gl") + # self._begin_dir = self.get_combobox_text("comboBox_begin_gl") + # self._end_dir = self.get_combobox_text("comboBox_end_gl") super().accept() diff --git a/src/View/Geometry/UndoCommand.py b/src/View/Geometry/UndoCommand.py index c0c9a7fc..c2de3232 100644 --- a/src/View/Geometry/UndoCommand.py +++ b/src/View/Geometry/UndoCommand.py @@ -28,8 +28,6 @@ from PyQt5.QtWidgets import ( from Model.Geometry import Reach from Model.Except import exception_message_box -from Meshing.Mage import MeshingWithMage - logger = logging.getLogger() @@ -238,7 +236,7 @@ class MeshingCommand(QUndoCommand): self._mesher = mesher self._command = command - self._profiles = reach.profiles.copy() + self._profiles = [p.copy() for p in reach.profiles] self._profiles.reverse() self._new_profiles = None @@ -262,7 +260,7 @@ class MeshingCommand(QUndoCommand): **self._data ) - self._new_profiles = self._reach.profiles.copy() + self._new_profiles = [p.copy() for p in self._reach.profiles] self._new_profiles.reverse() else: self._reach.purge() diff --git a/src/View/Geometry/Window.py b/src/View/Geometry/Window.py index 30eec639..e3e1ac9b 100644 --- a/src/View/Geometry/Window.py +++ b/src/View/Geometry/Window.py @@ -32,6 +32,7 @@ from PyQt5.QtGui import ( from PyQt5.QtCore import ( QModelIndex, Qt, QSettings, pyqtSlot, QItemSelectionModel, QCoreApplication, QSize, + QItemSelection, QItemSelectionRange, ) from PyQt5.QtWidgets import ( QApplication, QMainWindow, QFileDialog, QCheckBox, @@ -49,9 +50,7 @@ from View.Tools.Plot.PamhyrCanvas import MplCanvas from View.SelectReach.Window import SelectReachWindow -from Meshing.Mage import ( - MeshingWithMage, MeshingWithMageMailleurTT -) +from Meshing.Internal import InternalMeshing from View.Geometry.Table import GeometryReachTableModel from View.Geometry.PlotXY import PlotXY @@ -300,6 +299,13 @@ class GeometryWindow(PamhyrWindow): self.tableView.model().blockSignals(False) def edit_meshing(self): + + rows = list( + set( + (i.row() for i in self.tableView.selectedIndexes()) + ) + ) + selected_rk = [self._reach.profile(r).rk for r in rows] try: dlg = MeshingDialog( reach=self._reach, @@ -310,8 +316,6 @@ class GeometryWindow(PamhyrWindow): data = { "step": dlg.space_step, "limites": [dlg.begin_cs, dlg.end_cs], - "directrices": [dlg.begin_dir, dlg.end_dir], - "lplan": dlg.lplan, "linear": dlg.linear, } self._edit_meshing(data) @@ -319,18 +323,31 @@ class GeometryWindow(PamhyrWindow): logger_exception(e) return + ind = [] + for i in range(self._reach.number_profiles): + if self._reach.profile(i).rk in selected_rk: + ind.append(i) + self.tableView.setFocus() + selection = self.tableView.selectionModel() + index = QItemSelection() + if len(ind) > 0: + for i in ind: + index.append(QItemSelectionRange( + self.tableView.model().index(i, 0) + )) + selection.select( + index, + QItemSelectionModel.Rows | + QItemSelectionModel.ClearAndSelect | + QItemSelectionModel.Select + ) + def _edit_meshing(self, data): try: - mesher = MeshingWithMageMailleurTT() + mesher = InternalMeshing() self._table.meshing(mesher, data) except Exception as e: logger_exception(e) - raise ExternFileMissingError( - module="mage", - filename="MailleurTT", - path=MeshingWithMageMailleurTT._path(), - src_except=e - ) def update_rk(self): try: diff --git a/src/View/Results/Table.py b/src/View/Results/Table.py index 45482c38..1b30de26 100644 --- a/src/View/Results/Table.py +++ b/src/View/Results/Table.py @@ -23,6 +23,8 @@ from numpy import sqrt from tools import timer, trace +from itertools import compress + from PyQt5.QtGui import ( QKeySequence, QColor ) @@ -51,8 +53,12 @@ class TableModel(PamhyrTableModel): self._lst = _river.reachs elif self._opt_data == "profile": self._lst = _river.reach(0).profiles + # self._lst = list(filter(lambda x: x.name[0:8] != 'interpol', + # _river.reach(0).profiles)) elif self._opt_data == "raw_data": self._lst = _river.reach(0).profiles + # self._lst = list(filter(lambda x: x.name[0:8] != 'interpol', + # _river.reach(0).profiles)) elif self._opt_data == "solver": self._lst = self._parent._solvers @@ -152,7 +158,9 @@ class TableModel(PamhyrTableModel): if self._opt_data == "reach": self._lst = _river.reachs elif self._opt_data == "profile" or self._opt_data == "raw_data": - self._lst = _river.reach(reach).profiles + # self._lst = _river.reach(reach).profiles + self._lst = list(compress(_river.reach(reach).profiles, + _river.reach(reach).profile_mask)) elif self._opt_data == "solver": self._lst = self._parent._solvers diff --git a/src/View/ui/MeshingOptions.ui b/src/View/ui/MeshingOptions.ui index de48342a..b7eab409 100644 --- a/src/View/ui/MeshingOptions.ui +++ b/src/View/ui/MeshingOptions.ui @@ -7,24 +7,14 @@ 0 0 520 - 341 + 245 Dialog - - - - Qt::Horizontal - - - QDialogButtonBox::Cancel|QDialogButtonBox::Ok - - - - + Parameters @@ -79,7 +69,7 @@ Spline - true + false @@ -88,6 +78,9 @@ Linear + + true + @@ -108,47 +101,14 @@ - - - - Guideline used for distance computation + + + + Qt::Horizontal + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok - - - - - true - - - Second guideline - - - - - - - true - - - First guideline - - - - - - - true - - - - - - - true - - - - diff --git a/tests_cases/Hogneau/hogneau.pamhyr b/tests_cases/Hogneau/hogneau.pamhyr index 9680f261..5488aa32 100644 Binary files a/tests_cases/Hogneau/hogneau.pamhyr and b/tests_cases/Hogneau/hogneau.pamhyr differ