generator/tube.py

   1
   2 # ==========
   3 # PARAMETRES
   4 # ==========
   5
   6 # methode de filling: par les generatrices ou par des ecailles de tortue
   7
   8 # Longueur
   9 longueur_tube = 6363.655
  10
  11 # Facteur d'amplification pour voir si le filling passe bien par les generatrices
  12 amplification_factor = 1
  13
  14 # teste si le maillage gibi passe par tous les points
  15 test_gibi = False
  16 maillage_gibi = '/home/PROJETS/OUVERT/edf_anode/from_edf/2008_11_28/Tube/tube.mail.med'
  17
  18 # MAILLAGE
  19 # ========
  20
  21 type_maillage = "regle" # "regle" ou "libre"
  22
  23 if type_maillage == "libre":
  24      type_algo = "NETGEN_2D" # "BLSURF" ou "NETGEN_2D"
  25
  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
  29 else:
  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
  34     petite_distance = 100
  35     nb_seg_petites_distances = 3
  36
  37 # nombre de segments dans l'epaisseur du tube
  38 nb_seg_ep= 3
  39
  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
  43 else:
  44     nb_seg_arc = 5 # methode tortue: partition en 8
  45
  46 # nombre de segments dans la longueur de transition
  47 nb_seg_transition = 4
  48
  49 # nombre de segments dans la longueur d'amortissement
  50 nb_seg_amortissement = 11
  51
  52 # ===========
  53
  54 import salome
  55 from geompy import *
  56 import smesh
  57
  58 import os, csv, math, pdb, time
  59
  60 # DEBUT DES FONCTIONS
  61 # ===================================================
  62
  63 # supprime recursivement un element dans une liste
  64 # @exemple :
  65 #    l=["a", " ", "b", " ", "c"]
  66 #    l = recursive_remove(l, " ")
  67 #    print l
  68 # => ["a", "b", "c"]
  69 def recursive_remove(l, txt):
  70     finished = 0
  71     while not finished:
  72         try:
  73             l.remove(txt)
  74         except ValueError:
  75             finished = 1
  76     return l
  77
  78
  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")
  84
  85     reader = csv.reader(file, delimiter=';', quoting=csv.QUOTE_NONE)
  86
  87     # Dictionnaire pour stocker les mesures sur une generatrice donnee
  88     d_generatrices = {}
  89
  90     # Liste des noms des generatrices
  91     l_noms_generatrices = []
  92
  93     # Liste des abscisses des points de mesures
  94     l_abscisses = []
  95
  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
 100         if i==0:
 101             for nom in row:
 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)
 108         else:
 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)
 121
 122     file.close()
 123
 124     return d_generatrices, l_noms_generatrices, l_abscisses
 125
 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")
 130
 131     reader = csv.reader(file, delimiter=';', quoting=csv.QUOTE_NONE)
 132
 133     # Liste des angles des generatrices
 134     l_angles = []
 135
 136     for row in reader:
 137         # On supprime les cases vides
 138         row = recursive_remove(row, "")
 139
 140         # si la ligne comporte 3 valeurs, on peut lire les angles
 141         if len(row) == 3:
 142             angle_str = row[2]
 143             angle = float(angle_str)
 144             l_angles.append(angle)
 145     return l_angles
 146
 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)]
 157     if closed_wire:
 158         if nb_generatrices%2 != 1:
 159             raise "L'algo ne fonctionne pour l'instant qu'avec un nombre de generatrices impair"
 160     else:
 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:
 165         if not closed_wire:
 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)
 172
 173
 174     # on fait un filling pour tous les 1ers arcs, tous les 2emes arcs, ... jusqu'au 4e arc
 175     l_faces_filling = []
 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)
 185     return face
 186
 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)]
 194     l_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)
 203     if not closed_wire:
 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)
 208     return face
 209
 210
 211 # FIN DES FONCTIONS
 212 # ===============================================
 213
 214 time0 = time.time()
 215
 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)
 219
 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)
 223
 224 # dictionnaire indiquant les angles en fonction du nom de la generatrice
 225 d_angles = {}
 226 for nom, angle in zip(l_noms_generatrices, l_angles):
 227     d_angles[nom] = angle
 228
 229 time1 = time.time()
 230
 231 print "Temps de lecture des fichiers: %.1f s."%(time1-time0)
 232
 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
 235
 236 # Points et vecteurs de base
 237
 238 P0 = MakeVertex(0, 0, 0)
 239
 240 Vx = MakeVectorDXDYDZ(1000, 0, 0)
 241 Vy = MakeVectorDXDYDZ(0, 1000, 0)
 242 Vz = MakeVectorDXDYDZ(0, 0, 1)
 243 plane_size = longueur_tube * 10
 244
 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")
 249
 250 # initialisation de l'epaisseur minimale pour l'algo de recherche de l'epaisseur minimale
 251 ep_min = r_ext
 252 # initialisation de la liste de generatrices
 253 l_generatrices = []
 254 # initialisation de la double liste d'arcs
 255 l_arcs_points = [[] for abscisse in l_abscisses]
 256 angle = 0
 257 # angle de la generatrice ou se situe l'epaisseur minimale
 258 angle_ep_mini = 0
 259 # Point ou se situe l'epaisseur minimale
 260 P_ep_mini = None
 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]
 267     l_points_ep = []
 268     j = 0
 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
 278         if ep < ep_min:
 279             ep_min = ep
 280             P_ep_mini = P_generatrice
 281             i_generatrice_ep_mini = i
 282             angle_ep_mini = angle
 283         j+=1
 284     # creation des generatrices
 285     generatrice = MakeInterpol(l_points_ep)
 286     addToStudy(generatrice, "Generatrice_%s"%nom_generatrice)
 287     l_generatrices.append(generatrice)
 288
 289 print "epaisseur minimale mesuree: ", ep_min
 290 addToStudy(P_ep_mini, "P_ep_mini")
 291
 292 # On cree un objet contenant tous les points pour voir si la face generee passe par tous les points
 293 l_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")
 298
 299 # methode par les generatrices
 300 # ============================
 301
 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
 316
 317
 318 # methode par les arcs avec filling arc par arc, 3 points par points (methode ecaille de tortue)
 319 # ==================================================================
 320
 321 if methode != "generatrices":
 322     face_int_milieu = MakeShellFromPoints(l_arcs_points)
 323     addToStudy(face_int_milieu, "face_int_milieu")
 324
 325     min_distance = MinDistance(face_ext_milieu, face_int_milieu)
 326     print "epaisseur minimale avec la methode ecaille de tortue: ", min_distance
 327
 328 # => La face suit a la fois les generatrices et les arcs. L'epaisseur minimale est respectee
 329
 330 # Partie complementaire de la face interieure du tube
 331 # ===================================================
 332
 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)
 338
 339 print "longueur d'amortissement: ", longueur_amortissement
 340
 341 # Longueur de transition entre tube deformé et longueur d'amortissement
 342 longueur_transition = longueur_amortissement/5.
 343 print "longueur de transition: ", longueur_transition
 344
 345
 346 # On cree un nuage de points definissant les contours de la face de transition
 347 r_int = r_ext - ep_nominale*amplification_factor
 348 l_faces = []
 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 = []
 353
 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 = []
 358     angle = 0
 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)
 364
 365     # contour bas
 366     l_arcs_points_transition.append(l_arcs_points_transition_base)
 367     # contour haut
 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)
 373     else:
 374         face_transition = MakeShellFromPoints(l_arcs_points_transition)
 375     addToStudy(face_transition, "face_transition")
 376     l_faces.append(face_transition)
 377
 378
 379     # On recupere le contour bas pour creer la face d'amortissement
 380
 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)
 387
 388 l_faces.append(face_int_milieu)
 389
 390 if methode == "generatrices":
 391     face_int = MakeSewing(l_faces, 0.1)
 392 else:
 393     face_int = MakeShell(l_faces)
 394 addToStudy(face_int, "face_int")
 395
 396
 397 # Creation du tube solide
 398 # =======================
 399
 400 # Face exterieure
 401
 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")
 410
 411 # Face bas
 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)
 416
 417 face_bas = MakeFaceWires([cercle_ext_bas, cercle_int_bas], 1)
 418 addToStudy(face_bas, "face_bas")
 419
 420 # Face haut
 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)
 426
 427 face_haut = MakeFaceWires([cercle_ext_haut, cercle_int_haut], 1)
 428 addToStudy(face_haut, "face_haut")
 429
 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")
 433
 434 tube = MakeSolid([shell_tube])
 435 addToStudy(tube, "tube")
 436
 437 time2 = time.time()
 438
 439 print "Temps de creation de la geometrie: %.1f s."%(time2-time1)
 440
 441
 442 # Partition pour que le maillage passe par les points de mesure
 443 # =============================================================
 444
 445 l_plans_abscisses = []
 446
 447
 448 if methode == "generatrices":
 449     l_abscisses_plan = [l_abscisses_transition[0]] + [l_abscisses[0], l_abscisses[-1]] + [l_abscisses_transition[1]]
 450 else:
 451     l_abscisses_plan = [l_abscisses_transition[0]] + l_abscisses[:] + [l_abscisses_transition[1]]
 452
 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)
 458
 459 P_axe_tube = MakeVertex(abscisse_base_tube, 0, 0)
 460 axe_tube = MakePrismVecH(P_axe_tube, Vx, h_tube)
 461
 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)
 465
 466 angle = 0
 467 l_plans_diag = []
 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)
 482
 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)
 488
 489 tous_plans = MakeCompound(l_plans_abscisses + l_plans_diag)
 490 addToStudy(tous_plans, "tous_plans")
 491
 492 plans_diag = MakeCompound(l_plans_diag)
 493 addToStudy(plans_diag, "plans_diag")
 494
 495 plans_abscisses = MakeCompound(l_plans_abscisses)
 496 addToStudy(plans_abscisses, "plans_abscisses")
 497
 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))
 505     else:
 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")
 510
 511     cylindre_part = MakePartition([cylindre, face_int])
 512     addToStudy(cylindre_part, "cylindre_part")
 513
 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")
 518
 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])
 524     else:
 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")
 530
 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")
 536
 537 else:
 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"])
 546     else:
 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"])
 550
 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")
 554
 555     tube_part_improved = MakeSolid([shell_tube])
 556     addToStudy(tube_part_improved, "tube_part_improved")
 557
 558 time3 = time.time()
 559
 560 print "Temps de partitionnement: %.1f s."%(time3-time2)
 561
 562 # Sous-geometries pour les sous-maillages
 563 # =======================================
 564
 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")
 571
 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")
 579
 580 # edges sur les arcs
 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")
 586
 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")
 592
 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")
 603
 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")
 612
 613 # on recupere les edges d'une generatrice
 614
 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")
 622
 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")
 634
 635
 636
 637 # Sous-geometries pour les groupes
 638 # ================================
 639
 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")
 649
 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")
 655
 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")
 661
 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")
 667
 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)
 675
 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")
 679
 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")
 686
 687 time4 = time.time()
 688
 689 print "Temps de recuperation des sous-geometries: %.1f s."%(time4-time3)
 690
 691 # MAILLAGE
 692 # ========
 693
 694 # on divise le nombre de segments par le nombre de plans suivant les abscisses
 695
 696 Maillage = smesh.Mesh(tube_part_improved, "Tube")
 697
 698 if type_maillage == "regle":
 699     algo1D = Maillage.Segment()
 700     algo1D.NumberOfSegments(nb_seg_generatrices)
 701     #algo1D.QuadraticMesh()
 702
 703     Maillage.Quadrangle()
 704
 705     Maillage.Hexahedron()
 706 else:
 707     # 30
 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)
 713     else:
 714          algo1D = Maillage.Segment()
 715          algo1D.LocalLength(average_length)
 716
 717          Maillage.Triangle(smesh.NETGEN_2D)
 718
 719
 720     #algo3D = Maillage.Tetrahedron(smesh.GHS3D)
 721     algo3D = Maillage.Tetrahedron()
 722
 723
 724 # hypotheses locales
 725 # ==================
 726
 727 # On impose finalement un maillage fin partout, seul moyen d'avoir plusieurs elements dans l'epaisseur
 728 if type_maillage == "libre":
 729     # 8
 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)
 735          pass
 736     else:
 737         # 8
 738         algo1D = Maillage.Segment(faces_extremites)
 739         algo1D.LocalLength(average_length_extremites)
 740
 741 if type_maillage == "regle":
 742     # dans l'epaisseur
 743     algo1D = Maillage.Segment(edges_ep_bas)
 744     algo1D.NumberOfSegments(nb_seg_ep)
 745     algo1D.Propagation()
 746
 747     # sur les arcs
 748     algo1D = Maillage.Segment(edges_arc_bas)
 749     algo1D.NumberOfSegments(nb_seg_arc)
 750     algo1D.Propagation()
 751
 752     # sur les longueurs d'amortissement
 753     algo1D = Maillage.Segment(edges_amortissement)
 754     algo1D.NumberOfSegments(nb_seg_amortissement)
 755     algo1D.Propagation()
 756
 757     # sur les longueurs de transition
 758     algo1D = Maillage.Segment(edges_transition)
 759     algo1D.NumberOfSegments(nb_seg_transition)
 760     algo1D.Propagation()
 761
 762     if methode == "tortue":
 763         algo1D = Maillage.Segment(edges_petite_distance)
 764         algo1D.NumberOfSegments(nb_seg_petites_distances)
 765         algo1D.Propagation()
 766
 767 Maillage.Compute()
 768
 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])
 775
 776 # on a deja cree le groupe geometrique
 777 #if methode != "generatrices":
 778 Maillage.Group(P_ep_mini_sub)
 779
 780 # conversion en quadratique (tres long à afficher)
 781 #Maillage.ConvertToQuadratic(1)
 782
 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])
 787
 788 id_p1 = Maillage.AddNode(x+h_tube, y, z)
 789 Maillage.MakeGroupByIds("P1", smesh.NODE, [id_p1])
 790
 791 # Groupes
 792 # =========
 793
 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)
 806
 807 if test_gibi:
 808     time4 = time.time()
 809     ([MaillageGibi], status) = smesh.smesh.CreateMeshesFromMED(maillage_gibi)
 810
 811     # on met le maillage gibi dans le meme axe que le maillage SALOME
 812     MaillageGibi.RotateObject(MaillageGibi, Vz, -math.pi/2., False)
 813
 814     V_trans = MakeVectorDXDYDZ(-longueur_amortissement-longueur_transition+l_abscisses[0], 0, 0)
 815     MaillageGibi.TranslateObject(MaillageGibi, V_trans, False)
 816
 817     MaillageInt = MaillageGibi
 818
 819     gibi_groupes = MaillageGibi.GetGroups()
 820
 821     # on determine le groupe correspondant a la face interieure
 822     group_int = None
 823     for groupe in gibi_groupes:
 824         name = groupe.GetName()
 825         if name.strip() == "SURF_INT":
 826             group_int = groupe
 827             break
 828     l_faces_int = group_int.GetIDs()
 829     l_nodes_ids = []
 830     for face_id in l_faces_int:
 831         l_nodes = MaillageGibi.GetElemNodes(face_id)
 832         for node in l_nodes:
 833             if node not in l_nodes_ids:
 834                 l_nodes_ids.append(node)
 835
 836 time5 = time.time()
 837
 838 print "Temps de generation du maillage: %.1f s."%(time5-time4)
 839
 840 # Verifions si le maillage passe par les points de mesure
 841 # =======================================================
 842
 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)
 852
 853 l_points_mesures_ids = MaillageTousPoints.GetNodesId()
 854
 855 # on ajoute les noeuds mailles de la face interieure
 856 if not test_gibi:
 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)
 864
 865 # on trouve les noeuds en double
 866 tolerance = 1e0
 867 coincident_nodes = MaillageTousPoints.FindCoincidentNodes(tolerance)
 868
 869 # nombre de points de mesure
 870 nb_points = len(l_points)
 871
 872 # nombre de noeuds en commun
 873 nb_coincident_nodes = len(coincident_nodes)
 874
 875 # nombre de points perdus
 876 nb_points_perdus = nb_points - nb_coincident_nodes
 877
 878 print "%i/%i points de mesure ont ete conserves dans le maillage"%(nb_coincident_nodes, nb_points)
 879 if nb_points_perdus:
 880     print "%i/%i points de mesure ont ete perdus dans le maillage"%(nb_points_perdus, nb_points)
 881
 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)
 888
 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)
 896
 897 # On merge les noeuds en double
 898 if coincident_nodes:
 899     MaillageTousPoints.MergeNodes(coincident_nodes)
 900
 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])
 905
 906 coincident_nodes_ep_mini = MaillageTousPoints.FindCoincidentNodes(tolerance)
 907 if coincident_nodes_ep_mini:
 908     print "Le point d'epaisseur minimale a ete conserve"
 909 else:
 910     print "Le point d'epaisseur minimale a ete perdu"
 911
 912 time6 = time.time()
 913
 914 print "Temps de verification du maillage: %.1f s."%(time6-time5)
 915
 916 print "Temps total: %.1f s."%(time6 - time0)
 917
 918 salome.sg.updateObjBrowser(0)
 919