1 # -*- coding: utf-8 -*-
2 # Copyright (C) 2014-2021 EDF R&D
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.
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.
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
18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 """Groupe de quadrangles de face transformé en face géométrique par filling"""
27 from .geomsmesh import geompy
28 from .geomsmesh import geomPublish
33 """produit scalaire"""
36 def quadranglesToShapeNoCorner(meshQuad, shapeFissureParams, centreFondFiss):
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))
49 isVecteurDefaut = False
50 if 'vecteurDefaut' in shapeFissureParams:
51 isVecteurDefaut = True
52 vecteurDefaut = shapeFissureParams['vecteurDefaut']
54 fillings = list() # les faces reconstituées, découpées selon les arêtes vives
55 noeuds_bords = list() #
56 bords_Partages = list() # contient a la fin les courbes correspondant aux arêtes vives
57 fillconts = list() # les faces reconstituées, sans découpage selon les arêtes vives
58 idFilToCont = list() # index face découpée vers face sans découpe
59 iface = 0 # index face découpée
60 icont = 0 # index face continue
64 allNodeIds = meshQuad.GetNodesId()
65 while len(allNodeIds):
66 logging.debug("len(allNodeIds): %s ", len(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 idStart = idNode # le noeud de coin
73 elemStart = elems[0] # l'élément quadrangle au coin
75 xyz = meshQuad.GetNodeXYZ(idStart)
76 logging.debug("idStart %s, coords %s", idStart, str(xyz))
78 nodelines = list() # on va constituer une liste de lignes de points
82 logging.debug("--- une ligne")
86 agauche = False # sens de parcours des 4 noeuds d'un quadrangle
90 ligneIncomplete = True # on commence une ligne de points
94 while ligneIncomplete: # compléter la ligne de points
95 nodeline.append(idNode)
96 allNodeIds.remove(idNode)
98 nodes = meshQuad.GetElemNodes(elem)
99 i = nodes.index(idNode) # repérer l'index du noeud courant (i) dans l'élément quadrangle (0 a 3)
100 if agauche: # déterminer le noeud suivant (j) et celui opposé (k) dans le quadrangle
118 isuiv = nodes[j] #noeud suivant
119 iapres = nodes[k] #noeud opposé
122 # précédent a trouver, dernière ligne : précédent au lieu de suivant
129 elems3 = meshQuad.GetNodeInverseElements(iprec)
130 if len(elems3) == 1: # autre coin
137 #print nodes, idNode, isuiv, iapres
138 elems1 = meshQuad.GetNodeInverseElements(isuiv)
139 elems2 = meshQuad.GetNodeInverseElements(iapres)
140 ligneIncomplete = False
142 if elems1.count(elem2) and elem2 != elem:
143 ligneIncomplete = True
147 if not ligneIncomplete:
148 nodeline.append(isuiv)
149 allNodeIds.remove(isuiv)
150 logging.debug("nodeline %s", nodeline)
151 logging.debug("elemline %s", elemline)
152 nodelines.append(nodeline)
153 logging.debug("nodelines = %s", nodelines)
154 longueur = [len(val) for val in nodelines]
155 logging.debug("longueur = %s", longueur)
156 # on a constitué une liste de lignes de points connexes
157 logging.debug("dimensions [%s, %s]", len(nodelines), len(nodeline))
159 # stockage des coordonnées dans un tableau numpy
160 mat = np.zeros((len(nodelines), len(nodeline), 3))
161 for i, ligne in enumerate(nodelines):
162 for j, nodeId in enumerate(ligne):
163 mat[i,j] = meshQuad.GetNodeXYZ(nodeId)
164 logging.debug("matrice de coordonnées: \n%s",mat)
165 logging.debug("dimensions %s", mat.shape)
167 # recherche d'angles supérieurs a un seuil sur une ligne : angle entre deux vecteurs successifs
168 cosmin = np.cos(pisur4) # TODO: angle reference en paramètre
169 vecx = mat[:, 1:, :] - mat[:, :-1, :] # vecteurs selon direction "x"
170 vx0 = vecx[:, :-1, :] # vecteurs amont
171 vx1 = vecx[:, 1:, :] # vecteurs aval
172 e = np.einsum('ijk,ijk->ij', vx0, vx1) # produit scalaire des vecteurs
173 f = np.apply_along_axis(mydot, 2, vx0) # normes carrées vecteurs amont
174 g = np.apply_along_axis(mydot, 2, vx1) # normes carrées vecteurs aval
175 h = e/(np.sqrt(f*g)) # cosinus
176 ruptureX = h < cosmin # True si angle > reference
177 logging.debug("matrice de rupture X: \n%s",ruptureX)
178 rupX = [x for x in range(len(nodeline)-2) if np.prod(ruptureX[:,x])]
179 logging.debug("colonnes de rupture: %s",rupX)
180 # recherche d'angles supérieurs a un seuil sur une colonne : angle entre deux vecteurs successifs
181 vecy = mat[ 1:, :, :] - mat[:-1, :, :] # vecteurs selon direction "y"
182 vy0 = vecy[:-1, :, :] # vecteurs amont
183 vy1 = vecy[ 1:, :, :] # vecteurs aval
184 e = np.einsum('ijk,ijk->ij', vy0, vy1) # produit scalaire des vecteurs
185 f = np.apply_along_axis(mydot, 2, vy0) # normes carrées vecteurs amont
186 g = np.apply_along_axis(mydot, 2, vy1) # normes carrées vecteurs aval
187 h = e/(np.sqrt(f*g)) # cosinus
188 ruptureY = h < cosmin # True si angle > reference
189 logging.debug("matrice de rupture Y: \n%s",ruptureY)
190 rupY = [x for x in range(len(nodelines)-2) if np.prod(ruptureY[x, :])]
191 logging.debug("lignes de rupture: %s",rupY)
192 if (len(rupX)*len(rupY)) > 0:
193 logging.critical("""Cas non traité: présence d'angles vifs dans 2 directions,
194 lors de la reconstitution des faces géométriques dans la zone remaillée""")
197 bordsPartages = list()
199 rupX.append(mat.shape[1]-1)
200 for i, index in enumerate(rupX):
205 mats.append(mat[:, imin:imax, :])
206 if imax == mat.shape[1] + 1:
210 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
212 rupY.append(mat.shape[0]-1)
213 for i, index in enumerate(rupY):
218 mats.append(mat[imin:imax, :, :])
219 if imax == mat.shape[0] + 1:
223 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
226 bordsPartages.append([0,0]) # les indices différents de 0 correspondent à des bords partagés
229 for nmat, amat in enumerate(mats):
230 logging.debug("dimensions matrice %s: %s", nmat, amat.shape)
231 nbLignes = amat.shape[1] # pas de rupture, ou rupture selon des colonnes: on transpose
232 nbCols = amat.shape[0]
233 if len(rupY) > 0 : # rupture selon des lignes: pas de transposition
234 nbLignes = amat.shape[0]
235 nbCols = amat.shape[1]
239 noeudsBords.append([])
241 for i in range(nbLignes):
243 for j in range(nbCols):
244 #logging.debug("point[%s,%s] = (%s, %s, %s)",i,j,amat[i,j,0], amat[i,j,1], amat[i,j,2])
245 if len(rupY) > 0 : # pas de transposition
246 node = geompy.MakeVertex(amat[i,j,0], amat[i,j,1], amat[i,j,2])
247 else: # transposition
248 node = geompy.MakeVertex(amat[j,i,0], amat[j,i,1], amat[j,i,2])
249 nodeList.append(node)
251 noeudsBords[0].append(node)
253 #geomPublish(initLog.debug, node, name )
254 if i == (nbLignes -1):
255 noeudsBords[2].append(node)
257 #geomPublish(initLog.debug, node, name )
259 noeudsBords[1].append(node)
261 #geomPublish(initLog.debug, node, name )
263 noeudsBords[3].append(node)
265 #geomPublish(initLog.debug, node, name )
267 curve = geompy.MakeInterpol(nodeList, False, False)
269 #geomPublish(initLog.debug, curve, name )
270 if len(curvconts) == 0 or len(curves) > 0: # éliminer les doublons de la surface sans découpe
271 curvconts.append(nodeList)
273 if bordsPartages[nmat][0] :
274 bordsPartages[nmat][0] = curves[0] # la première ligne est un bord partagé
276 bordsPartages[nmat][0] = None
277 if bordsPartages[nmat][1] :
278 bordsPartages[nmat][1] = curves[-1] # la dernière ligne est un bord partagé
280 bordsPartages[nmat][1] = None
281 filling = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
282 # --- test orientation filling
283 vertex = geompy.MakeVertexOnSurface(filling, 0.5, 0.5)
284 normal = geompy.GetNormal(filling, vertex)
286 if centreFondFiss is not None:
287 logging.debug("orientation filling a l'aide du centre de fond de fissure")
288 vecteurDefaut = geompy.MakeVector(centreFondFiss, vertex)
290 if not isVecteurDefaut:
294 pointExplicite = False
295 if 'pointIn_x' in shapeFissureParams:
296 pointExplicite = True
297 pointIn_x = shapeFissureParams['pointIn_x']
298 if 'pointIn_y' in shapeFissureParams:
299 pointExplicite = True
300 pointIn_y = shapeFissureParams['pointIn_y']
301 if 'pointIn_z' in shapeFissureParams:
302 pointExplicite = True
303 pointIn_z = shapeFissureParams['pointIn_z']
305 cdg = geompy.MakeVertex(pointIn_x, pointIn_y, pointIn_z)
306 logging.debug("orientation filling par point intérieur %s", (pointIn_x, pointIn_y, pointIn_z))
307 vecteurDefaut = geompy.MakeVector(cdg, vertex)
309 if 'convexe' in shapeFissureParams:
310 isConvexe = shapeFissureParams['convexe']
311 logging.debug("orientation filling par indication de convexité %s", isConvexe)
312 cdg = geompy.MakeCDG(filling)
314 vecteurDefaut = geompy.MakeVector(cdg, vertex)
316 vecteurDefaut = geompy.MakeVector(vertex, cdg)
318 if vecteurDefaut is not None:
319 geomPublish(initLog.debug, normal, "normFillOrig%d"%iface)
320 geomPublish(initLog.debug, vecteurDefaut, "fromInterieur%d"%iface)
321 if ( geompy.GetAngleRadians(vecteurDefaut, normal) > pisur2 ):
322 filling = geompy.ChangeOrientation(filling)
323 geomPublish(initLog.debug, filling, "filling%d"%iface )
324 #geompy.ExportBREP(filling, "filling.brep")
326 fillings.append(filling)
327 noeuds_bords.append(noeudsBords)
328 idFilToCont.append(icont)
329 bords_Partages += bordsPartages
330 logging.debug("bords_Partages = %s", bords_Partages)
332 # --- reconstruction des faces continues à partir des listes de noeuds
333 # les courbes doivent suivre la courbure pour éviter les oscillations
334 if icont == iface - 1: # pas de découpe, on garde la même face
335 fillcont = fillings[-1]
337 nbLignes = len(curvconts[0])
339 for i in range(nbLignes):
340 nodes = [curvconts[j][i] for j in range(len(curvconts))]
341 curve = geompy.MakeInterpol(nodes, False, False)
343 fillcont = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
344 geomPublish(initLog.debug, fillcont, "filcont%d"%icont )
345 fillconts.append(fillcont)
347 # --- loop while there are remaining nodes
349 return fillings, noeuds_bords, bords_Partages, fillconts, idFilToCont