1 # -*- coding: utf-8 -*-
4 from geomsmesh import geompy
12 # -----------------------------------------------------------------------------
13 # --- groupe de quadrangles de face transformé en face géométrique par filling
15 def quadranglesToShapeNoCorner(meshQuad, shapeFissureParams, centreFondFiss):
17 groupe de quadrangles de face transformée en faces géométriques par filling
18 on part de quadrangles définissant une zone a 4 cotés (convexe), et on reconstitue n lignes de p points.
19 Ces n lignes de p points sont transformées en n courbes géométriques,
20 à partir desquelles on reconstitue une surface géométrique.
21 Il peut y avoir plusieurs faces géométriques reconstituées, si on fournit des groupes de quadrangles non connexes.
22 On détecte les angles vifs, pour conserver des arêtes vives délimitant des faces connexes.
23 @param meshQuad : maillages constitué de quadrangles constituant une ou plusieurs zones convexes
24 @return (fillings, noeuds_Bords) : liste de geomObject, listes des bords (bord = liste ordonnée de noeuds (geomObject))
28 isVecteurDefaut = False
29 if shapeFissureParams.has_key('vecteurDefaut'):
30 isVecteurDefaut = True
31 vecteurDefaut = shapeFissureParams['vecteurDefaut']
33 fillings = [] # les faces reconstituées, découpées selon les arêtes vives
35 bords_Partages = [] # contient a la fin les courbes correspondant aux arêtes vives
36 fillconts = [] # les faces reconstituées, sans découpage selon les arêtes vives
37 idFilToCont = [] # index face découpée vers face sans découpe
38 iface = 0 # index face découpée
39 icont = 0 # index face continue
41 allNodeIds = meshQuad.GetNodesId()
42 while len(allNodeIds):
44 for idNode in nodeIds: # rechercher un coin
45 elems = meshQuad.GetNodeInverseElements(idNode)
47 # un coin: un noeud, un element quadrangle
50 idStart = idNode # le noeud de coin
51 elemStart = elem # l'élément quadrangle au coin
52 xyz = meshQuad.GetNodeXYZ(idStart)
53 logging.debug("idStart %s, coords %s", idStart, str(xyz))
55 nodelines =[] # on va constituer une liste de lignes de points
59 logging.debug("--- une ligne")
63 agauche = False # sens de parcours des 4 noeuds d'un quadrangle
67 ligneIncomplete = True # on commence une ligne de points
71 while ligneIncomplete: # compléter la ligne de points
72 nodeline.append(idNode)
73 allNodeIds.remove(idNode)
75 nodes = meshQuad.GetElemNodes(elem)
76 i = nodes.index(idNode) # repérer l'index du noeud courant (i) dans l'élément quadrangle (0 a 3)
77 if agauche: # déterminer le noeud suivant (j) et celui opposé (k) dans le quadrangle
95 isuiv = nodes[j] #noeud suivant
96 iapres = nodes[k] #noeud opposé
99 # précédent a trouver, dernière ligne : précédent au lieu de suivant
106 elems3 = meshQuad.GetNodeInverseElements(iprec)
107 if len(elems3) == 1: # autre coin
114 #print nodes, idNode, isuiv, iapres
115 elems1 = meshQuad.GetNodeInverseElements(isuiv)
116 elems2 = meshQuad.GetNodeInverseElements(iapres)
117 ligneIncomplete = False
119 if elems1.count(elem2) and elem2 != elem:
120 ligneIncomplete = True
124 if not ligneIncomplete:
125 nodeline.append(isuiv)
126 allNodeIds.remove(isuiv)
127 logging.debug("nodeline %s", nodeline)
128 logging.debug("elemline %s", elemline)
129 nodelines.append(nodeline)
130 logging.debug("nodelines = %s", nodelines)
131 longueur = [len(val) for val in nodelines]
132 logging.debug("longueur = %s", longueur)
133 # on a constitué une liste de lignes de points connexes
134 logging.debug("dimensions [%s, %s]", len(nodelines), len(nodeline))
136 # stockage des coordonnées dans un tableau numpy
137 mat = np.zeros((len(nodelines), len(nodeline), 3))
138 for i, ligne in enumerate(nodelines):
139 for j, nodeId in enumerate(ligne):
140 mat[i,j] = meshQuad.GetNodeXYZ(nodeId)
141 logging.debug("matrice de coordonnées: \n%s",mat)
142 logging.debug("dimensions %s", mat.shape)
144 # recherche d'angles supérieurs a un seuil sur une ligne : angle entre deux vecteurs successifs
145 cosmin = math.cos(math.pi/4.) # TODO: angle reference en paramètre
146 vecx = mat[:, 1:, :] - mat[:, :-1, :] # vecteurs selon direction "x"
147 vx0 = vecx[:, :-1, :] # vecteurs amont
148 vx1 = vecx[:, 1:, :] # vecteurs aval
149 e = np.einsum('ijk,ijk->ij', vx0, vx1) # produit scalaire des vecteurs
150 f = np.apply_along_axis(mydot, 2, vx0) # normes carrées vecteurs amont
151 g = np.apply_along_axis(mydot, 2, vx1) # normes carrées vecteurs aval
152 h = e/(np.sqrt(f*g)) # cosinus
153 ruptureX = h < cosmin # True si angle > reference
154 logging.debug("matrice de rupture X: \n%s",ruptureX)
155 rupX = filter(lambda x: np.prod(ruptureX[:,x]), range(len(nodeline)-2))
156 logging.debug("colonnes de rupture: %s",rupX)
157 # recherche d'angles supérieurs a un seuil sur une colonne : angle entre deux vecteurs successifs
158 vecy = mat[ 1:, :, :] - mat[:-1, :, :] # vecteurs selon direction "y"
159 vy0 = vecy[:-1, :, :] # vecteurs amont
160 vy1 = vecy[ 1:, :, :] # vecteurs aval
161 e = np.einsum('ijk,ijk->ij', vy0, vy1) # produit scalaire des vecteurs
162 f = np.apply_along_axis(mydot, 2, vy0) # normes carrées vecteurs amont
163 g = np.apply_along_axis(mydot, 2, vy1) # normes carrées vecteurs aval
164 h = e/(np.sqrt(f*g)) # cosinus
165 ruptureY = h < cosmin # True si angle > reference
166 logging.debug("matrice de rupture Y: \n%s",ruptureY)
167 rupY = filter(lambda x: np.prod(ruptureY[x, :]), range(len(nodelines)-2))
168 logging.debug("lignes de rupture: %s",rupY)
169 if (len(rupX)*len(rupY)) > 0:
170 logging.critical("""Cas non traité: présence d'angles vifs dans 2 directions,
171 lors de la reconstitution des faces géométriques dans la zone remaillée""")
176 rupX.append(mat.shape[1]-1)
177 for i, index in enumerate(rupX):
182 mats.append(mat[:, imin:imax, :])
183 if imax == mat.shape[1] + 1:
187 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
189 rupY.append(mat.shape[0]-1)
190 for i, index in enumerate(rupY):
195 mats.append(mat[imin:imax, :, :])
196 if imax == mat.shape[0] + 1:
200 bordsPartages.append([imin,ifin]) # les indices différents de 0 correspondent à des bords partagés
203 bordsPartages.append([0,0]) # les indices différents de 0 correspondent à des bords partagés
206 for nmat, amat in enumerate(mats):
207 logging.debug("dimensions matrice %s: %s", nmat, amat.shape)
208 nbLignes = amat.shape[1] # pas de rupture, ou rupture selon des colonnes: on transpose
209 nbCols = amat.shape[0]
210 if len(rupY) > 0 : # rupture selon des lignes: pas de transposition
211 nbLignes = amat.shape[0]
212 nbCols = amat.shape[1]
216 noeudsBords.append([])
218 for i in range(nbLignes):
220 for j in range(nbCols):
221 #logging.debug("point[%s,%s] = (%s, %s, %s)",i,j,amat[i,j,0], amat[i,j,1], amat[i,j,2])
222 if len(rupY) > 0 : # pas de transposition
223 node = geompy.MakeVertex(amat[i,j,0], amat[i,j,1], amat[i,j,2])
224 else: # transposition
225 node = geompy.MakeVertex(amat[j,i,0], amat[j,i,1], amat[j,i,2])
226 nodeList.append(node)
228 noeudsBords[0].append(node)
230 #geompy.addToStudy( node, name )
231 if i == (nbLignes -1):
232 noeudsBords[2].append(node)
234 #geompy.addToStudy( node, name )
236 noeudsBords[1].append(node)
238 #geompy.addToStudy( node, name )
240 noeudsBords[3].append(node)
242 #geompy.addToStudy( node, name )
244 curve = geompy.MakeInterpol(nodeList, False, False)
246 #geompy.addToStudy( curve, name )
247 if len(curvconts) == 0 or len(curves) > 0: # éliminer les doublons de la surface sans découpe
248 curvconts.append(nodeList)
250 if bordsPartages[nmat][0] :
251 bordsPartages[nmat][0] = curves[0] # la première ligne est un bord partagé
253 bordsPartages[nmat][0] = None
254 if bordsPartages[nmat][1] :
255 bordsPartages[nmat][1] = curves[-1] # la dernière ligne est un bord partagé
257 bordsPartages[nmat][1] = None
258 filling = geompy.MakeFilling(geompy.MakeCompound(curves), 2, 5, 0.0001, 0.0001, 0, GEOM.FOM_Default, True)
259 # --- test orientation filling
260 vertex = geompy.MakeVertexOnSurface(filling, 0.5, 0.5)
261 normal = geompy.GetNormal(filling, vertex)
263 if centreFondFiss is not None:
264 logging.debug("orientation filling a l'aide du centre de fond de fissure")
265 vecteurDefaut = geompy.MakeVector(centreFondFiss, vertex)
267 if not isVecteurDefaut:
271 pointExplicite = False
272 if shapeFissureParams.has_key('pointIn_x'):
273 pointExplicite = True
274 pointIn_x = shapeFissureParams['pointIn_x']
275 if shapeFissureParams.has_key('pointIn_y'):
276 pointExplicite = True
277 pointIn_y = shapeFissureParams['pointIn_y']
278 if shapeFissureParams.has_key('pointIn_z'):
279 pointExplicite = True
280 pointIn_z = shapeFissureParams['pointIn_z']
282 cdg = geompy.MakeVertex(pointIn_x, pointIn_y, pointIn_z)
283 logging.debug("orientation filling par point intérieur %s", (pointIn_x, pointIn_y, pointIn_z))
284 vecteurDefaut = geompy.MakeVector(cdg, vertex)
286 if shapeFissureParams.has_key('convexe'):
287 isConvexe = shapeFissureParams['convexe']
288 logging.debug("orientation filling par indication de convexité %s", isConvexe)
289 cdg = geompy.MakeCDG(filling)
291 vecteurDefaut = geompy.MakeVector(cdg, vertex)
293 vecteurDefaut = geompy.MakeVector(vertex, cdg)
295 if vecteurDefaut is not None:
296 geompy.addToStudy(normal, "normFillOrig%d"%iface)
297 geompy.addToStudy(vecteurDefaut, "fromInterieur%d"%iface)
298 if geompy.GetAngleRadians(vecteurDefaut, normal) > math.pi/2.0:
299 filling = geompy.ChangeOrientation(filling)
300 geompy.addToStudy( filling, "filling%d"%iface )
301 #geompy.ExportBREP(filling, "filling.brep")
303 fillings.append(filling)
304 noeuds_bords.append(noeudsBords)
305 idFilToCont.append(icont)
306 bords_Partages += bordsPartages
307 logging.debug("bords_Partages = %s", bords_Partages)
308 pass # --- loop on mats
309 # --- reconstruction des faces continues à partir des listes de noeuds
310 # les courbes doivent suivre la courbure pour éviter les oscillations
311 if icont == iface - 1: # pas de découpe, on garde la même 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 geompy.addToStudy( 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