src/Tools/blocFissure/gmu/quadranglesToShape.py

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