1 # -*- coding: utf-8 -*-
2 # Copyright (C) 2014-2019 CEA/DEN, 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
22 from .geomsmesh import geompy
23 from .geomsmesh import geomPublish
24 from .geomsmesh import geomPublishInFather
33 # -----------------------------------------------------------------------------
34 # --- groupe de quadrangles de face transformé en face géométrique par filling
36 def quadranglesToShape(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 = [] # les faces reconstituées, découpées selon les arêtes vives
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
62 allNodeIds = meshQuad.GetNodesId()
63 while len(allNodeIds):
65 for idNode in nodeIds: # rechercher un coin
66 elems = meshQuad.GetNodeInverseElements(idNode)
68 # un coin: un noeud, un element quadrangle
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))
76 nodelines =[] # on va constituer une liste de lignes de points
80 logging.debug("--- une ligne")
84 agauche = False # sens de parcours des 4 noeuds d'un quadrangle
88 ligneIncomplete = True # on commence une ligne de points
92 while ligneIncomplete: # compléter la ligne de points
93 nodeline.append(idNode)
94 allNodeIds.remove(idNode)
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
116 isuiv = nodes[j] #noeud suivant
117 iapres = nodes[k] #noeud opposé
120 # précédent a trouver, dernière ligne : précédent au lieu de suivant
127 elems3 = meshQuad.GetNodeInverseElements(iprec)
128 if len(elems3) == 1: # autre coin
135 #print nodes, idNode, isuiv, iapres
136 elems1 = meshQuad.GetNodeInverseElements(isuiv)
137 elems2 = meshQuad.GetNodeInverseElements(iapres)
138 ligneIncomplete = False
140 if elems1.count(elem2) and elem2 != elem:
141 ligneIncomplete = True
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)
152 # on a constitué une liste de lignes de points connexes
153 logging.debug("dimensions [%s, %s]", len(nodelines), len(nodeline))
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)
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""")
195 rupX.append(mat.shape[1]-1)
196 for i, index in enumerate(rupX):
201 mats.append(mat[:, imin:imax, :])
202 if imax == mat.shape[1] + 1:
206 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
208 rupY.append(mat.shape[0]-1)
209 for i, index in enumerate(rupY):
214 mats.append(mat[imin:imax, :, :])
215 if imax == mat.shape[0] + 1:
219 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
222 bordsPartages.append([0,0]) # les indices différents de 0 correspondent à des bords partagés
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]
235 noeudsBords.append([])
237 for i in range(nbLignes):
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)
247 noeudsBords[0].append(node)
249 #geomPublish(initLog.debug, node, name )
250 if i == (nbLignes -1):
251 noeudsBords[2].append(node)
253 #geomPublish(initLog.debug, node, name )
255 noeudsBords[1].append(node)
257 #geomPublish(initLog.debug, node, name )
259 noeudsBords[3].append(node)
261 #geomPublish(initLog.debug, node, name )
263 curve = geompy.MakeInterpol(nodeList, False, False)
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)
269 if bordsPartages[nmat][0] :
270 bordsPartages[nmat][0] = curves[0] # la première ligne est un bord partagé
272 bordsPartages[nmat][0] = None
273 if bordsPartages[nmat][1] :
274 bordsPartages[nmat][1] = curves[-1] # la dernière ligne est un bord partagé
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)
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)
286 if not isVecteurDefaut:
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']
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)
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)
310 vecteurDefaut = geompy.MakeVector(cdg, vertex)
312 vecteurDefaut = geompy.MakeVector(vertex, cdg)
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")
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]
332 nbLignes = len(curvconts[0])
334 for i in range(nbLignes):
335 nodes = [curvconts[j][i] for j in range(len(curvconts))]
336 curve = geompy.MakeInterpol(nodes, False, False)
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)
342 pass # --- loop while there are remaining nodes
344 return fillings, noeuds_bords, bords_Partages, fillconts, idFilToCont