Aster/Cata/cataSTA10/Macro/post_gp_ops.py

   1 #@ MODIF post_gp_ops Macro  DATE 16/11/2009   AUTEUR COURTOIS M.COURTOIS
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2006  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20
  21 EnumTypes = (list, tuple)
  22
  23 # -----------------------------------------------------------------------------
  24 def post_gp_ops(self, **args):
  25    """
  26       Corps de la macro POST_GP
  27    """
  28    import pdb
  29    macro = 'POST_GP'
  30    ier=0
  31    from Accas import _F
  32    from Utilitai.Utmess       import UTMESS
  33    from Utilitai.Table        import Table, merge, Colonne
  34    from Cata_Utils.t_fonction import t_fonction
  35    from Cata.cata             import evol_noli
  36    import aster
  37
  38    # ----- On importe les definitions des commandes a utiliser dans la macro
  39    CALC_THETA    = self.get_cmd('CALC_THETA')
  40    CALC_G        = self.get_cmd('CALC_G')
  41    POST_ELEM     = self.get_cmd('POST_ELEM')
  42    POST_RELEVE_T = self.get_cmd('POST_RELEVE_T')
  43    CREA_TABLE    = self.get_cmd('CREA_TABLE')
  44    DEFI_LIST_ENTI= self.get_cmd('DEFI_LIST_ENTI')
  45    CALC_ELEM     = self.get_cmd('CALC_ELEM')
  46    RECU_FONCTION = self.get_cmd('RECU_FONCTION')
  47    DEFI_GROUP    = self.get_cmd('DEFI_GROUP')
  48    FIN           = self.get_cmd('FIN')
  49
  50    # ----- Comptage, commandes + déclaration concept sortant
  51    self.set_icmd(1)
  52    self.DeclareOut('result', self.sd)
  53    self.DeclareOut('tabresult', self['TABL_RESU'])
  54    info = self['INFO']
  55
  56    # 0. ----- Type de calcul
  57    identification = self['IDENTIFICATION'] != None
  58    if identification:
  59       # 0.1. --- identification : on boucle sur les valeurs de TEMP.
  60       #          Pour chaque couple (T, Kjc(T)), on évalue les Ki, Kmoy et
  61       #          les valeurs de Gpmax, DeltaLmax, inst.max correspondantes.
  62       mccalc = self['IDENTIFICATION']
  63       l_crit = mccalc['KJ_CRIT']
  64       lv_ident = []
  65       l_temp = mccalc['TEMP']
  66    else:
  67       # 0.2. --- prédiction : on ne fait qu'une itération.
  68       #          Il faut un RESU_THER (sinon on utilise la température du
  69       #          premier Gpcrit et cà n'a pas trop d'intéret).
  70       #          A chaque instant, on regarde à quelle température est le
  71       #          fond d'entaille et on compare Gpmax à cet instant au Gpcrit.
  72       mccalc = self['PREDICTION']
  73       l_crit = mccalc['GP_CRIT']
  74       lv_pred = []
  75       l_temp = mccalc['TEMP'][0]
  76
  77    if not type(l_temp) in EnumTypes:
  78       l_temp = [l_temp,]
  79    if not type(l_crit) in EnumTypes:
  80       l_crit = [l_crit,]
  81
  82
  83    # Maillage associe au modele
  84    __MAIL = aster.getvectjev( self['MODELE'].nom.ljust(8) + '.MODELE    .LGRF        ' )
  85    nom_maillage = __MAIL[0].strip()
  86
  87    maya = self.get_concept(nom_maillage)
  88
  89    # Excitation
  90    args={}
  91    if self['EXCIT']:args={'EXCIT'   : self['EXCIT'].List_F()}
  92
  93    # 1. ----- calcul de G-theta
  94
  95    # Cas 2D
  96    if self['THETA_2D'] is not None:
  97       is_2D = True
  98    else:
  99       is_2D = False
 100
 101    if is_2D:
 102       nbcour = len(self['THETA_2D'])
 103       nb_tranches = 1
 104       ep_z = 1
 105       l_tab = []
 106       for occ in self['THETA_2D']:
 107          dMC = occ.cree_dict_valeurs(occ.mc_liste)
 108
 109          __theta = CALC_THETA(MODELE=self['MODELE'],
 110                               DIRECTION=self['DIRECTION'],
 111                               THETA_2D=_F(GROUP_NO=dMC['GROUP_NO'],
 112                                           MODULE=1.0,
 113                                           R_INF=dMC['R_INF'],
 114                                           R_SUP=dMC['R_SUP']),)
 115
 116          __gtheta = CALC_G(THETA=_F(THETA=__theta),
 117                            RESULTAT=self['RESULTAT'],
 118                            TOUT_ORDRE='OUI',
 119                            SYME_CHAR=self['SYME_CHAR'],
 120                            COMP_ELAS=self['COMP_ELAS'].List_F(),
 121                            **args)
 122
 123          tab = __gtheta.EXTR_TABLE()
 124
 125          # une Table par couronne
 126          l_tab.append(tab)
 127
 128    else:
 129       #Cas 3D
 130       nbcour = len(self['THETA_3D'])
 131       nb_tranches = self['NB_TRANCHES']
 132       l_tab = []
 133       l_noeuds_fissure, pas = getFondFissInfo(self['FOND_FISS'])
 134       nb_noeuds_fissure = len(l_noeuds_fissure)
 135
 136       for occ in self['THETA_3D']:
 137          dMC = occ.cree_dict_valeurs(occ.mc_liste)
 138
 139          # on met les mots-clés facultatifs dans des dictionnaires
 140          dpar_theta = {}
 141          if self['DIRECTION'] is not None:
 142             dpar_theta['DIRECTION'] = self['DIRECTION']
 143
 144          __gtheta = CALC_G(
 145                            THETA=_F(R_INF=dMC['R_INF'],
 146                                     R_SUP=dMC['R_SUP'],
 147                                     MODULE=1.0,
 148                                     FOND_FISS=self['FOND_FISS'],
 149                                     **dpar_theta),
 150                            RESULTAT=self['RESULTAT'],
 151                            TOUT_ORDRE='OUI',
 152                            SYME_CHAR=self['SYME_CHAR'],
 153                            COMP_ELAS=self['COMP_ELAS'].List_F(),
 154                            LISSAGE=self['LISSAGE'].List_F(),
 155                            **args)
 156
 157
 158          tab = __gtheta.EXTR_TABLE()
 159
 160
 161          # une Table par couronne
 162          l_tab.append(tab)
 163
 164
 165
 166    # 2. ----- Calcul de l'energie élastique en exploitant les groupes de
 167    #          mailles fournis par la procedure de maillage
 168
 169    l_copo_tot = [grma.strip() for grma in self['GROUP_MA']]
 170    nbcop_tot = len(l_copo_tot)
 171    nbcop = nbcop_tot/nb_tranches
 172
 173    if self['LIST_EP_COPEAUX'] is not None:
 174       l_ep_copeaux_tot = self['LIST_EP_COPEAUX']
 175
 176    l_t_enel = []
 177
 178    if self['TRAC_COMP']=='OUI':
 179       # prise en compte de la traction-compression dans le calcul de l'energie
 180       resu2=CALC_ELEM(OPTION=('EQUI_ELNO_SIGM'),
 181                      RESULTAT=self['RESULTAT'],
 182                      )
 183
 184       # indices des mailles du dernier group_ma
 185       # (pour avoir le nombre de mailles par tranche)
 186       l_mailles_last_gm = maya.GROUPEMA.get()[l_copo_tot[-1].ljust(8)]
 187
 188       # initialisation des concepts reutilises dans la boucle
 189       # on suppose que chaque tranche a le meme nombre de mailles
 190
 191       kk = 0
 192
 193       E_el = [None]*len(l_mailles_last_gm)*nb_tranches
 194
 195       T_el = [None]*len(l_mailles_last_gm)*nb_tranches
 196
 197       # on recupere les sd en dehors de la boucle
 198       maya_GROUPEMA = maya.GROUPEMA.get()
 199       maya_NOMMAI = maya.NOMMAI.get()
 200       maya_CONNEX = maya.CONNEX.get()
 201       maya_NOMNOE = maya.NOMNOE.get()
 202
 203    # liste des tables tb_Gpmax repartie aux noeuds
 204    l_tb_Gpmax_noeuds = []
 205
 206    # Charges
 207    args={}
 208    if self['EXCIT']:args={'CHARGE': [charg['CHARGE'] for charg in self['EXCIT']]}
 209    for i in range(0,nb_tranches):
 210       l_copo = l_copo_tot[i*nbcop:(i+1)*nbcop]
 211
 212       if info >= 2 and not is_2D:
 213          print "<I> Calcul de la tranche %i"%(i+1)
 214
 215       if self['TRAC_COMP']=='OUI':
 216
 217          # l_copo est une liste commulative de mailles
 218          # il faut lancer POST_ELEM sur chaque maille de chaque copeau
 219          # puis regarder la trace de SIEF_ELGA sur ce copeau
 220
 221          # on fera attention a ne pas lancer POST_ELEM sur une maille qui
 222          # a deja ete calculee en stockant son resultat pour la maille en question
 223          d_groupma={}
 224          d_nomma={}
 225
 226          # indices des mailles des copeaux
 227          for group_ma in l_copo:
 228             d_groupma[group_ma] = maya_GROUPEMA[group_ma.ljust(8)]
 229
 230          # le dernier copeau contient tous les elements
 231          # on calcule l energie de chaque element de ce copeau
 232          last_copo = l_copo[-1]
 233
 234          d_ener = {}
 235
 236          d_nomma = {}
 237
 238          for k, id_elem in enumerate(d_groupma[last_copo]):
 239
 240
 241             # les id des elements dans Aster commencent a 1
 242             # la liste python commence a 0
 243             elem = maya_NOMMAI[id_elem-1]
 244             d_nomma[id_elem]=elem
 245
 246             E_el[kk] = POST_ELEM(MODELE=self['MODELE'],
 247                                  RESULTAT=self['RESULTAT'],
 248                                  TOUT_ORDRE='OUI',
 249                                  ENER_ELAS=_F(MAILLE=elem),
 250                                  TITRE='Energie élastique',
 251                                  **args)
 252
 253             T_el[kk] = E_el[kk].EXTR_TABLE()
 254
 255             l_enel = T_el[kk].TOTALE.values()
 256
 257             # signe de la trace <=> signe de la composante VMIS_SG du tenseur EQUI_ELNO_SIGM,
 258             # mais E_enel est par element => on fait une moyenne sur les noeuds de l'element
 259
 260             list_no = []
 261             for ind_no in maya_CONNEX[id_elem] :
 262                nomnoe = maya_NOMNOE[ind_no-1]
 263                if nomnoe not in list_no :
 264                   list_no.append(nomnoe)
 265
 266             # print "liste des noeuds de la maille ", id_elem, ": ", list_no
 267
 268             l_inst = T_el[kk].INST.values()
 269             nb_inst = len(l_inst)
 270
 271             T_noeuds = Table()
 272             T_noeuds['INST']=l_inst
 273
 274             # pour chaque noeud de l'element on recupere sa trace
 275             for noeud in list_no:
 276
 277                __VM=RECU_FONCTION(RESULTAT=resu2,
 278                                     TOUT_INST='OUI',
 279                                     NOM_CHAM='EQUI_ELNO_SIGM',
 280                                     NOM_CMP='VMIS_SG',
 281                                     MAILLE=elem,
 282                                     NOEUD=noeud);
 283
 284                T_noeuds[noeud]=__VM.Ordo()
 285
 286             T_noeuds.fromfunction('VM_MAIL', moyenne, list_no)
 287
 288             l_VM_MAIL = T_noeuds.VM_MAIL.values()
 289
 290             for j, vm in enumerate(l_VM_MAIL):
 291                if vm < 0:
 292                   l_enel[j]=-l_enel[j]
 293
 294             del T_el[kk]['TOTALE']
 295             T_el[kk][elem]=l_enel
 296
 297             if k==0:
 298
 299                # Table de l'energie elastique sur le GROUP_MA
 300                T_el_gm = Table()
 301                T_el_gm['NUME_ORDRE'] = T_el[kk].NUME_ORDRE.values()
 302                T_el_gm['INST'] = T_el[kk].INST.values()
 303                T_el_gm['LIEU'] = [last_copo]*nb_inst
 304                T_el_gm['ENTITE'] = ['GROUP_MA']*nb_inst
 305
 306             T_el_gm[elem]=l_enel
 307
 308             kk+=1
 309
 310          # sommation sur les mailles du group_ma:
 311          l_nomma = d_nomma.values()
 312          T_el_gm.fromfunction('TOTALE', mysum, l_nomma)
 313
 314          # Table totale
 315          t_enel=Table(titr="Energie élastique")
 316          t_enel['NUME_ORDRE']=T_el_gm.NUME_ORDRE.values()
 317          t_enel['INST']=T_el_gm.INST.values()
 318          t_enel['LIEU']=T_el_gm.LIEU.values()
 319          t_enel['ENTITE']=T_el_gm.ENTITE.values()
 320          t_enel['TOTALE']=T_el_gm.TOTALE.values()
 321
 322          # t_enel ne contient jusqu'ici que l'energie elastique du dernier copeau
 323
 324          # calcul de l'energie elastique pour les autres copeaux
 325          T_el_sub = T_el_gm.copy()
 326
 327          for k in range(len(l_copo)-2,-1,-1):
 328             group_ma = l_copo[k]
 329             T_el_sub = T_el_sub.copy()
 330             del T_el_sub['LIEU']
 331             del T_el_sub['TOTALE']
 332             T_el_sub['LIEU']=[group_ma]*nb_inst
 333             l_id_elem = d_groupma[group_ma]
 334             l_nom_elem = []
 335
 336             for id_elem, nom_elem in d_nomma.items():
 337                if not id_elem in l_id_elem:
 338                   # colonne a supprimer
 339                   del T_el_sub[nom_elem]
 340                   del d_nomma[id_elem]
 341                else:
 342                   l_nom_elem.append(nom_elem)
 343
 344             T_el_sub.fromfunction('TOTALE', sum_and_check, l_nom_elem)
 345
 346             # Table de l'energie elastique sur le GROUP_MA
 347             T_el_gm_k = Table()
 348             T_el_gm_k['NUME_ORDRE'] =T_el_sub.NUME_ORDRE.values()
 349             T_el_gm_k['INST'] = T_el_sub.INST.values()
 350             T_el_gm_k['LIEU'] = [group_ma]*nb_inst
 351             T_el_gm_k['ENTITE'] = ['GROUP_MA']*nb_inst
 352             T_el_gm_k['TOTALE'] = T_el_sub.TOTALE.values()
 353
 354             # contribution du group_ma a la table totale:
 355             t_enel = merge(t_enel, T_el_gm_k)
 356
 357          t_enel.sort('NUME_ORDRE')
 358
 359       else:
 360          # si self['TRAC_COMP']!='OUI'
 361          # calcul classique de l'energie elastique
 362
 363          __ener = POST_ELEM(MODELE=self['MODELE'],
 364                                  RESULTAT=self['RESULTAT'],
 365                                  TOUT_ORDRE='OUI',
 366                                  ENER_ELAS=_F(GROUP_MA=l_copo),
 367                                  TITRE='Energie élastique',
 368                                  **args)
 369
 370          t_enel = __ener.EXTR_TABLE()
 371
 372       # 2.1. ----- Indice de chaque copeau et deltaL
 373       d_icop = dict(zip(l_copo, range(1, nbcop + 1)))
 374
 375       l_lieu = [grma.strip() for grma in t_enel.LIEU.values()]
 376
 377       l_icop = [d_icop[grma] for grma in l_lieu]
 378
 379       t_enel['ICOP'] = l_icop
 380
 381       l_numord = list(set(t_enel.NUME_ORDRE.values()))
 382       l_numord.sort()
 383
 384       if self['PAS_ENTAILLE'] is not None:
 385          t_enel.fromfunction('DELTAL', fDL, 'ICOP', { 'pascop' : self['PAS_ENTAILLE'] })
 386       else:
 387          l_ep_copeaux_tranche = l_ep_copeaux_tot[i*nbcop:(i+1)*nbcop]
 388          t_enel['DELTAL'] = l_ep_copeaux_tranche*len(l_numord)
 389
 390       # 2.2. ----- Calcul de Gp fonction de Ener.Totale et de deltaL
 391       if is_2D:
 392          t_enel.fromfunction('GP', fGp_Etot, ('TOTALE', 'DELTAL'),
 393             {'syme'   : self['SYME_CHAR'] != 'SANS',
 394              'R'      : self['RAYON_AXIS'],})
 395       else:
 396          ep_tranche = largeur_tranche(nom_maillage, l_noeuds_fissure, pas, i)
 397          #print "ep_tranche %i: "%i, ep_tranche
 398          t_enel.fromfunction('GP', fGp_Etot, ('TOTALE', 'DELTAL'),
 399                {'syme'   : self['SYME_CHAR'] != 'SANS',
 400                 'R'      : ep_tranche })
 401
 402       #if info >= 2:
 403       #   print "Table de l'énergie élastique: ", t_enel
 404
 405       l_t_enel.append(t_enel)
 406       # 2.3. ----- Tableau de Gp = f(icop) pour chaque instant
 407       if info >= 2:
 408          tGp_t_icop = t_enel['INST', 'DELTAL', 'GP']
 409          tGp_t_icop.titr = "Gp à chaque instant en fonction de la distance au " \
 410                            "fond d'entaille"
 411          tGp_t_icop.ImprTabCroise()
 412
 413       # 2.4. ----- Table Gpmax
 414       ttmp = t_enel['NUME_ORDRE', 'INST', 'ICOP', 'DELTAL', 'GP']
 415
 416       for j in l_numord:
 417          tj = ttmp.NUME_ORDRE == j
 418
 419          ## pour tester le comportement de Gpmax quand GP est identiquement nul
 420          #del tj['GP']
 421          #tj['GP']=[0]*len(tj.GP.values())
 422
 423          if self['CRIT_MAXI_GP'] == 'ABSOLU':
 424             t = tj.GP.MAXI()
 425          else:
 426             t = MaxRelatif(tj, 'GP')
 427
 428
 429          # cas GP identiquement nul: plusieurs max de GP
 430          # on prend le DELTAL minimum
 431          if len(t.GP.values())>1:
 432             t = t.DELTAL.MINI()
 433
 434          if j == 1:
 435             tb_Gpmax_i = t
 436          else:
 437             tb_Gpmax_i = tb_Gpmax_i | t
 438
 439
 440       tb_Gpmax_i.Renomme('GP', 'GPMAX')
 441       tb_Gpmax_i.Renomme('ICOP', 'ICOPMAX')
 442       tb_Gpmax_i.Renomme('DELTAL', 'DELTALMAX')
 443       tb_Gpmax_i.titr = 'Gpmax à chaque instant'
 444
 445       # On transfert Gpmax au noeud sommet à gauche et au milieu (si cas quadratique)
 446       # sauf pour la dernière tranche où on transfère également au noeud sommet de droite.
 447       # Tout cela pour pouvoir avoir une table complète avec les G et les Gpmax.
 448       if not is_2D:
 449          tb_Gpmax_i['NUME_TRANCHE']=[i+1]*len(tb_Gpmax_i['GPMAX'])
 450          if i==0:
 451             l_inst = tb_Gpmax_i.INST.values()
 452             nb_inst = len(l_inst)
 453          if pas==1:
 454             tb_Gpmax_i_noeuds = tb_Gpmax_i.copy()
 455             tb_Gpmax_i_noeuds['NOEUD']=[l_noeuds_fissure[i]]*nb_inst
 456             l_tb_Gpmax_noeuds.append(tb_Gpmax_i_noeuds)
 457          else:
 458             tb_Gpmax_i_noeuds_1 = tb_Gpmax_i.copy()
 459             tb_Gpmax_i_noeuds_1['NOEUD'] = [l_noeuds_fissure[pas*i]]*nb_inst
 460             l_tb_Gpmax_noeuds.append(tb_Gpmax_i_noeuds_1)
 461             tb_Gpmax_i_noeuds_2 = tb_Gpmax_i.copy()
 462             tb_Gpmax_i_noeuds_2['NOEUD'] = [l_noeuds_fissure[pas*i+1]]*nb_inst
 463             l_tb_Gpmax_noeuds.append(tb_Gpmax_i_noeuds_2)
 464          if i==nb_tranches-1:
 465             tb_Gpmax_i_noeuds_3 = tb_Gpmax_i.copy()
 466             tb_Gpmax_i_noeuds_3['NOEUD'] = [l_noeuds_fissure[-1]]*nb_inst
 467             l_tb_Gpmax_noeuds.append(tb_Gpmax_i_noeuds_3)
 468
 469
 470
 471       if i == 0:
 472          tb_Gpmax = tb_Gpmax_i
 473       else:
 474          tb_Gpmax = merge(tb_Gpmax, tb_Gpmax_i)
 475
 476    # FIN BOUCLE SUR LES TRANCHES
 477
 478    if not is_2D:
 479       tb_Gpmax_noeuds = Table(para=tb_Gpmax.para+['NOEUD'])
 480       for j, tb in enumerate(l_tb_Gpmax_noeuds):
 481          if j==0:
 482             tb_Gpmax_noeuds = tb
 483          else:
 484             tb_Gpmax_noeuds = merge(tb_Gpmax_noeuds, tb)
 485
 486
 487    # 2.5. ----- extraction de la température en fond d'entaille
 488    #voir le cas 3D => THETA_3D, mais qu'en est-il des tranches?
 489    if self['RESU_THER']:
 490       #sur un seul noeud ou sur tous les noeuds du fond d'entaille?
 491
 492       if is_2D:
 493          grno_fond = self['THETA_2D'][0]['GROUP_NO']
 494       else:
 495          grma_fond = self['THETA_3D'][0]['GROUP_MA']
 496          grno_fond = "GRNOFOND"
 497          DEFI_GROUP(reuse =maya,
 498                     MAILLAGE=maya,
 499                     CREA_GROUP_NO=_F(GROUP_MA=grma_fond,
 500                                      NOM=grno_fond,),);
 501
 502       l_ordres = DEFI_LIST_ENTI(VALE=l_numord)
 503       __relev = POST_RELEVE_T(ACTION=_F(RESULTAT=self['RESU_THER'],
 504                                         OPERATION='EXTRACTION',
 505                                         INTITULE='Temperature',
 506                                         NOM_CHAM='TEMP',
 507                                         LIST_ORDRE=l_ordres,
 508                                         NOM_CMP='TEMP',
 509                                         GROUP_NO=grno_fond,),)
 510
 511       t_relev = __relev.EXTR_TABLE()['NUME_ORDRE', 'NOEUD', 'TEMP']
 512
 513    # 3. ----- boucle sur les mots-clés facteurs
 514    #          opérations dépendant de la température
 515    MATER = self['MATER']
 516    flag_mat = True
 517
 518    for iocc, TEMP in enumerate(l_temp):
 519       # 3.0. ----- Temperature fonction du temps : si on n'a pas de RESU_THER,
 520       #            on prend la température d'identification.
 521       if not self['RESU_THER']:
 522          l_rows = [{'NUME_ORDRE' : j, 'TEMP' : TEMP} for j in l_numord]
 523          t_relev = Table(rows=l_rows, para=('NUME_ORDRE', 'TEMP'), typ=('R', 'R'))
 524          flag_mat = True
 525
 526       # 3.1. ----- extrait du matériau E(TEMP) et NU(TEMP) (si nécessaire)
 527       if flag_mat:
 528          t_relev.fromfunction('YOUNG', CallRCVALE, 'TEMP',
 529                { 'para' : 'E', 'MATER' : MATER })
 530          t_relev.fromfunction('NU', CallRCVALE, 'TEMP',
 531                { 'para' : 'NU', 'MATER' : MATER })
 532          flag_mat = False
 533
 534       # 3.2. ----- paramètres
 535       dict_constantes = {
 536          'YOUNG' : CallRCVALE(TEMP, 'E', MATER),
 537          'NU'    : CallRCVALE(TEMP, 'NU', MATER),
 538       }
 539
 540       # 3.3. ----- calcul de Kj(G)
 541       l_tabi = []
 542       for k, tab in enumerate(l_tab):
 543          #tab: table de la couronne k
 544
 545          # calcul de Kj(G) = K_i
 546          new_para = 'K_%d' % (k + 1)
 547          if is_2D:
 548             # fusion avec TEMP, E et nu
 549             tab = merge(tab, t_relev, 'NUME_ORDRE')
 550             tab.fromfunction(new_para, fKj, ('G', 'YOUNG', 'NU'))
 551             # renomme G en G_i
 552             tab.Renomme('G', 'G_%d' % (k + 1))
 553          else:
 554             if self['RESU_THER']:
 555                tab=merge(tab, t_relev, ['NUME_ORDRE', 'NOEUD'])
 556             else:
 557                tab=mergeLineInTable(tab, t_relev, nb_noeuds_fissure)
 558
 559             # en 3D, le paramètre R n'intervient pas
 560             tab.fromfunction(new_para, fKj, ('G_LOCAL', 'YOUNG', 'NU'))
 561             tab.Renomme('G_LOCAL', 'G_%d' % (k + 1))
 562
 563          l_tabi.append(tab)
 564
 565       # 3.4 ----- Table des Gi, Ki sur les differentes couronnes + Kmoyen
 566       if is_2D:
 567          tabK_G = l_tabi[0]['NUME_ORDRE']
 568          for tab in l_tabi:
 569             tabK_G = merge(tabK_G, tab, 'NUME_ORDRE')
 570       else:
 571          tabK_G=l_tabi[0]
 572          for i in range(1,len(l_tabi)):
 573             tabK_G = merge(tabK_G, l_tabi[i], ['NUME_ORDRE', 'NOEUD'])
 574       tabK_G.titr = 'G et K sur les differentes couronnes + moyennes'
 575       tabK_G.fromfunction('GMOY', moyenne_positive, ['G_%d' % (k + 1) for k in range(nbcour)])
 576       tabK_G.fromfunction('KMOY', moyenne, ['K_%d' % (k + 1) for k in range(nbcour)])
 577
 578       # 3.5. ----- Contribution à la table globale
 579
 580       if is_2D:
 581          tabres = merge(tabK_G, tb_Gpmax, 'NUME_ORDRE')
 582          tabres['OCCURRENCE'] = [iocc + 1] * len(l_numord)
 583       else:
 584          # tb_Gpmax est une table sur les tranches,
 585          # on l'ajoute dans la table aux noeuds tabres avec la convention:
 586          # au 1er noeud et noeud milieu de la tranche on affecte la valeur de la tranche
 587          # sauf pour la derniere tranche où on affecte la valeur sur les 3 noeuds de la tranche
 588          tabres = tabK_G
 589          tabres = merge(tabK_G, tb_Gpmax_noeuds, ['NUME_ORDRE', 'NOEUD'])
 590          tabres['OCCURRENCE'] = [iocc + 1] * len(l_numord) * nb_noeuds_fissure
 591       #if info >= 2:
 592       #   print tabres
 593
 594       # 3.5.1. --- Table globale
 595       if iocc == 0:
 596          tabl_glob = tabres
 597       else:
 598          tabl_glob = merge(tabl_glob, tabres)
 599       tabl_glob.titr = 'G, K sur les differentes couronnes, Gmoy, Kmoy et ' \
 600                        'Gpmax fonctions du temps'
 601
 602       # 3.6. ----- traitement selon identification / prédiction
 603       d_para = {
 604          'INTERPOL' : ['LIN', 'LIN'],
 605          'NOM_PARA' : 'INST',
 606          'PROL_DROITE' : 'CONSTANT',
 607          'PROL_GAUCHE' : 'CONSTANT',
 608       }
 609
 610       # 3.6.1. --- identification
 611       if identification:
 612          KJ_CRIT = l_crit[iocc]
 613          # on verifie que KJ_CRIT soit compris dans l'intervalle [KMOY_min, KMOY_max]
 614          valkmoy = tabres.KMOY.values()
 615          if not (min(valkmoy) <= KJ_CRIT <= max(valkmoy)):
 616 #                               'constant utilisé).')
 617             UTMESS('A','RUPTURE0_1')
 618          if is_2D:
 619             # définition des fonctions pour faire les interpolations
 620             d_para.update({ 'NOM_RESU' : 'DELTALMAX' })
 621             # DeltaMax en fonction du temps
 622             fdL = t_fonction(tb_Gpmax.INST.values(), tb_Gpmax.DELTALMAX.values(), d_para)
 623             # Gpmax fonction du temps
 624             d_para.update({ 'NOM_RESU' : 'GPMAX' })
 625             fGp = t_fonction(tb_Gpmax.INST.values(), tb_Gpmax.GPMAX.values(), d_para)
 626
 627             d_para.update({ 'NOM_PARA' : 'KMOY',
 628                            'NOM_RESU' : 'INST', })
 629             valkmoy = tabres.KMOY.values()
 630             # temps en fonction de Kmoy
 631             finv = t_fonction(valkmoy, tabres.INST.values(), d_para)
 632
 633             # valeurs à mettre dans la table
 634             # temps correspondant à KJ_CRIT
 635             ti   = finv(KJ_CRIT)
 636             # GP correspondant au temps critique
 637             Gpi  = fGp(ti)
 638             d_ident = {
 639                'KJ_CRIT'   : KJ_CRIT,
 640                'INST'      : ti,
 641                'GPMAX'     : Gpi,
 642                'KGPMAX'    : fKj(Gpi, **dict_constantes),
 643                'DELTALMAX' : fdL(ti),
 644             }
 645             lv_ident.append(d_ident)
 646          else:
 647
 648             l_i_noeuds_sommets = range(0,len(l_noeuds_fissure),pas)
 649             t_noeud_Kcrit = Table(para=tabres.para)
 650
 651             # On determine quel noeud sommet maximise KMOY au cours du temps:
 652
 653             row_KMOY_max = tabres.KMOY.MAXI()
 654             noeud_KMOY_max = row_KMOY_max.NOEUD.values()[0]
 655
 656             # avec le noeud ou KJ_CRIT est atteint, on regarde GP a gauche et a droite.
 657             # le GP le plus grand correspond au GPmax
 658             # on garde la tranche ou GP est le plus grand
 659
 660             id_noeud_KMOY_max = list(l_noeuds_fissure).index(noeud_KMOY_max)
 661             if id_noeud_KMOY_max==0:
 662                # "Gpmax sur la 1ere tranche"
 663                nume_tranche_Gpmax = 1
 664             elif id_noeud_KMOY_max==(len(l_noeuds_fissure)-1):
 665                # "Gpmax sur la derniere tranche"
 666                nume_tranche_Gpmax = nb_tranches
 667             else:
 668                # "Gpmax sur une tranche intermediaire"
 669                Gpi_tot = Table(para=tb_Gpmax.para)
 670                Gpi_gauche = tb_Gpmax.NUME_TRANCHE==(id_noeud_KMOY_max/pas)
 671                Gpi_tot.append(Gpi_gauche.rows[0])
 672                Gpi_droite = tb_Gpmax.NUME_TRANCHE==(id_noeud_KMOY_max/pas+1)
 673                Gpi_tot.append(Gpi_droite.rows[0])
 674                Gpi_tab = Gpi_tot.GPMAX.MAXI()
 675                nume_tranche_Gpmax = Gpi_tab.NUME_TRANCHE.values()[0]
 676
 677             # tb_Gpmax_TrancheCrit est une table de la meme nature que la table 2D tb_Gpmax
 678             # i.e. valeurs sur la seule tranche qui nous interesse (ou on sait
 679             # que KJ_CRIT sera atteint)
 680
 681             tb_Gpmax_TrancheCrit = tb_Gpmax.NUME_TRANCHE==nume_tranche_Gpmax
 682
 683             # avec le noeud ou KJ_CRIT est atteint, on determine le temps
 684             # critique par interpolation
 685             tabres_NoeudCrit = tabres.NOEUD==noeud_KMOY_max
 686
 687             # la suite est idem 2D, seuls les noms des tables changent
 688
 689             # définition des fonctions pour faire les interpolations
 690             d_para.update({ 'NOM_RESU' : 'DELTALMAX' })
 691             # DeltaMax en fonction du temps
 692             fdL = t_fonction(tb_Gpmax_TrancheCrit.INST.values(),
 693                               tb_Gpmax_TrancheCrit.DELTALMAX.values(),
 694                               d_para)
 695
 696             # Gpmax fonction du temps
 697             d_para.update({ 'NOM_RESU' : 'GPMAX' })
 698             fGp = t_fonction(tb_Gpmax_TrancheCrit.INST.values(),
 699                               tb_Gpmax_TrancheCrit.GPMAX.values(),
 700                               d_para)
 701
 702
 703             d_para.update({ 'NOM_PARA' : 'KMOY',
 704                            'NOM_RESU' : 'INST', })
 705             valkmoy = tabres_NoeudCrit.KMOY.values()
 706             # temps en fonction de Kmoy
 707             finv = t_fonction(valkmoy, tabres_NoeudCrit.INST.values(), d_para)
 708
 709             # valeurs à mettre dans la table
 710             # temps correspondant a KJ_CRIT
 711             ti   = finv(KJ_CRIT)
 712             # GP correspondant au temps critique
 713             Gpi  = fGp(ti)
 714             # par rapport a 2D, on ajoute 'NUME_TRANCHE'
 715             d_ident = {
 716                'KJ_CRIT'      : KJ_CRIT,
 717                'INST'         : ti,
 718                'NUME_TRANCHE' : int(nume_tranche_Gpmax),
 719                'GPMAX'        : Gpi,
 720                'KGPMAX'       : fKj(Gpi, **dict_constantes),
 721                'DELTALMAX'    : fdL(ti),
 722             }
 723             lv_ident.append(d_ident)
 724
 725       # 3.6.2. --- prédiction
 726       else:
 727          pass
 728
 729    # 4. ----- Construction de la table résultat si demandée
 730    # 4.1. --- identification
 731    if identification:
 732       if is_2D:
 733          para_tab_ident=('KJ_CRIT', 'INST', 'GPMAX', 'KGPMAX', 'DELTALMAX')
 734       else:
 735          para_tab_ident=('KJ_CRIT', 'INST', 'NUME_TRANCHE', 'GPMAX', 'KGPMAX', 'DELTALMAX')
 736       tab_ident = Table(rows=lv_ident,
 737                         para=para_tab_ident,
 738                         typ = ('R')*len(para_tab_ident),
 739                         titr='Identification aux valeurs de tenacités critiques')
 740       dprod_result = tab_ident.dict_CREA_TABLE()
 741       if info >= 2:
 742          print tab_ident
 743
 744    # 4.2. --- prédiction
 745    else:
 746       # définition de la fonction GPcrit = f(TEMP)
 747       d_para.update({ 'NOM_PARA' : 'TEMP',
 748                       'NOM_RESU' : 'GP_CRIT', })
 749       fGpc = t_fonction(mccalc['TEMP'], mccalc['GP_CRIT'], d_para)
 750
 751       # en 3D, GPMAX et DELTALMAX ne sont pas dans tab_glob: les recuperer de tb_Gpmax
 752
 753       if is_2D:
 754          tab_pred = tabl_glob['NUME_ORDRE', 'INST', 'TEMP', 'DELTALMAX', 'GPMAX']
 755       else:
 756          tab_pred = tb_Gpmax['NUME_ORDRE', 'INST', 'NUME_TRANCHE', 'DELTALMAX', 'GPMAX']
 757          # on recupere TEMP de tabl_glob
 758          tab_temp_tranche=Table(para=['NUME_ORDRE', 'NUME_TRANCHE', 'TEMP'])
 759          # on fait une moyenne de la temperature sur les noeuds d'une tranche
 760          for ordre in l_numord:
 761             tabl_glob_ORDRE_i = tabl_glob.NUME_ORDRE==ordre
 762             temp_ORDRE_i = tabl_glob_ORDRE_i.TEMP.values()
 763             for i_tranche in range(nb_tranches):
 764                l_temp_noeuds_tranche_i = temp_ORDRE_i[i_tranche : i_tranche+pas+1]
 765                temp_tranche_i = moyenne(*l_temp_noeuds_tranche_i)
 766                d = {'NUME_ORDRE': ordre, 'NUME_TRANCHE': i_tranche+1, 'TEMP': temp_tranche_i}
 767                tab_temp_tranche.append(d)
 768          tab_pred = merge(tab_pred, tab_temp_tranche, ['NUME_ORDRE', 'NUME_TRANCHE'])
 769
 770       tab_pred.fromfunction('GP_CRIT', fGpc, 'TEMP')
 771       tab_pred.fromfunction('PREDICTION', crit, ('GP_CRIT', 'GPMAX'))
 772       tab_pred.titr = 'Comparaison Gpmax à Gpcrit(T)'
 773       dprod_result = tab_pred.dict_CREA_TABLE()
 774
 775    # 9. ----- création de la table_sdaster résultat
 776    dprod = tabl_glob.dict_CREA_TABLE()
 777    result = CREA_TABLE(**dprod)
 778    tabresult = CREA_TABLE(**dprod_result)
 779
 780
 781
 782 # -----------------------------------------------------------------------------
 783 def CallRCVALE(TEMP, para, MATER):
 784    """Fonction appelant RCVALE et retourne la valeur d'un paramètre.
 785    """
 786    valres, flag_ok = MATER.RCVALE('ELAS', 'TEMP', TEMP, para)
 787    assert list(flag_ok).count('OK') != 1, \
 788          'Erreur lors de la récupération des valeurs du matériau.'
 789    return valres
 790
 791 # -----------------------------------------------------------------------------
 792 def fKj(G, YOUNG, NU):
 793    """Calcul de Kj à partir de G (formule d'Irwin)
 794    """
 795    Kj=(abs(G * YOUNG / (1.0 - NU**2)))**0.5
 796    return Kj
 797
 798 # -----------------------------------------------------------------------------
 799 def fDL(ICOP, pascop):
 800    """DeltaL = numéro copeau * pas d'entaille
 801    """
 802    return ICOP * pascop
 803
 804 # -----------------------------------------------------------------------------
 805 def fGp_Etot(TOTALE, DELTAL, R, syme=False):
 806    """Gp(Etotale, K), deltal pris dans le context global.
 807       ICOP      : numéro du copeau,
 808       DELTAL    : liste des epaisseurs des copeaux
 809       R         : rayon en axisymetrique,
 810                   longueur de l'élément 1D situé sur le front d'entaille si modèle 3D.
 811       syme      : True s'il y a symétrie.
 812    """
 813    import types
 814    fact_syme = 1.
 815    if syme:
 816       fact_syme = 2.
 817    Gp_Etot = fact_syme * TOTALE / (DELTAL * R )
 818    return Gp_Etot
 819
 820 # -----------------------------------------------------------------------------
 821 def MaxRelatif(table, nom_para):
 822    """Extrait le dernier maxi du champ `nom_para` de la table.
 823    """
 824    l_val = getattr(table, nom_para).values()
 825    l_val.reverse()
 826    Vlast = l_val[0]
 827    for val in l_val:
 828       if val < Vlast:
 829          break
 830       Vlast = val
 831    return getattr(table, nom_para) == Vlast
 832
 833 # -----------------------------------------------------------------------------
 834 def crit(GP_CRIT, GPMAX):
 835    """Retourne 1 quand GP_CRIT > GPMAX
 836    """
 837    if GPMAX > GP_CRIT:
 838       return 1
 839    else:
 840       return 0
 841
 842 # -----------------------------------------------------------------------------
 843 def moyenne(*args):
 844    """Fonction moyenne
 845    """
 846    return sum(args)/len(args)
 847
 848 def moyenne_positive(*args):
 849    """Fonction moyenne
 850    """
 851    return sum([abs(a) for a in args])/len(args)
 852
 853 def mysum(*args):
 854    """Fonction sum.
 855       La fonction sum ne peut pas etre appelee sur une liste de parametre
 856       d'une table via fromfunction
 857    """
 858    return sum(args)
 859
 860 def sum_and_check(*args):
 861    """Fonction sum.
 862       Verifie si la somme est positive.
 863       Si la somme est negative, on la met egale a zero
 864    """
 865    somme = sum(args)
 866    if somme<0:
 867       somme=0
 868    return somme
 869
 870 # On recupere des infos sur le fond de fissure
 871 def getFondFissInfo(fondfiss):
 872    # >FONFISS .FOND      .NOEU        <
 873    # >FONFISS .FOND      .TYPE        <
 874    import aster
 875    l_noeuds_fissure = aster.getvectjev(fondfiss.nom.ljust(8)+'.FOND      .NOEU        ')
 876    type_mailles = aster.getvectjev(fondfiss.nom.ljust(8)+'.FOND      .TYPE        ')
 877    if (type_mailles[0].strip() == 'SEG3' ):
 878       pas = 2
 879    else:
 880       pas = 1
 881    return l_noeuds_fissure, pas
 882
 883 ########################################################################
 884 # determination de la distance min entre 2 points consécutifs de la ligne de coupe
 885
 886 def largeur_tranche(nom_maillage, l_noms_noeuds_fissure, pas, i_tranche):
 887    # >MA      .COORDO    .VALE        <
 888    from math import sqrt
 889    import aster
 890
 891    # tuple des noms des noeuds du maillage
 892    t_noms_noeuds_maillage = aster.getvectjev(nom_maillage.ljust(8)+'.NOMNOE')
 893    # on convertit en liste pour utiliser la methode index
 894    # qui est plus optimal qu'une boucle sur les indices du tuple
 895    l_noms_noeuds_maillage = list(t_noms_noeuds_maillage)
 896
 897    l_numeros_noeuds_fissure = []
 898    for i in range(0,len(l_noms_noeuds_fissure),pas):
 899       nom = l_noms_noeuds_fissure[i]
 900       index = l_noms_noeuds_maillage.index(nom)
 901       l_numeros_noeuds_fissure.append(index)
 902
 903    coor1=aster.getvectjev(nom_maillage.ljust(8)+'.COORDO    .VALE        ',
 904                         l_numeros_noeuds_fissure[i_tranche]*3,3)
 905    coor2=aster.getvectjev(nom_maillage.ljust(8)+'.COORDO    .VALE        ',
 906                         l_numeros_noeuds_fissure[i_tranche+1]*3,3)
 907
 908    d=sqrt( (coor1[0]-coor2[0])**2+(coor1[1]-coor2[1])**2+(coor1[2]-coor2[2])**2)
 909    return d
 910
 911 def mergeLineInTable(multiTable, lineTable, nb_noeuds):
 912    # on ajoute a la table multiTable les colonnes de lineTable
 913    # pour chaque nume_ordre autant de fois qu'il y a de nb_noeuds
 914    from Utilitai.Table      import Table, merge
 915
 916    l_ordre = lineTable.NUME_ORDRE
 917    l_para = lineTable.copy().para
 918    l_para.remove('NUME_ORDRE')
 919    for i, ordre in enumerate(l_ordre):
 920       multiTable_i = multiTable.NUME_ORDRE==ordre
 921       row_i = lineTable.rows[i]
 922       for para in l_para:
 923          valeur_i = row_i[para]
 924          multiTable_i[para] = [valeur_i] * nb_noeuds
 925       if i==0:
 926          newTable=multiTable_i
 927       else:
 928          newTable = merge(newTable, multiTable_i)
 929
 930    return newTable