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