# 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, interp 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 # we never modify original profiles, we work on copies m_profiles = self.st_to_m(reach, limites) new_profiles = self.interpolate_transversal_step(m_profiles, step) for new_profiles2 in new_profiles: for p in new_profiles2: p.hard_purge_aligned_points() return new_profiles def st_to_m(self, reach, limites): gl = reach.compute_guidelines() profiles = [reach.profile(i).copy() for i in range(limites[0], limites[1]+1) ] # 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]) return profiles def interpolate_transversal_step(self, profiles, step): # calcul number of intermediate profiles np = [] for i in range(len(profiles)-1): np.append(int((profiles[i+1].rk - profiles[i].rk) / step) - 1) if np[-1] < 0: np[-1] = 0 new_profiles = [] for i in range(len(profiles)-1): new_profiles2 = [] d = 1.0/float(np[i]+1) for j in range(np[i]): p = profiles[i].copy() # RATIO entre les deux sections initiales dj = float(j+1)*d pt1 = profiles[i].points pt2 = profiles[i+1].points for k in range(len(pt1)): p.points[k].x = pt1[k].x + \ dj*(pt2[k].x - pt1[k].x) p.points[k].y = pt1[k].y + \ dj*(pt2[k].y - pt1[k].y) p.points[k].z = pt1[k].z + \ dj*(pt2[k].z - pt1[k].z) p.rk = profiles[i].rk + \ dj*(profiles[i+1].rk - profiles[i].rk) p.name = f'interpol{p.rk}' new_profiles2.append(p) new_profiles.append(new_profiles2) return new_profiles 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, begin_rk, end_rk, 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 [] nprof = reach.number_profiles if origin >= nprof or origin < 0: logger.warning( f"Origin outside reach" ) return reach.get_rk() # orientation is the orintation of the bief: # 1: up -> downstream # 2: down -> upstream # else: keep current orientation old_orientation = sign(reach.profiles[-1].rk - reach.profiles[0].rk) if orientation == 1: sgn = 1.0 elif orientation == 2: sgn = -1.0 else: sgn = old_orientation new_rk = [0.0]*nprof new_rk[origin] = origin_value frictions = reach._parent.frictions old_begin_rk = list(map(lambda x: x * old_orientation, begin_rk)) old_end_rk = list(map(lambda x: x * old_orientation, end_rk)) old_rk = list(map(lambda x: x * old_orientation, reach.get_rk())) 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])) new_rk[i] = new_rk[i-1] + (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])) new_rk[i] = new_rk[i+1] - (sgn * (d1 + d2) / 2) new_begin_rk = interp(old_begin_rk, old_rk, new_rk) new_end_rk = interp(old_end_rk, old_rk, new_rk) if abs(sgn - old_orientation) > 0.1: # orientation changed return new_rk, new_end_rk, new_begin_rk else: return new_rk, new_begin_rk, new_end_rk