src/Tools/blocFissure/gmu/quadranglesToShapeNoCorner.py

   1 # -*- coding: utf-8 -*-
   2 # Copyright (C) 2014-2021  EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20 """Groupe de quadrangles de face transformé en face géométrique par filling"""
  21
  22 import logging
  23 import numpy as np
  24
  25 import GEOM
  26
  27 from .geomsmesh import geompy
  28 from .geomsmesh import geomPublish
  29
  30 from . import initLog
  31
  32 def mydot(a):
  33   """produit scalaire"""
  34   return np.dot(a,a)
  35
  36 def quadranglesToShapeNoCorner(meshQuad, shapeFissureParams, centreFondFiss):
  37   """
  38   groupe de quadrangles de face transformée en faces géométriques par filling
  39   on part de quadrangles définissant une zone a 4 cotés (convexe), et on reconstitue n lignes de p points.
  40   Ces n lignes de p points sont transformées en n courbes géométriques,
  41   à partir desquelles on reconstitue une surface géométrique.
  42   Il peut y avoir plusieurs faces géométriques reconstituées, si on fournit des groupes de quadrangles non connexes.
  43   On détecte les angles vifs, pour conserver des arêtes vives délimitant des faces connexes.
  44   @param meshQuad : maillages constitué de quadrangles constituant une ou plusieurs zones convexes
  45   @return (fillings, noeuds_Bords) : liste de geomObject, listes des bords (bord = liste ordonnée de noeuds (geomObject))
  46   """
  47   logging.info("start")
  48
  49   isVecteurDefaut = False
  50   if 'vecteurDefaut' in shapeFissureParams:
  51     isVecteurDefaut = True
  52     vecteurDefaut = shapeFissureParams['vecteurDefaut']
  53
  54   fillings = list()       # les faces reconstituées, découpées selon les arêtes vives
  55   noeuds_bords = list()   #
  56   bords_Partages = list() # contient a la fin les courbes correspondant aux arêtes vives
  57   fillconts = list()      # les faces reconstituées, sans découpage selon les arêtes vives
  58   idFilToCont = list()    # index face découpée vers face sans découpe
  59   iface = 0           # index face découpée
  60   icont = 0           # index face continue
  61   pisur2 = np.pi/2.0
  62   pisur4 = np.pi/4.0
  63
  64   allNodeIds = meshQuad.GetNodesId()
  65   while len(allNodeIds):
  66     logging.debug("len(allNodeIds): %s ", len(allNodeIds))
  67     nodeIds = allNodeIds
  68     for idNode in nodeIds: # rechercher un coin
  69       elems = meshQuad.GetNodeInverseElements(idNode)
  70       if ( len(elems) == 1 ):
  71         # un coin: un noeud, un element quadrangle
  72         idStart = idNode # le noeud de coin
  73         elemStart = elems[0] # l'élément quadrangle au coin
  74         break
  75     xyz = meshQuad.GetNodeXYZ(idStart)
  76     logging.debug("idStart %s, coords %s", idStart, str(xyz))
  77
  78     nodelines = list() # on va constituer une liste de lignes de points
  79     nextLine = True
  80     ligneFinale = False
  81     while nextLine:
  82       logging.debug("--- une ligne")
  83       idNode = idStart
  84       elem = elemStart
  85       if ligneFinale:
  86         agauche = False      # sens de parcours des 4 noeuds d'un quadrangle
  87         nextLine = False
  88       else:
  89         agauche = True
  90       ligneIncomplete = True # on commence une ligne de points
  91       debutLigne = True
  92       nodeline = list()
  93       elemline = list()
  94       while ligneIncomplete: # compléter la ligne de points
  95         nodeline.append(idNode)
  96         allNodeIds.remove(idNode)
  97         elemline.append(elem)
  98         nodes = meshQuad.GetElemNodes(elem)
  99         i = nodes.index(idNode) # repérer l'index du noeud courant (i) dans l'élément quadrangle (0 a 3)
 100         if agauche:             # déterminer le noeud suivant (j) et celui opposé (k) dans le quadrangle
 101           if i < 3:
 102             j = i+1
 103           else:
 104             j = 0
 105           if j < 3:
 106             k = j+1
 107           else:
 108             k = 0
 109         else:
 110           if i > 0:
 111             j = i -1
 112           else:
 113             j = 3
 114           if j > 0:
 115             k = j -1
 116           else:
 117             k = 3
 118         isuiv = nodes[j]   #noeud suivant
 119         iapres = nodes[k]  #noeud opposé
 120         if debutLigne:
 121           debutLigne = False
 122           # précédent a trouver, dernière ligne : précédent au lieu de suivant
 123           if agauche:
 124             if i > 0:
 125               iprec = nodes[i -1]
 126             else:
 127               iprec = nodes[3]
 128             idStart = iprec
 129             elems3 = meshQuad.GetNodeInverseElements(iprec)
 130             if len(elems3) == 1: # autre coin
 131               ligneFinale = True
 132             else:
 133               for elem3 in elems3:
 134                 if elem3 != elem:
 135                   elemStart = elem3
 136                   break
 137         #print nodes, idNode, isuiv, iapres
 138         elems1 = meshQuad.GetNodeInverseElements(isuiv)
 139         elems2 = meshQuad.GetNodeInverseElements(iapres)
 140         ligneIncomplete = False
 141         for elem2 in elems2:
 142           if elems1.count(elem2) and elem2 != elem:
 143             ligneIncomplete = True
 144             idNode = isuiv
 145             elem = elem2
 146             break
 147         if not  ligneIncomplete:
 148           nodeline.append(isuiv)
 149           allNodeIds.remove(isuiv)
 150       logging.debug("nodeline %s", nodeline)
 151       logging.debug("elemline %s", elemline)
 152       nodelines.append(nodeline)
 153     logging.debug("nodelines = %s", nodelines)
 154     longueur = [len(val) for val in nodelines]
 155     logging.debug("longueur = %s", longueur)
 156     # on a constitué une liste de lignes de points connexes
 157     logging.debug("dimensions [%s, %s]", len(nodelines),  len(nodeline))
 158
 159     # stockage des coordonnées dans un tableau numpy
 160     mat = np.zeros((len(nodelines), len(nodeline), 3))
 161     for i, ligne in enumerate(nodelines):
 162       for j, nodeId in enumerate(ligne):
 163         mat[i,j] = meshQuad.GetNodeXYZ(nodeId)
 164     logging.debug("matrice de coordonnées: \n%s",mat)
 165     logging.debug("dimensions %s", mat.shape)
 166
 167     # recherche d'angles supérieurs a un seuil sur une ligne : angle entre deux vecteurs successifs
 168     cosmin = np.cos(pisur4)                # TODO: angle reference en paramètre
 169     vecx = mat[:, 1:,  :] - mat[:, :-1, :] # vecteurs selon direction "x"
 170     vx0 = vecx[:, :-1, :]                  # vecteurs amont
 171     vx1 = vecx[:, 1:,  :]                  # vecteurs aval
 172     e = np.einsum('ijk,ijk->ij', vx0, vx1) # produit scalaire des vecteurs
 173     f = np.apply_along_axis(mydot, 2, vx0) # normes carrées vecteurs amont
 174     g = np.apply_along_axis(mydot, 2, vx1) # normes carrées vecteurs aval
 175     h = e/(np.sqrt(f*g))                   # cosinus
 176     ruptureX = h < cosmin                  # True si angle > reference
 177     logging.debug("matrice de rupture X: \n%s",ruptureX)
 178     rupX = [x for x in range(len(nodeline)-2) if np.prod(ruptureX[:,x])]
 179     logging.debug("colonnes de rupture: %s",rupX)
 180     # recherche d'angles supérieurs a un seuil sur une colonne : angle entre deux vecteurs successifs
 181     vecy = mat[ 1:, :, :] - mat[:-1, :, :] # vecteurs selon direction "y"
 182     vy0 = vecy[:-1, :, :]                  # vecteurs amont
 183     vy1 = vecy[ 1:, :, :]                  # vecteurs aval
 184     e = np.einsum('ijk,ijk->ij', vy0, vy1) # produit scalaire des vecteurs
 185     f = np.apply_along_axis(mydot, 2, vy0) # normes carrées vecteurs amont
 186     g = np.apply_along_axis(mydot, 2, vy1) # normes carrées vecteurs aval
 187     h = e/(np.sqrt(f*g))                   # cosinus
 188     ruptureY = h < cosmin                  # True si angle > reference
 189     logging.debug("matrice de rupture Y: \n%s",ruptureY)
 190     rupY = [x for x in range(len(nodelines)-2) if np.prod(ruptureY[x, :])]
 191     logging.debug("lignes de rupture: %s",rupY)
 192     if (len(rupX)*len(rupY)) > 0:
 193       logging.critical("""Cas non traité: présence d'angles vifs dans 2 directions,
 194       lors de la reconstitution des faces géométriques dans la zone remaillée""")
 195
 196     mats = list()
 197     bordsPartages = list()
 198     if (len(rupX)> 0):
 199       rupX.append(mat.shape[1]-1)
 200       for i, index in enumerate(rupX):
 201         imax = index+2
 202         imin = 0
 203         if i > 0:
 204           imin = rupX[i-1] + 1
 205         mats.append(mat[:, imin:imax, :])
 206         if imax == mat.shape[1] + 1:
 207           ifin = 0
 208         else:
 209           ifin = imax
 210         bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
 211     elif (len(rupY)> 0):
 212       rupY.append(mat.shape[0]-1)
 213       for i, index in enumerate(rupY):
 214         imax = index+2
 215         imin = 0
 216         if i > 0:
 217           imin = rupY[i-1] + 1
 218         mats.append(mat[imin:imax, :, :])
 219         if imax == mat.shape[0] + 1:
 220           ifin = 0
 221         else:
 222           ifin = imax
 223         bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
 224     else:
 225       mats.append(mat)
 226       bordsPartages.append([0,0])         # les indices différents de 0 correspondent à des bords partagés
 227
 228     curvconts = list()
 229     for nmat, amat in enumerate(mats):
 230       logging.debug("dimensions matrice %s: %s", nmat, amat.shape)
 231       nbLignes = amat.shape[1] # pas de rupture, ou rupture selon des colonnes: on transpose
 232       nbCols = amat.shape[0]
 233       if len(rupY) > 0 :       # rupture selon des lignes: pas de transposition
 234         nbLignes = amat.shape[0]
 235         nbCols = amat.shape[1]
 236       curves = list()
 237       noeudsBords = list()
 238       for i in range(4):
 239         noeudsBords.append([])
 240       k = 0
 241       for i in range(nbLignes):
 242         nodeList = list()
 243         for j in range(nbCols):
 244           #logging.debug("point[%s,%s] = (%s, %s, %s)",i,j,amat[i,j,0], amat[i,j,1], amat[i,j,2])
 245           if len(rupY) > 0 : # pas de transposition
 246             node = geompy.MakeVertex(amat[i,j,0], amat[i,j,1], amat[i,j,2])
 247           else:              # transposition
 248             node = geompy.MakeVertex(amat[j,i,0], amat[j,i,1], amat[j,i,2])
 249           nodeList.append(node)
 250           if i == 0:
 251             noeudsBords[0].append(node)
 252             #name = "bord0_%d"%k
 253             #geomPublish(initLog.debug,  node, name )
 254           if i == (nbLignes -1):
 255             noeudsBords[2].append(node)
 256             #name = "bord2_%d"%k
 257             #geomPublish(initLog.debug,  node, name )
 258           if j == 0:
 259             noeudsBords[1].append(node)
 260             #name = "bord1_%d"%k
 261             #geomPublish(initLog.debug,  node, name )
 262           if j == (nbCols -1):
 263             noeudsBords[3].append(node)
 264             #name = "bord3_%d"%k
 265             #geomPublish(initLog.debug,  node, name )
 266             k += 1
 267         curve = geompy.MakeInterpol(nodeList, False, False)
 268         #name = "curve_%d"%i
 269         #geomPublish(initLog.debug,  curve, name )
 270         if len(curvconts) == 0 or len(curves) > 0: # éliminer les doublons de la surface sans découpe
 271           curvconts.append(nodeList)
 272         curves.append(curve)
 273       if bordsPartages[nmat][0] :
 274         bordsPartages[nmat][0] = curves[0]  # la première ligne est un bord partagé
 275       else:
 276         bordsPartages[nmat][0] = None
 277       if bordsPartages[nmat][1] :
 278         bordsPartages[nmat][1] = curves[-1] # la dernière ligne est un bord partagé
 279       else:
 280         bordsPartages[nmat][1] = None
 281       filling = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
 282       # --- test orientation filling
 283       vertex = geompy.MakeVertexOnSurface(filling, 0.5, 0.5)
 284       normal = geompy.GetNormal(filling, vertex)
 285
 286       if centreFondFiss is not None:
 287         logging.debug("orientation filling a l'aide du centre de fond de fissure")
 288         vecteurDefaut = geompy.MakeVector(centreFondFiss, vertex)
 289
 290       if not isVecteurDefaut:
 291         pointIn_x = 0.0
 292         pointIn_y = 0.0
 293         pointIn_z = 0.0
 294         pointExplicite = False
 295         if 'pointIn_x' in shapeFissureParams:
 296           pointExplicite = True
 297           pointIn_x = shapeFissureParams['pointIn_x']
 298         if 'pointIn_y' in shapeFissureParams:
 299           pointExplicite = True
 300           pointIn_y = shapeFissureParams['pointIn_y']
 301         if 'pointIn_z' in shapeFissureParams:
 302           pointExplicite = True
 303           pointIn_z = shapeFissureParams['pointIn_z']
 304         if pointExplicite:
 305           cdg = geompy.MakeVertex(pointIn_x, pointIn_y, pointIn_z)
 306           logging.debug("orientation filling par point intérieur %s", (pointIn_x, pointIn_y, pointIn_z))
 307           vecteurDefaut = geompy.MakeVector(cdg, vertex)
 308
 309       if 'convexe' in shapeFissureParams:
 310         isConvexe = shapeFissureParams['convexe']
 311         logging.debug("orientation filling par indication de convexité %s", isConvexe)
 312         cdg = geompy.MakeCDG(filling)
 313         if isConvexe:
 314           vecteurDefaut = geompy.MakeVector(cdg, vertex)
 315         else:
 316           vecteurDefaut = geompy.MakeVector(vertex, cdg)
 317
 318       if vecteurDefaut is not None:
 319         geomPublish(initLog.debug, normal, "normFillOrig%d"%iface)
 320         geomPublish(initLog.debug, vecteurDefaut, "fromInterieur%d"%iface)
 321         if ( geompy.GetAngleRadians(vecteurDefaut, normal) > pisur2 ):
 322           filling = geompy.ChangeOrientation(filling)
 323       geomPublish(initLog.debug,  filling, "filling%d"%iface )
 324       #geompy.ExportBREP(filling, "filling.brep")
 325       iface = iface+1
 326       fillings.append(filling)
 327       noeuds_bords.append(noeudsBords)
 328       idFilToCont.append(icont)
 329       bords_Partages += bordsPartages
 330       logging.debug("bords_Partages = %s", bords_Partages)
 331     # --- loop on mats
 332     # --- reconstruction des faces continues à partir des listes de noeuds
 333     #     les courbes doivent suivre la courbure pour éviter les oscillations
 334     if icont == iface - 1: # pas de découpe, on garde la même face
 335       fillcont = fillings[-1]
 336     else:
 337       nbLignes = len(curvconts[0])
 338       curves = list()
 339       for i in range(nbLignes):
 340         nodes = [curvconts[j][i] for j in range(len(curvconts))]
 341         curve = geompy.MakeInterpol(nodes, False, False)
 342         curves.append(curve)
 343       fillcont = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
 344     geomPublish(initLog.debug,  fillcont, "filcont%d"%icont )
 345     fillconts.append(fillcont)
 346     icont = icont+1
 347     # --- loop while there are remaining nodes
 348
 349   return fillings, noeuds_bords, bords_Partages, fillconts, idFilToCont