src/Tools/blocFissure/gmu/quadranglesToShape.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 """Remarque : cette focntion n'est jamais appelée ????"""
  21
  22 import logging
  23 import math
  24 import numpy as np
  25
  26 import GEOM
  27
  28 from .geomsmesh import geompy
  29 from .geomsmesh import geomPublish
  30 from .geomsmesh import geomPublishInFather
  31
  32 from . import initLog
  33
  34 def mydot(a):
  35   """produit scalaire"""
  36   return np.dot(a,a)
  37
  38
  39 def quadranglesToShape(meshQuad, shapeFissureParams, centreFondFiss):
  40   """
  41   groupe de quadrangles de face transformée en faces géométriques par filling
  42   on part de quadrangles définissant une zone a 4 cotés (convexe), et on reconstitue n lignes de p points.
  43   Ces n lignes de p points sont transformées en n courbes géométriques,
  44   à partir desquelles on reconstitue une surface géométrique.
  45   Il peut y avoir plusieurs faces géométriques reconstituées, si on fournit des groupes de quadrangles non connexes.
  46   On détecte les angles vifs, pour conserver des arêtes vives délimitant des faces connexes.
  47   @param meshQuad : maillages constitué de quadrangles constituant une ou plusieurs zones convexes
  48   @return (fillings, noeuds_Bords) : liste de geomObject, listes des bords (bord = liste ordonnée de noeuds (geomObject))
  49   """
  50   logging.info("start")
  51
  52   isVecteurDefaut = False
  53   if 'vecteurDefaut' in shapeFissureParams:
  54     isVecteurDefaut = True
  55     vecteurDefaut = shapeFissureParams['vecteurDefaut']
  56
  57   fillings = list()       # les faces reconstituées, découpées selon les arêtes vives
  58   noeuds_bords = list()   #
  59   bords_Partages = list() # contient a la fin les courbes correspondant aux arêtes vives
  60   fillconts = list()      # les faces reconstituées, sans découpage selon les arêtes vives
  61   idFilToCont = list()    # index face découpée vers face sans découpe
  62   iface = 0           # index face découpée
  63   icont = 0           # index face continue
  64
  65   allNodeIds = meshQuad.GetNodesId()
  66   while 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         elem = elems[0]
  73         break
  74     idStart = idNode # le noeud de coin
  75     elemStart = elem # l'élément quadrangle au coin
  76     xyz = meshQuad.GetNodeXYZ(idStart)
  77     logging.debug("idStart %s, coords %s", idStart, str(xyz))
  78
  79     nodelines = list() # on va constituer une liste de lignes de points
  80     nextLine = True
  81     ligneFinale = False
  82     while nextLine:
  83       logging.debug("--- une ligne")
  84       idNode = idStart
  85       elem = elemStart
  86       if ligneFinale:
  87         agauche = False      # sens de parcours des 4 noeuds d'un quadrangle
  88         nextLine = False
  89       else:
  90         agauche = True
  91       ligneIncomplete = True # on commence une ligne de points
  92       debutLigne = True
  93       nodeline = list()
  94       elemline = list()
  95       while ligneIncomplete: # compléter la ligne de points
  96         nodeline.append(idNode)
  97         allNodeIds.remove(idNode)
  98         elemline.append(elem)
  99         nodes = meshQuad.GetElemNodes(elem)
 100         i = nodes.index(idNode) # repérer l'index du noeud courant (i) dans l'élément quadrangle (0 a 3)
 101         if agauche:             # déterminer le noeud suivant (j) et celui opposé (k) dans le quadrangle
 102           if i < 3:
 103             j = i+1
 104           else:
 105             j = 0
 106           if j < 3:
 107             k = j+1
 108           else:
 109             k = 0
 110         else:
 111           if i > 0:
 112             j = i -1
 113           else:
 114             j = 3
 115           if j > 0:
 116             k = j -1
 117           else:
 118             k = 3
 119         isuiv = nodes[j]   #noeud suivant
 120         iapres = nodes[k]  #noeud opposé
 121         if debutLigne:
 122           debutLigne = False
 123           # précédent a trouver, dernière ligne : précédent au lieu de suivant
 124           if agauche:
 125             if i > 0:
 126               iprec = nodes[i -1]
 127             else:
 128               iprec = nodes[3]
 129             idStart = iprec
 130             elems3 = meshQuad.GetNodeInverseElements(iprec)
 131             if len(elems3) == 1: # autre coin
 132               ligneFinale = True
 133             else:
 134               for elem3 in elems3:
 135                 if elem3 != elem:
 136                   elemStart = elem3
 137                   break
 138         #print nodes, idNode, isuiv, iapres
 139         elems1 = meshQuad.GetNodeInverseElements(isuiv)
 140         elems2 = meshQuad.GetNodeInverseElements(iapres)
 141         ligneIncomplete = False
 142         for elem2 in elems2:
 143           if elems1.count(elem2) and elem2 != elem:
 144             ligneIncomplete = True
 145             idNode = isuiv
 146             elem = elem2
 147             break
 148         if not  ligneIncomplete:
 149           nodeline.append(isuiv)
 150           allNodeIds.remove(isuiv)
 151       logging.debug("nodeline %s", nodeline)
 152       logging.debug("elemline %s", elemline)
 153       nodelines.append(nodeline)
 154
 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 = list()
 196     bordsPartages = list()
 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 = list()
 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 = list()
 236       noeudsBords = list()
 237       for i in range(4):
 238         noeudsBords.append([])
 239       k = 0
 240       for i in range(nbLignes):
 241         nodeList = list()
 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       # --- loop on mats
 330     # --- reconstruction des faces continues à partir des listes de noeuds
 331     #     les courbes doivent suivre la courbure pour éviter les oscillations
 332     if icont == iface - 1: # pas de découpe, on garde la même face
 333       fillcont = fillings[-1]
 334     else:
 335       nbLignes = len(curvconts[0])
 336       curves = list()
 337       for i in range(nbLignes):
 338         nodes = [curvconts[j][i] for j in range(len(curvconts))]
 339         curve = geompy.MakeInterpol(nodes, False, False)
 340         curves.append(curve)
 341       fillcont = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
 342     geomPublish(initLog.debug,  fillcont, "filcont%d"%icont )
 343     fillconts.append(fillcont)
 344     icont = icont+1
 345     # --- loop while there are remaining nodes
 346
 347   return fillings, noeuds_bords, bords_Partages, fillconts, idFilToCont