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