6 # methode de filling: par les generatrices ou par des ecailles de tortue
9 longueur_tube = 6363.655
11 # Facteur d'amplification pour voir si le filling passe bien par les generatrices
12 amplification_factor = 1
14 # teste si le maillage gibi passe par tous les points
16 maillage_gibi = '/home/PROJETS/OUVERT/edf_anode/from_edf/2008_11_28/Tube/tube.mail.med'
21 type_maillage = "regle" # "regle" ou "libre"
23 if type_maillage == "libre":
24 type_algo = "NETGEN_2D" # "BLSURF" ou "NETGEN_2D"
26 if methode == "generatrices":
27 # nombre de segments dans la longueur des generatrices dans la zone de sous-epaisseur
28 nb_seg_generatrices = 150 # methode generatrices
30 # nombre de segments dans la longueur des generatrices dans la zone de sous-epaisseur
31 nb_seg_generatrices = 5 # methode tortue
32 # distance entre 2 abscisses de points de mesure au dessous de laquelle on discretise avec nb_seg_petites_distances
33 # au lieu de nb_seg_generatrices
35 nb_seg_petites_distances = 3
37 # nombre de segments dans l'epaisseur du tube
40 # nombre de segments dans l'arc du tube entre 2 generatrices
41 if methode == "generatrices":
42 nb_seg_arc = 20 # methode generatrices: partition en 2
44 nb_seg_arc = 5 # methode tortue: partition en 8
46 # nombre de segments dans la longueur de transition
49 # nombre de segments dans la longueur d'amortissement
50 nb_seg_amortissement = 11
58 import os, csv, math, pdb, time
61 # ===================================================
63 # supprime recursivement un element dans une liste
65 # l=["a", " ", "b", " ", "c"]
66 # l = recursive_remove(l, " ")
69 def recursive_remove(l, txt):
79 ## lit les valeurs des epaisseurs sur les generatrices a partir d'un fichier csv
80 # @return d_generatrices : dictionnaire dont la cle est le nom de la generatrice et la valeur est la liste des epaisseur sur cette generatrice
81 # @return l_abscisses: liste des abscisses des points de mesures
82 def read_generatrice(filename):
83 file = open(filename, "rb")
85 reader = csv.reader(file, delimiter=';', quoting=csv.QUOTE_NONE)
87 # Dictionnaire pour stocker les mesures sur une generatrice donnee
90 # Liste des noms des generatrices
91 l_noms_generatrices = []
93 # Liste des abscisses des points de mesures
96 for i, row in enumerate(reader):
97 # On supprime les cases vides
98 row = recursive_remove(row, "")
99 # Sur la premiere ligne, on releve les noms des generatrices
102 if nom not in ["Abscisse", "Longueur"]:
103 # on initialise chaque generatrice avec une liste vide
104 d_generatrices[nom] = []
105 l_noms_generatrices.append(nom)
106 # nombre de generatrices trouvees:
107 nb_generatrices = len(d_generatrices)
109 # sur les lignes suivantes, on releve les mesures des epaisseurs
110 for j, nom_generatrice in enumerate(l_noms_generatrices):
111 # la liste des epaisseurs commence a 1
112 epaisseur_str_fr = row[j+1]
113 # on convertit le format decimal francais par le format anglais
114 epaisseur_str = epaisseur_str_fr.replace(",", ".")
115 epaisseur = float(epaisseur_str)
116 d_generatrices[nom_generatrice].append(epaisseur)
117 # on ajoute la valeur de l'abscisse
118 abscisse_str = row[nb_generatrices+1]
119 abscisse = float(abscisse_str)
120 l_abscisses.append(abscisse)
124 return d_generatrices, l_noms_generatrices, l_abscisses
126 ## lit les valeurs des angles definissant les generatrices a partir d'un fichier csv
127 # @return l_angles : liste des angles definissant les generatrices
128 def read_angles(filename):
129 file = open(filename, "rb")
131 reader = csv.reader(file, delimiter=';', quoting=csv.QUOTE_NONE)
133 # Liste des angles des generatrices
137 # On supprime les cases vides
138 row = recursive_remove(row, "")
140 # si la ligne comporte 3 valeurs, on peut lire les angles
143 angle = float(angle_str)
144 l_angles.append(angle)
147 ## Cree une face a partir d'un nuage de points
148 # @param l_arcs_points : liste de points sur un contour
149 # @param closed_wire: flag pour savoir si le contour est deja ferme ou s'il faut le fermer
150 # @return face : une face passant par tous les points
151 # @warning: TODO: completer l'algo pour qu'il fonctionne si le nombre de generatrices est impair!!!
152 def MakeShellFromPoints(l_arcs_points, closed_wire = False):
153 # on cree les arcs pour chaque quart de cercle
154 nb_generatrices = len(l_arcs_points[0])
155 nb_arcs_filling = nb_generatrices/2
156 l_arcs_filling = [[] for i in range(nb_arcs_filling)]
158 if nb_generatrices%2 != 1:
159 raise "L'algo ne fonctionne pour l'instant qu'avec un nombre de generatrices impair"
161 if nb_generatrices%2 != 0:
162 raise "L'algo ne fonctionne pour l'instant qu'avec un nombre de generatrices pair"
163 # Creation des arcs a partir des points
164 for arc_points in l_arcs_points:
166 # Pour cloturer le contour
167 arc_points.append(arc_points[0])
168 for i in range(nb_arcs_filling):
169 # 3 points a la meme abscisse sur 3 generatrices consecutives
170 arc = MakeArc(arc_points[2*i], arc_points[1+2*i], arc_points[2+2*i])
171 l_arcs_filling[i].append(arc)
174 # on fait un filling pour tous les 1ers arcs, tous les 2emes arcs, ... jusqu'au 4e arc
176 for i, arc_filling in enumerate(l_arcs_filling):
177 # On fait un filling entre 2 quarts de cercles
178 for j in range(0, len(arc_filling) - 1):
179 l_quart_arc = [arc_filling[j], arc_filling[j+1]]
180 compound_quart_arcs = MakeCompound(l_quart_arc)
181 quart_face_filling = MakeFilling(compound_quart_arcs, 0, 10, 1e-05, 1e-05, 0)
182 #addToStudy(quart_face_filling, "quart_face_filling_%i"%(j+1))
183 l_faces_filling.append(quart_face_filling)
184 face = MakeShell(l_faces_filling)
187 ## Cree une face a partir d'un nuage de points par un filling sur les generatrices
188 # @param l_arcs_points : liste de points sur un contour
189 # @param closed_wire: flag pour savoir si le contour est deja ferme ou s'il faut le fermer
190 # @return face : une face passant par tous les points
191 def MakeFillingFromPoints(l_arcs_points, closed_wire = False):
192 nb_generatrices = len(l_arcs_points[0])
193 l_points_generatrices = [[] for i in range(nb_generatrices)]
195 # Creation des generatrices a partir des points
196 for arc_points in l_arcs_points:
197 for i, point in enumerate(arc_points):
198 # on ajoute le point dans la generatrice correspondante
199 l_points_generatrices[i].append(point)
200 for points_generatrice in l_points_generatrices:
201 generatrice_i = MakeInterpol(points_generatrice)
202 l_generatrices.append(generatrice_i)
204 # Pour cloturer le contour
205 l_generatrices.append(l_generatrices[0])
206 compound_generatrices = MakeCompound(l_generatrices)
207 face = MakeFilling(compound_generatrices, 0, 10, 1e-05, 1e-05, 0)
212 # ===============================================
216 # lecture des mesures dans le fichier csv
217 mesures_filename = os.path.join(dir_name, "mesure-transposee.csv")
218 d_generatrices, l_noms_generatrices, l_abscisses = read_generatrice(mesures_filename)
220 # lecture des angles dans le fichier csv
221 angles_filename = os.path.join(dir_name, "mesure-angles.csv")
222 l_angles = read_angles(angles_filename)
224 # dictionnaire indiquant les angles en fonction du nom de la generatrice
226 for nom, angle in zip(l_noms_generatrices, l_angles):
227 d_angles[nom] = angle
231 print "Temps de lecture des fichiers: %.1f s."%(time1-time0)
233 # Rq: pour conserver le point de plus faible epaisseur, il faut que la couture de la face de filling
234 # se situe sur la generatrice ou se situe ce point
236 # Points et vecteurs de base
238 P0 = MakeVertex(0, 0, 0)
240 Vx = MakeVectorDXDYDZ(1000, 0, 0)
241 Vy = MakeVectorDXDYDZ(0, 1000, 0)
242 Vz = MakeVectorDXDYDZ(0, 0, 1)
243 plane_size = longueur_tube * 10
245 P_ext = MakeVertex(0, r_ext, 0)
246 cercle_ext = MakeRevolution(P_ext, Vx, 2*math.pi)
247 face_ext_milieu = MakePrismVecH(cercle_ext, Vx, l_abscisses[-1])
248 addToStudy(face_ext_milieu, "face_ext_milieu")
250 # initialisation de l'epaisseur minimale pour l'algo de recherche de l'epaisseur minimale
252 # initialisation de la liste de generatrices
254 # initialisation de la double liste d'arcs
255 l_arcs_points = [[] for abscisse in l_abscisses]
257 # angle de la generatrice ou se situe l'epaisseur minimale
259 # Point ou se situe l'epaisseur minimale
261 # indice de la generatrice ou se situe le point d'epaisseur minimale
262 i_generatrice_ep_mini = None
263 # Creation des generatrices
264 for i, nom_generatrice in enumerate(l_noms_generatrices):
265 angle += d_angles[nom_generatrice]
266 l_ep = d_generatrices[nom_generatrice]
269 for ep, abscisse in zip(l_ep, l_abscisses):
270 P_generatrice_tmp = MakeVertex(abscisse, r_ext - ep*amplification_factor, 0)
271 P_generatrice = MakeRotation(P_generatrice_tmp, Vx, math.radians(angle))
272 #addToStudy(P_generatrice, "P_generatrice_%i"%(i+1))
273 # pour la methode par les generatrices
274 l_points_ep.append(P_generatrice)
275 # pour la methode par les arcs
276 l_arcs_points[j].append(P_generatrice)
277 # Test sur l'epaisseur minimale
280 P_ep_mini = P_generatrice
281 i_generatrice_ep_mini = i
282 angle_ep_mini = angle
284 # creation des generatrices
285 generatrice = MakeInterpol(l_points_ep)
286 addToStudy(generatrice, "Generatrice_%s"%nom_generatrice)
287 l_generatrices.append(generatrice)
289 print "epaisseur minimale mesuree: ", ep_min
290 addToStudy(P_ep_mini, "P_ep_mini")
292 # On cree un objet contenant tous les points pour voir si la face generee passe par tous les points
294 for arcs_points in l_arcs_points:
295 l_points += arcs_points
296 Tous_Points = MakeCompound(l_points)
297 addToStudy(Tous_Points, "Tous_Points")
299 # methode par les generatrices
300 # ============================
302 # Pour s'assurer de passer par le point d'epaisseur minimale,
303 # on decalle la generatrice de recollement de la face et on l'ajoute a la fin pour fermer la face
304 l_generatrices = l_generatrices[i_generatrice_ep_mini:] + l_generatrices[:i_generatrice_ep_mini] + [l_generatrices[i_generatrice_ep_mini]]
305 # Creation de la face
306 compound_generatrices = MakeCompound(l_generatrices)
307 addToStudy(compound_generatrices, "compound_generatrices")
308 min_compound_generatrices = MinDistance(face_ext_milieu, compound_generatrices)
309 print "epaisseur minimale entre les generatrices et la face exterieure: ", min_compound_generatrices
310 if methode == "generatrices":
311 face_int_milieu_tmp = MakeFilling(compound_generatrices, 0, 10, 1e-05, 1e-05, 0)
312 face_int_milieu = ChangeOrientation(face_int_milieu_tmp)
313 addToStudy(face_int_milieu, "face_int_milieu")
314 min_distance = MinDistance(face_ext_milieu, face_int_milieu)
315 print "epaisseur minimale avec la methode generatrices: ", min_distance
318 # methode par les arcs avec filling arc par arc, 3 points par points (methode ecaille de tortue)
319 # ==================================================================
321 if methode != "generatrices":
322 face_int_milieu = MakeShellFromPoints(l_arcs_points)
323 addToStudy(face_int_milieu, "face_int_milieu")
325 min_distance = MinDistance(face_ext_milieu, face_int_milieu)
326 print "epaisseur minimale avec la methode ecaille de tortue: ", min_distance
328 # => La face suit a la fois les generatrices et les arcs. L'epaisseur minimale est respectee
330 # Partie complementaire de la face interieure du tube
331 # ===================================================
333 # calcul de la longueur d'amortissement a partir de la formule de U4.PC.10-D
334 r_moyen = r_ext - ep_nominale/2.
335 l_amor_1 = 3/2.*math.sqrt(r_moyen**3/ep_nominale)
336 l_amor_2 = 3*math.sqrt(r_moyen*ep_nominale)
337 longueur_amortissement = max(l_amor_1, l_amor_2)
339 print "longueur d'amortissement: ", longueur_amortissement
341 # Longueur de transition entre tube deformé et longueur d'amortissement
342 longueur_transition = longueur_amortissement/5.
343 print "longueur de transition: ", longueur_transition
346 # On cree un nuage de points definissant les contours de la face de transition
347 r_int = r_ext - ep_nominale*amplification_factor
349 l_abscisses_transition = []
350 # boucle pour traiter en meme temps le prolongement en bas et en haut
351 for i_abscisse, coef_translation in zip([0, -1], [-1, 1]):
352 l_arcs_points_transition = []
354 # On cree les points sur le cercle de la face de transition, en vis-a-vis des points du premier cercle de la face int
355 abscisse_transition = l_abscisses[i_abscisse] + coef_translation*longueur_transition
356 l_abscisses_transition.append(abscisse_transition)
357 l_arcs_points_transition_base = []
359 for nom_generatrice in l_noms_generatrices:
360 angle += d_angles[nom_generatrice]
361 P_transition_base_tmp = MakeVertex(abscisse_transition, r_int, 0)
362 P_transition_base = MakeRotation(P_transition_base_tmp, Vx, math.radians(angle))
363 l_arcs_points_transition_base.append(P_transition_base)
366 l_arcs_points_transition.append(l_arcs_points_transition_base)
368 l_arcs_points_transition.append(l_arcs_points[i_abscisse])
369 if methode == "generatrices":
370 face_transition = MakeFillingFromPoints(l_arcs_points_transition)
371 if coef_translation == -1:
372 face_transition = ChangeOrientation(face_transition)
374 face_transition = MakeShellFromPoints(l_arcs_points_transition)
375 addToStudy(face_transition, "face_transition")
376 l_faces.append(face_transition)
379 # On recupere le contour bas pour creer la face d'amortissement
381 P_transition = MakeVertex(abscisse_transition, 0, 0)
382 l_edge_transition = GetShapesOnPlaneWithLocation(face_transition, ShapeType["EDGE"], Vx, P_transition, GEOM.ST_ON)
383 wire_bas_transition = MakeWire(l_edge_transition)
384 face_amortissement = MakePrismVecH(wire_bas_transition, Vx, coef_translation*longueur_amortissement)
385 addToStudy(face_amortissement, "face_amortissement")
386 l_faces.append(face_amortissement)
388 l_faces.append(face_int_milieu)
390 if methode == "generatrices":
391 face_int = MakeSewing(l_faces, 0.1)
393 face_int = MakeShell(l_faces)
394 addToStudy(face_int, "face_int")
397 # Creation du tube solide
398 # =======================
402 h_tube = l_abscisses[-1] - l_abscisses[0] + 2*(longueur_amortissement+longueur_transition)
403 abscisse_base_tube = l_abscisses[0]-(longueur_amortissement+longueur_transition)
404 cercle_ext_bas = MakeTranslation(cercle_ext, abscisse_base_tube, 0, 0)
405 face_ext_tmp = MakePrismVecH(cercle_ext_bas, Vx, h_tube)
406 # on tourne la face, pour ne pas avoir l'edge de couture
407 face_ext = MakeRotation(face_ext_tmp, Vx, math.radians(l_angles[0]))
408 #face_ext = face_ext_tmp
409 addToStudy(face_ext, "face_ext")
412 P_bas_ext = MakeTranslation(P_ext, abscisse_base_tube, 0, 0)
413 cercle_int_bas = CreateGroup(face_int, ShapeType["EDGE"])
414 l_cercle_int_bas = GetShapesOnPlaneWithLocation(face_int, ShapeType["EDGE"], Vx, P_bas_ext, GEOM.ST_ON)
415 UnionList(cercle_int_bas, l_cercle_int_bas)
417 face_bas = MakeFaceWires([cercle_ext_bas, cercle_int_bas], 1)
418 addToStudy(face_bas, "face_bas")
421 P_haut_ext = MakeTranslation(P_bas_ext, h_tube, 0, 0)
422 cercle_ext_haut = MakeTranslation(cercle_ext_bas, h_tube, 0, 0)
423 cercle_int_haut = CreateGroup(face_int, ShapeType["EDGE"])
424 l_cercle_int_haut = GetShapesOnPlaneWithLocation(face_int, ShapeType["EDGE"], Vx, P_haut_ext, GEOM.ST_ON)
425 UnionList(cercle_int_haut, l_cercle_int_haut)
427 face_haut = MakeFaceWires([cercle_ext_haut, cercle_int_haut], 1)
428 addToStudy(face_haut, "face_haut")
430 l_faces_tube = [face_int, face_ext, face_bas, face_haut]
431 shell_tube = MakeShell(l_faces_tube)
432 addToStudy(shell_tube, "shell_tube")
434 tube = MakeSolid([shell_tube])
435 addToStudy(tube, "tube")
439 print "Temps de creation de la geometrie: %.1f s."%(time2-time1)
442 # Partition pour que le maillage passe par les points de mesure
443 # =============================================================
445 l_plans_abscisses = []
448 if methode == "generatrices":
449 l_abscisses_plan = [l_abscisses_transition[0]] + [l_abscisses[0], l_abscisses[-1]] + [l_abscisses_transition[1]]
451 l_abscisses_plan = [l_abscisses_transition[0]] + l_abscisses[:] + [l_abscisses_transition[1]]
453 # un plan par abscisse
454 for abscisse in l_abscisses_plan:
455 P_plan_part = MakeVertex(abscisse, 0, 0)
456 plan_part = MakePlane(P_plan_part, Vx, plane_size)
457 l_plans_abscisses.append(plan_part)
459 P_axe_tube = MakeVertex(abscisse_base_tube, 0, 0)
460 axe_tube = MakePrismVecH(P_axe_tube, Vx, h_tube)
462 # pour rabotter les plans de partition
463 P_bas = MakeVertex(abscisse_base_tube, 0, 0)
464 cylindre_int = MakeCylinder(P_bas, Vx, r_int/2., h_tube)
468 # un plan sur toutes les generatrices
469 # Rq: si on cree un plan toutes les 2 generatrices uniquement,
470 # le maillage ne passera pas par les points de mesure de l'autre generatrice sur 2
471 for i, nom_generatrice in enumerate(l_noms_generatrices):
472 angle += d_angles[nom_generatrice]
473 if (methode != "generatrices") or (methode == "generatrices" and (i==i_generatrice_ep_mini)):
474 # TODO: lorsque MakePartition fonctionnera mieux (!),
475 # supprimer le if ci-dessus pour toujours partitionner par les plans des generatrices
476 P_vec_plan_tmp = MakeVertex(0, r_int, 0)
477 P_vec_plan = MakeRotation(P_vec_plan_tmp, Vx, math.radians(angle))
478 V_plan = MakeVector(P0, P_vec_plan)
479 plan_diag_tmp = MakePrismVecH(axe_tube, V_plan, 2*r_ext)
480 plan_diag = MakeCut(plan_diag_tmp, cylindre_int)
481 l_plans_diag.append(plan_diag)
483 # TODO: lorsque MakePartition fonctionnera mieux (!), supprimer ce bloc
484 # car on aura partitionne par toutes les generatrices dans le bloc precedent.
485 if methode == "generatrices":
486 plan_oppose = MakeRotation(l_plans_diag[-1], Vx, math.pi)
487 l_plans_diag.append(plan_oppose)
489 tous_plans = MakeCompound(l_plans_abscisses + l_plans_diag)
490 addToStudy(tous_plans, "tous_plans")
492 plans_diag = MakeCompound(l_plans_diag)
493 addToStudy(plans_diag, "plans_diag")
495 plans_abscisses = MakeCompound(l_plans_abscisses)
496 addToStudy(plans_abscisses, "plans_abscisses")
498 if type_maillage == "regle":
499 # Partion d'un tube plein par la face_int
500 cylindre_tmp = MakeCylinder(P_bas, Vx, r_ext, h_tube)
501 # le cylindre ainsi cree a sa ligne de couture sur Z
502 # => on la decalle sur l'edge de couture de la face interieure: -pi/2 + angle_ep_mini
503 if methode == "generatrices":
504 cylindre = MakeRotation(cylindre_tmp, Vx, -math.pi/2. + math.radians(angle_ep_mini))
506 # en methode ecailles de tortue, la reparation plante si on tourne le cylindre
507 # mais reussi si le cylindre reste avec des edges de couture non paralleles!
508 cylindre = cylindre_tmp
509 addToStudy(cylindre, "cylindre")
511 cylindre_part = MakePartition([cylindre, face_int])
512 addToStudy(cylindre_part, "cylindre_part")
514 # on recupere le solide correspondant au tube
515 P_tube = MakeVertex(abscisse_base_tube, (r_ext+r_int)/2., 0)
516 tube = GetBlockNearPoint(cylindre_part, P_tube)
517 addToStudy(tube, "tube")
519 if methode == "generatrices":
520 # partition par plans diag puis plans abscisses
521 tube_part_tmp = MakePartition([tube], [plans_diag])
522 addToStudy(tube_part_tmp, "tube_part_tmp")
523 tube_part = MakePartition([tube_part_tmp], [plans_abscisses])
525 # partition par plans abscisses puis plans diag
526 tube_part_tmp = MakePartition([tube], [plans_abscisses])
527 addToStudy(tube_part_tmp, "tube_part_tmp")
528 tube_part = MakePartition([tube_part_tmp], [plans_diag])
529 addToStudy(tube_part, "tube_part")
531 tube_part_improved = CheckAndImprove(tube_part)
532 if not tube_part_improved:
533 print "pas de reparation effectuee"
534 tube_part_improved = tube_part
535 addToStudy(tube_part_improved, "tube_part_improved")
538 # on partitionne les faces du bas et du haut pour pouvoir imposer le nombre de segments dans l'epaisseur
539 face_bas_part = MakePartition([face_bas], [tous_plans], Limit=ShapeType["FACE"])
540 face_haut_part = MakePartition([face_haut], [tous_plans], Limit=ShapeType["FACE"])
541 # pour le maillage libre, on partitionne uniquement avec les points de mesure
542 # pour qu'ils soient contenus dans le maillage
543 if methode == "generatrices":
544 # pour la methode generatrices, il faut partitionner la face interieure
545 face_int_part = MakePartition([face_int], [tous_plans], Limit=ShapeType["FACE"])
547 # pour la methode tortue, il suffit de partitionner par les points
548 # (en fait, seuls manquent les points au milieu des arcs)
549 face_int_part = MakePartition([face_int], [Tous_Points], Limit=ShapeType["FACE"])
551 l_faces_tube = [face_int_part, face_ext, face_bas_part, face_haut_part]
552 shell_tube = MakeShell(l_faces_tube)
553 addToStudy(shell_tube, "shell_tube")
555 tube_part_improved = MakeSolid([shell_tube])
556 addToStudy(tube_part_improved, "tube_part_improved")
560 print "Temps de partitionnement: %.1f s."%(time3-time2)
562 # Sous-geometries pour les sous-maillages
563 # =======================================
565 # edges dans l'epaisseur
566 l_edges_bas = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["EDGE"], Vx, P_bas, GEOM.ST_ON)
567 edges_bas = CreateGroup(tube_part_improved, ShapeType["EDGE"])
568 UnionList(edges_bas, l_edges_bas)
569 edges_ep_bas = GetEdgesByLength (edges_bas, 0, r_ext-r_int + 1e-1)
570 addToStudyInFather(tube_part_improved, edges_ep_bas, "edges_ep_bas")
572 if type_maillage == "libre":
573 # on recupere les faces bas et haut
574 l_faces_bas = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_bas, GEOM.ST_ON)
575 l_faces_haut = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_haut_ext, GEOM.ST_ON)
576 faces_extremites = CreateGroup(tube_part_improved, ShapeType["FACE"])
577 UnionList(faces_extremites, l_faces_bas + l_faces_haut)
578 addToStudyInFather(tube_part_improved, faces_extremites, "faces_extremites")
581 edges_arc_bas = CreateGroup(edges_bas, ShapeType["EDGE"])
582 l_edges_ep_bas = SubShapeAllSorted(edges_ep_bas, ShapeType["EDGE"])
583 UnionList(edges_arc_bas, l_edges_bas)
584 DifferenceList(edges_arc_bas, l_edges_ep_bas)
585 addToStudyInFather(tube_part_improved, edges_arc_bas, "edges_arc_bas")
587 # on recupere la face interieure
588 l_face_int = GetShapesOnCylinder(tube_part_improved, ShapeType["FACE"], Vx, r_ext - ep_min, GEOM.ST_IN)
589 sub_face_int = CreateGroup(tube_part_improved, ShapeType["FACE"])
590 UnionList(sub_face_int, l_face_int)
591 addToStudyInFather(tube_part_improved, sub_face_int, "SURF_INT")
593 # on recupere les edges d'amortissement
594 P_bas_int = MakeTranslation(P_bas, 0, r_int, 0)
595 P_edge_tmp = MakeRotation(P_bas_int, Vx, math.radians(angle_ep_mini))
596 P_edge_amortissement_1 = MakeTranslation(P_edge_tmp, longueur_amortissement/2., 0, 0)
597 P_edge_amortissement_2 = MakeTranslation(P_edge_tmp, h_tube-longueur_amortissement/2., 0, 0)
598 edge_amortissement_1 = GetEdgeNearPoint(tube_part_improved, P_edge_amortissement_1)
599 edge_amortissement_2 = GetEdgeNearPoint(tube_part_improved, P_edge_amortissement_2)
600 edges_amortissement = CreateGroup(tube_part_improved, ShapeType["EDGE"])
601 UnionList(edges_amortissement, [edge_amortissement_1, edge_amortissement_2])
602 addToStudyInFather(tube_part_improved, edges_amortissement, "edges_amortissement")
604 # on recupere les edges de transition
605 P_edge_transition_1 = MakeTranslation(P_edge_tmp, longueur_transition/2. + longueur_amortissement, 0, 0)
606 P_edge_transition_2 = MakeTranslation(P_edge_tmp, h_tube-(longueur_transition/2. + longueur_amortissement), 0, 0)
607 edge_transition_1 = GetEdgeNearPoint(tube_part_improved, P_edge_transition_1)
608 edge_transition_2 = GetEdgeNearPoint(tube_part_improved, P_edge_transition_2)
609 edges_transition = CreateGroup(tube_part_improved, ShapeType["EDGE"])
610 UnionList(edges_transition, [edge_transition_1, edge_transition_2])
611 addToStudyInFather(tube_part_improved, edges_transition, "edges_transition")
613 # on recupere les edges d'une generatrice
615 axe_generatrice_tmp = MakeTranslation(Vx, 0, r_int, 0)
616 axe_generatrice = MakeRotation(axe_generatrice_tmp, Vx, math.radians(angle_ep_mini))
617 l_edges_generatrice = GetShapesOnCylinder(tube_part_improved, ShapeType["EDGE"], axe_generatrice, ep_nominale, GEOM.ST_IN)
618 edges_generatrice = CreateGroup(tube_part_improved, ShapeType["EDGE"])
619 UnionList(edges_generatrice, l_edges_generatrice)
620 DifferenceList(edges_generatrice, [edge_amortissement_1, edge_amortissement_2]+[edge_transition_1, edge_transition_2])
621 addToStudyInFather(tube_part_improved, edges_generatrice, "edges_generatrice")
623 # on recupere les edges d'une generatrice dont la longueur est inferieure a petite_distance
624 # pour pouvoir imposer moins de segments
625 if methode != "generatrices":
626 edges_petite_distance = CreateGroup(tube_part_improved, ShapeType["EDGE"])
627 l_petite_distance = []
628 for edge in l_edges_generatrice:
629 length = BasicProperties(edge)[0]
630 if length <= petite_distance:
631 l_petite_distance.append(edge)
632 UnionList(edges_petite_distance, l_petite_distance)
633 addToStudyInFather(tube_part_improved, edges_petite_distance, "edges_petite_distance")
637 # Sous-geometries pour les groupes
638 # ================================
640 # on recupere la face interieure sans les zones de transition et d'amortissement
641 P_abs_first = MakeVertex(l_abscisses[0], 0, 0)
642 l_face_int_inf = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_abs_first, GEOM.ST_IN)
643 P_abs_last = MakeVertex(l_abscisses[-1], 0, 0)
644 l_face_int_sup = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_abs_last, GEOM.ST_OUT)
645 sub_face_int_milieu = CreateGroup(tube_part_improved, ShapeType["FACE"])
646 UnionList(sub_face_int_milieu, l_face_int)
647 DifferenceList(sub_face_int_milieu, l_face_int_inf + l_face_int_sup)
648 addToStudyInFather(tube_part_improved, sub_face_int_milieu, "FaceIntM")
650 # on recupere la face exterieure
651 l_face_ext = GetShapesOnCylinder(tube_part_improved, ShapeType["FACE"], Vx, r_ext, GEOM.ST_ON)
652 sub_face_ext = CreateGroup(tube_part_improved, ShapeType["FACE"])
653 UnionList(sub_face_ext, l_face_ext)
654 addToStudyInFather(tube_part_improved, sub_face_ext, "SURF_EXT")
656 # on recupere la face a l'extremite amont du tube
657 l_face_bas = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_bas, GEOM.ST_ON)
658 sub_face_bas = CreateGroup(tube_part_improved, ShapeType["FACE"])
659 UnionList(sub_face_bas, l_face_bas)
660 addToStudyInFather(tube_part_improved, sub_face_bas, "CLGV")
662 # on recupere la face a l'extremite aval du tube
663 l_face_haut = GetShapesOnPlaneWithLocation(tube_part_improved, ShapeType["FACE"], Vx, P_haut_ext, GEOM.ST_ON)
664 sub_face_haut = CreateGroup(tube_part_improved, ShapeType["FACE"])
665 UnionList(sub_face_haut, l_face_haut)
666 addToStudyInFather(tube_part_improved, sub_face_haut, "EXTUBE")
668 # On recupere les edges communs a face_int et face_haut
669 l_edge_int = GetShapesOnCylinderIDs(tube_part_improved, ShapeType["EDGE"], Vx, r_ext - ep_min, GEOM.ST_IN)
670 l_edge_haut = GetShapesOnPlaneWithLocationIDs(tube_part_improved, ShapeType["EDGE"], Vx, P_haut_ext, GEOM.ST_ON)
671 l_edge_bord_int_haut = []
672 for id_edge in l_edge_int:
673 if id_edge in l_edge_haut:
674 l_edge_bord_int_haut.append(id_edge)
676 edge_bord_int_haut = CreateGroup(tube_part_improved, ShapeType["EDGE"])
677 UnionIDs(edge_bord_int_haut, l_edge_bord_int_haut)
678 addToStudyInFather(tube_part_improved, edge_bord_int_haut, "BORDTU")
680 # on recupere le point d'epaisseur minimale
681 # avec la methode tortue le maillage passe forcement par le point d'epaisseur minimale
682 #if methode != "generatrices":
683 x, y, z = PointCoordinates(P_ep_mini)
684 P_ep_mini_sub = GetPoint(tube_part_improved, x, y, z, 1e-5)
685 addToStudyInFather(tube_part_improved, P_ep_mini_sub, "P_ep_mini")
689 print "Temps de recuperation des sous-geometries: %.1f s."%(time4-time3)
694 # on divise le nombre de segments par le nombre de plans suivant les abscisses
696 Maillage = smesh.Mesh(tube_part_improved, "Tube")
698 if type_maillage == "regle":
699 algo1D = Maillage.Segment()
700 algo1D.NumberOfSegments(nb_seg_generatrices)
701 #algo1D.QuadraticMesh()
703 Maillage.Quadrangle()
705 Maillage.Hexahedron()
708 average_length = h_tube/(nb_seg_generatrices + 2*(nb_seg_amortissement+nb_seg_transition))/5.
709 # BLSURF est un algo 1D/2D
710 if type_algo == "BLSURF":
711 algo2D = Maillage.Triangle(algo=smesh.BLSURF)
712 algo2D.SetPhySize(average_length)
714 algo1D = Maillage.Segment()
715 algo1D.LocalLength(average_length)
717 Maillage.Triangle(smesh.NETGEN_2D)
720 #algo3D = Maillage.Tetrahedron(smesh.GHS3D)
721 algo3D = Maillage.Tetrahedron()
727 # On impose finalement un maillage fin partout, seul moyen d'avoir plusieurs elements dans l'epaisseur
728 if type_maillage == "libre":
730 average_length_extremites = average_length/8.
731 # BLSURF est un aglo 1D/2D
732 if type_algo == "BLSURF":
733 #algo2D = Maillage.Triangle(geom=faces_extremites, algo=smesh.BLSURF)
734 #algo2D.SetPhySize(average_length_extremites)
738 algo1D = Maillage.Segment(faces_extremites)
739 algo1D.LocalLength(average_length_extremites)
741 if type_maillage == "regle":
743 algo1D = Maillage.Segment(edges_ep_bas)
744 algo1D.NumberOfSegments(nb_seg_ep)
748 algo1D = Maillage.Segment(edges_arc_bas)
749 algo1D.NumberOfSegments(nb_seg_arc)
752 # sur les longueurs d'amortissement
753 algo1D = Maillage.Segment(edges_amortissement)
754 algo1D.NumberOfSegments(nb_seg_amortissement)
757 # sur les longueurs de transition
758 algo1D = Maillage.Segment(edges_transition)
759 algo1D.NumberOfSegments(nb_seg_transition)
762 if methode == "tortue":
763 algo1D = Maillage.Segment(edges_petite_distance)
764 algo1D.NumberOfSegments(nb_seg_petites_distances)
769 # on fait passer le maillage par le point de plus faible epaisseur
770 #if methode == "generatrices":
771 #x, y, z = PointCoordinates(P_ep_mini)
772 #id_node = Maillage.MeshToPassThroughAPoint(x, y, z)
773 ## on cree le groupe avec le point de plus faible epaisseur
774 #Maillage.MakeGroupByIds("P_ep_mini", smesh.NODE, [id_node])
776 # on a deja cree le groupe geometrique
777 #if methode != "generatrices":
778 Maillage.Group(P_ep_mini_sub)
780 # conversion en quadratique (tres long à afficher)
781 #Maillage.ConvertToQuadratic(1)
783 # on ajoute deux points aux extremites de l'axe du tube
784 x, y, z = PointCoordinates(P_bas)
785 id_p2 = Maillage.AddNode(x, y, z)
786 Maillage.MakeGroupByIds("P2", smesh.NODE, [id_p2])
788 id_p1 = Maillage.AddNode(x+h_tube, y, z)
789 Maillage.MakeGroupByIds("P1", smesh.NODE, [id_p1])
794 Maillage.Group(edge_bord_int_haut)
795 Maillage.Group(sub_face_int)
796 Maillage.Group(sub_face_int, "PEAUINT")
797 Maillage.Group(sub_face_ext)
798 Maillage.Group(sub_face_ext, "PEAUEXT")
799 Maillage.Group(sub_face_bas)
800 Maillage.Group(sub_face_bas, "FACE1")
801 Maillage.Group(sub_face_haut)
802 Maillage.Group(sub_face_haut, "FACE2")
803 group_nodes_int_milieu = Maillage.GroupOnGeom(sub_face_int_milieu, "NodesInt", smesh.NODE)
804 Maillage.GroupOnGeom(tube_part_improved, "VOL_TUBE", smesh.VOLUME)
805 Maillage.GroupOnGeom(tube_part_improved, "COUDE", smesh.VOLUME)
809 ([MaillageGibi], status) = smesh.smesh.CreateMeshesFromMED(maillage_gibi)
811 # on met le maillage gibi dans le meme axe que le maillage SALOME
812 MaillageGibi.RotateObject(MaillageGibi, Vz, -math.pi/2., False)
814 V_trans = MakeVectorDXDYDZ(-longueur_amortissement-longueur_transition+l_abscisses[0], 0, 0)
815 MaillageGibi.TranslateObject(MaillageGibi, V_trans, False)
817 MaillageInt = MaillageGibi
819 gibi_groupes = MaillageGibi.GetGroups()
821 # on determine le groupe correspondant a la face interieure
823 for groupe in gibi_groupes:
824 name = groupe.GetName()
825 if name.strip() == "SURF_INT":
828 l_faces_int = group_int.GetIDs()
830 for face_id in l_faces_int:
831 l_nodes = MaillageGibi.GetElemNodes(face_id)
833 if node not in l_nodes_ids:
834 l_nodes_ids.append(node)
838 print "Temps de generation du maillage: %.1f s."%(time5-time4)
840 # Verifions si le maillage passe par les points de mesure
841 # =======================================================
843 # on cree un maillage avec les points de mesure
844 MaillageTousPoints = smesh.Mesh(Tous_Points, "Tous_Points")
845 ## BUG: smesh ne peut pas creer un maillage de points!
846 #MaillageTousPoints.Compute()
847 # => On ajoute les points un par un...
848 l_tous_points = SubShapeAllSorted(Tous_Points, ShapeType["VERTEX"])
849 for point in l_tous_points:
850 x, y, z = PointCoordinates(point)
851 MaillageTousPoints.AddNode(x, y, z)
853 l_points_mesures_ids = MaillageTousPoints.GetNodesId()
855 # on ajoute les noeuds mailles de la face interieure
857 l_nodes_ids = group_nodes_int_milieu.GetIDs()
858 MaillageInt = Maillage
859 for node in l_nodes_ids:
860 # on recupere les coordonnees depuis le maillage global
861 x, y, z = MaillageInt.GetNodeXYZ(node)
862 # on ajoute ce noeud dans le maillage de points
863 MaillageTousPoints.AddNode(x, y, z)
865 # on trouve les noeuds en double
867 coincident_nodes = MaillageTousPoints.FindCoincidentNodes(tolerance)
869 # nombre de points de mesure
870 nb_points = len(l_points)
872 # nombre de noeuds en commun
873 nb_coincident_nodes = len(coincident_nodes)
875 # nombre de points perdus
876 nb_points_perdus = nb_points - nb_coincident_nodes
878 print "%i/%i points de mesure ont ete conserves dans le maillage"%(nb_coincident_nodes, nb_points)
880 print "%i/%i points de mesure ont ete perdus dans le maillage"%(nb_points_perdus, nb_points)
882 # affichage des points de mesure conserves
883 group_kept_nodes = MaillageTousPoints.CreateEmptyGroup(smesh.NODE, "Kept_measure_points")
884 l_id_coincident_nodes = []
885 for l_id in coincident_nodes:
886 l_id_coincident_nodes += l_id
887 group_kept_nodes.Add(l_id_coincident_nodes)
889 # affichage des points de mesure perdus
890 group_lost_nodes = MaillageTousPoints.CreateEmptyGroup(smesh.NODE, "Lost_measure_points")
891 l_id_lost_points = []
892 for id_point in l_points_mesures_ids:
893 if id_point not in l_id_coincident_nodes:
894 l_id_lost_points.append(id_point)
895 group_lost_nodes.Add(l_id_lost_points)
897 # On merge les noeuds en double
899 MaillageTousPoints.MergeNodes(coincident_nodes)
901 # on regarde si le point d'epaisseur minimale fait partie des points garde
902 x_mini, y_mini, z_mini = PointCoordinates(P_ep_mini)
903 id_p_ep_mini = MaillageTousPoints.AddNode(x_mini, y_mini, z_mini)
904 MaillageTousPoints.MakeGroupByIds("P_ep_mini", smesh.NODE, [id_p_ep_mini])
906 coincident_nodes_ep_mini = MaillageTousPoints.FindCoincidentNodes(tolerance)
907 if coincident_nodes_ep_mini:
908 print "Le point d'epaisseur minimale a ete conserve"
910 print "Le point d'epaisseur minimale a ete perdu"
914 print "Temps de verification du maillage: %.1f s."%(time6-time5)
916 print "Temps total: %.1f s."%(time6 - time0)
918 salome.sg.updateObjBrowser(0)