src/Tools/blocFissure/gmu/quadranglesToShapeNoCorner.py

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