1 # -*- coding: utf-8 -*-
4 from .geomsmesh import geompy
5 from .geomsmesh import geomPublish
6 from .geomsmesh import geomPublishInFather
15 # -----------------------------------------------------------------------------
16 # --- groupe de quadrangles de face transforme en face geometrique par filling
18 def quadranglesToShape(meshQuad, shapeFissureParams, centreFondFiss):
20 groupe de quadrangles de face transformee en faces geometriques par filling
21 on part de quadrangles definissant une zone a 4 cotes (convexe), et on reconstitue n lignes de p points.
22 Ces n lignes de p points sont transformees en n courbes geometriques,
23 a partir desquelles on reconstitue une surface geometrique.
24 Il peut y avoir plusieurs faces geometriques reconstituees, si on fournit des groupes de quadrangles non connexes.
25 On detecte les angles vifs, pour conserver des aretes vives delimitant des faces connexes.
26 @param meshQuad : maillages constitue de quadrangles constituant une ou plusieurs zones convexes
27 @return (fillings, noeuds_Bords) : liste de geomObject, listes des bords (bord = liste ordonnee de noeuds (geomObject))
31 isVecteurDefaut = False
32 if 'vecteurDefaut' in shapeFissureParams:
33 isVecteurDefaut = True
34 vecteurDefaut = shapeFissureParams['vecteurDefaut']
36 fillings = [] # les faces reconstituees, decoupees selon les aretes vives
38 bords_Partages = [] # contient a la fin les courbes correspondant aux aretes vives
39 fillconts = [] # les faces reconstituees, sans decoupage selon les aretes vives
40 idFilToCont = [] # index face decoupee vers face sans decoupe
41 iface = 0 # index face decoupee
42 icont = 0 # index face continue
44 allNodeIds = meshQuad.GetNodesId()
45 while len(allNodeIds):
47 for idNode in nodeIds: # rechercher un coin
48 elems = meshQuad.GetNodeInverseElements(idNode)
50 # un coin: un noeud, un element quadrangle
53 idStart = idNode # le noeud de coin
54 elemStart = elem # l'element quadrangle au coin
55 xyz = meshQuad.GetNodeXYZ(idStart)
56 logging.debug("idStart %s, coords %s", idStart, str(xyz))
58 nodelines =[] # on va constituer une liste de lignes de points
62 logging.debug("--- une ligne")
66 agauche = False # sens de parcours des 4 noeuds d'un quadrangle
70 ligneIncomplete = True # on commence une ligne de points
74 while ligneIncomplete: # completer la ligne de points
75 nodeline.append(idNode)
76 allNodeIds.remove(idNode)
78 nodes = meshQuad.GetElemNodes(elem)
79 i = nodes.index(idNode) # reperer l'index du noeud courant (i) dans l'element quadrangle (0 a 3)
80 if agauche: # determiner le noeud suivant (j) et celui oppose (k) dans le quadrangle
98 isuiv = nodes[j] #noeud suivant
99 iapres = nodes[k] #noeud oppose
102 # precedent a trouver, derniere ligne : precedent au lieu de suivant
109 elems3 = meshQuad.GetNodeInverseElements(iprec)
110 if len(elems3) == 1: # autre coin
117 #print nodes, idNode, isuiv, iapres
118 elems1 = meshQuad.GetNodeInverseElements(isuiv)
119 elems2 = meshQuad.GetNodeInverseElements(iapres)
120 ligneIncomplete = False
122 if elems1.count(elem2) and elem2 != elem:
123 ligneIncomplete = True
127 if not ligneIncomplete:
128 nodeline.append(isuiv)
129 allNodeIds.remove(isuiv)
130 logging.debug("nodeline %s", nodeline)
131 logging.debug("elemline %s", elemline)
132 nodelines.append(nodeline)
134 # on a constitue une liste de lignes de points connexes
135 logging.debug("dimensions [%s, %s]", len(nodelines), len(nodeline))
137 # stockage des coordonnees dans un tableau numpy
138 mat = np.zeros((len(nodelines), len(nodeline), 3))
139 for i, ligne in enumerate(nodelines):
140 for j, nodeId in enumerate(ligne):
141 mat[i,j] = meshQuad.GetNodeXYZ(nodeId)
142 logging.debug("matrice de coordonnees: \n%s",mat)
143 logging.debug("dimensions %s", mat.shape)
145 # recherche d'angles superieurs a un seuil sur une ligne : angle entre deux vecteurs successifs
146 cosmin = math.cos(math.pi/4.) # TODO: angle reference en parametre
147 vecx = mat[:, 1:, :] - mat[:, :-1, :] # vecteurs selon direction "x"
148 vx0 = vecx[:, :-1, :] # vecteurs amont
149 vx1 = vecx[:, 1:, :] # vecteurs aval
150 e = np.einsum('ijk,ijk->ij', vx0, vx1) # produit scalaire des vecteurs
151 f = np.apply_along_axis(mydot, 2, vx0) # normes carrees vecteurs amont
152 g = np.apply_along_axis(mydot, 2, vx1) # normes carrees vecteurs aval
153 h = e/(np.sqrt(f*g)) # cosinus
154 ruptureX = h < cosmin # True si angle > reference
155 logging.debug("matrice de rupture X: \n%s",ruptureX)
156 rupX = [x for x in range(len(nodeline)-2) if np.prod(ruptureX[:,x])]
157 logging.debug("colonnes de rupture: %s",rupX)
158 # recherche d'angles superieurs a un seuil sur une colonne : angle entre deux vecteurs successifs
159 vecy = mat[ 1:, :, :] - mat[:-1, :, :] # vecteurs selon direction "y"
160 vy0 = vecy[:-1, :, :] # vecteurs amont
161 vy1 = vecy[ 1:, :, :] # vecteurs aval
162 e = np.einsum('ijk,ijk->ij', vy0, vy1) # produit scalaire des vecteurs
163 f = np.apply_along_axis(mydot, 2, vy0) # normes carrees vecteurs amont
164 g = np.apply_along_axis(mydot, 2, vy1) # normes carrees vecteurs aval
165 h = e/(np.sqrt(f*g)) # cosinus
166 ruptureY = h < cosmin # True si angle > reference
167 logging.debug("matrice de rupture Y: \n%s",ruptureY)
168 rupY = [x for x in range(len(nodelines)-2) if np.prod(ruptureY[x, :])]
169 logging.debug("lignes de rupture: %s",rupY)
170 if (len(rupX)*len(rupY)) > 0:
171 logging.critical("""Cas non traite: presence d'angles vifs dans 2 directions,
172 lors de la reconstitution des faces geometriques dans la zone remaillee""")
177 rupX.append(mat.shape[1]-1)
178 for i, index in enumerate(rupX):
183 mats.append(mat[:, imin:imax, :])
184 if imax == mat.shape[1] + 1:
188 bordsPartages.append([imin,ifin]) # les indices differents de 0 correspondent a des bords partages
190 rupY.append(mat.shape[0]-1)
191 for i, index in enumerate(rupY):
196 mats.append(mat[imin:imax, :, :])
197 if imax == mat.shape[0] + 1:
201 bordsPartages.append([imin,ifin]) # les indices differents de 0 correspondent a des bords partages
204 bordsPartages.append([0,0]) # les indices differents de 0 correspondent a des bords partages
207 for nmat, amat in enumerate(mats):
208 logging.debug("dimensions matrice %s: %s", nmat, amat.shape)
209 nbLignes = amat.shape[1] # pas de rupture, ou rupture selon des colonnes: on transpose
210 nbCols = amat.shape[0]
211 if len(rupY) > 0 : # rupture selon des lignes: pas de transposition
212 nbLignes = amat.shape[0]
213 nbCols = amat.shape[1]
217 noeudsBords.append([])
219 for i in range(nbLignes):
221 for j in range(nbCols):
222 #logging.debug("point[%s,%s] = (%s, %s, %s)",i,j,amat[i,j,0], amat[i,j,1], amat[i,j,2])
223 if len(rupY) > 0 : # pas de transposition
224 node = geompy.MakeVertex(amat[i,j,0], amat[i,j,1], amat[i,j,2])
225 else: # transposition
226 node = geompy.MakeVertex(amat[j,i,0], amat[j,i,1], amat[j,i,2])
227 nodeList.append(node)
229 noeudsBords[0].append(node)
231 #geomPublish(initLog.debug, node, name )
232 if i == (nbLignes -1):
233 noeudsBords[2].append(node)
235 #geomPublish(initLog.debug, node, name )
237 noeudsBords[1].append(node)
239 #geomPublish(initLog.debug, node, name )
241 noeudsBords[3].append(node)
243 #geomPublish(initLog.debug, node, name )
245 curve = geompy.MakeInterpol(nodeList, False, False)
247 #geomPublish(initLog.debug, curve, name )
248 if len(curvconts) == 0 or len(curves) > 0: # eliminer les doublons de la surface sans decoupe
249 curvconts.append(nodeList)
251 if bordsPartages[nmat][0] :
252 bordsPartages[nmat][0] = curves[0] # la premiere ligne est un bord partage
254 bordsPartages[nmat][0] = None
255 if bordsPartages[nmat][1] :
256 bordsPartages[nmat][1] = curves[-1] # la derniere ligne est un bord partage
258 bordsPartages[nmat][1] = None
259 filling = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
260 # --- test orientation filling
261 vertex = geompy.MakeVertexOnSurface(filling, 0.5, 0.5)
262 normal = geompy.GetNormal(filling, vertex)
264 if centreFondFiss is not None:
265 logging.debug("orientation filling a l'aide du centre de fond de fissure")
266 vecteurDefaut = geompy.MakeVector(centreFondFiss, vertex)
268 if not isVecteurDefaut:
272 pointExplicite = False
273 if 'pointIn_x' in shapeFissureParams:
274 pointExplicite = True
275 pointIn_x = shapeFissureParams['pointIn_x']
276 if 'pointIn_y' in shapeFissureParams:
277 pointExplicite = True
278 pointIn_y = shapeFissureParams['pointIn_y']
279 if 'pointIn_z' in shapeFissureParams:
280 pointExplicite = True
281 pointIn_z = shapeFissureParams['pointIn_z']
283 cdg = geompy.MakeVertex(pointIn_x, pointIn_y, pointIn_z)
284 logging.debug("orientation filling par point interieur %s", (pointIn_x, pointIn_y, pointIn_z))
285 vecteurDefaut = geompy.MakeVector(cdg, vertex)
287 if 'convexe' in shapeFissureParams:
288 isConvexe = shapeFissureParams['convexe']
289 logging.debug("orientation filling par indication de convexite %s", isConvexe)
290 cdg = geompy.MakeCDG(filling)
292 vecteurDefaut = geompy.MakeVector(cdg, vertex)
294 vecteurDefaut = geompy.MakeVector(vertex, cdg)
296 if vecteurDefaut is not None:
297 geomPublish(initLog.debug, normal, "normFillOrig%d"%iface)
298 geomPublish(initLog.debug, vecteurDefaut, "fromInterieur%d"%iface)
299 if geompy.GetAngleRadians(vecteurDefaut, normal) > math.pi/2.0:
300 filling = geompy.ChangeOrientation(filling)
301 geomPublish(initLog.debug, filling, "filling%d"%iface )
302 #geompy.ExportBREP(filling, "filling.brep")
304 fillings.append(filling)
305 noeuds_bords.append(noeudsBords)
306 idFilToCont.append(icont)
307 bords_Partages += bordsPartages
308 pass # --- loop on mats
309 # --- reconstruction des faces continues a partir des listes de noeuds
310 # les courbes doivent suivre la courbure pour eviter les oscillations
311 if icont == iface - 1: # pas de decoupe, on garde la meme face
312 fillcont = fillings[-1]
314 nbLignes = len(curvconts[0])
316 for i in range(nbLignes):
317 nodes = [curvconts[j][i] for j in range(len(curvconts))]
318 curve = geompy.MakeInterpol(nodes, False, False)
320 fillcont = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
321 geomPublish(initLog.debug, fillcont, "filcont%d"%icont )
322 fillconts.append(fillcont)
324 pass # --- loop while there are remaining nodes
326 return fillings, noeuds_bords, bords_Partages, fillconts, idFilToCont