Aster/Cata/cataSTA9/Macro/post_gp_ops.py

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