src/Tools/blocFissure/gmu/quadranglesToShape.py

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