Aster/Cata/cataSTA9/Macro/post_gp_ops.py

   1 #@ MODIF post_gp_ops Macro  DATE 18/11/2009   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
 542       # 3.3. ----- calcul de Kj(G)
 543       l_tabi = []
 544       for k, tab in enumerate(l_tab):
 545          #tab: table de la couronne k
 546
 547          # calcul de Kj(G) = K_i
 548          new_para = 'K_%d' % (k + 1)
 549          if is_2D:
 550             # fusion avec TEMP, E et nu
 551             tab = merge(tab, t_relev, 'NUME_ORDRE')
 552             tab.fromfunction(new_para, fKj, ('G', 'YOUNG', 'NU'))
 553             # renomme G en G_i
 554             tab.Renomme('G', 'G_%d' % (k + 1))
 555          else:
 556             if self['RESU_THER']:
 557                tab=merge(tab, t_relev, ['NUME_ORDRE', 'NOEUD'])
 558             else:
 559                tab=mergeLineInTable(tab, t_relev, nb_noeuds_fissure)
 560
 561             # en 3D, le paramètre R n'intervient pas
 562             tab.fromfunction(new_para, fKj, ('G_LOCAL', 'YOUNG', 'NU'))
 563             tab.Renomme('G_LOCAL', 'G_%d' % (k + 1))
 564
 565          l_tabi.append(tab)
 566
 567       # 3.4 ----- Table des Gi, Ki sur les differentes couronnes + Kmoyen
 568       if is_2D:
 569          tabK_G = l_tabi[0]['NUME_ORDRE']
 570          for tab in l_tabi:
 571             tabK_G = merge(tabK_G, tab, 'NUME_ORDRE')
 572       else:
 573          tabK_G=l_tabi[0]
 574          for i in range(1,len(l_tabi)):
 575             tabK_G = merge(tabK_G, l_tabi[i], ['NUME_ORDRE', 'NOEUD'])
 576       tabK_G.titr = 'G et K sur les differentes couronnes + moyennes'
 577       tabK_G.fromfunction('GMOY', moyenne_positive, ['G_%d' % (k + 1) for k in range(nbcour)])
 578       tabK_G.fromfunction('KMOY', moyenne, ['K_%d' % (k + 1) for k in range(nbcour)])
 579
 580       # 3.5. ----- Contribution à la table globale
 581
 582       if is_2D:
 583          tabres = merge(tabK_G, tb_Gpmax, 'NUME_ORDRE')
 584          tabres['OCCURRENCE'] = [iocc + 1] * len(l_numord)
 585       else:
 586          # tb_Gpmax est une table sur les tranches,
 587          # on l'ajoute dans la table aux noeuds tabres avec la convention:
 588          # au 1er noeud et noeud milieu de la tranche on affecte la valeur de la tranche
 589          # sauf pour la derniere tranche où on affecte la valeur sur les 3 noeuds de la tranche
 590          tabres = tabK_G
 591          tabres = merge(tabK_G, tb_Gpmax_noeuds, ['NUME_ORDRE', 'NOEUD'])
 592          tabres['OCCURRENCE'] = [iocc + 1] * len(l_numord) * nb_noeuds_fissure
 593       #if info >= 2:
 594       #   print tabres
 595
 596       # 3.5.1. --- Table globale
 597       if iocc == 0:
 598          tabl_glob = tabres
 599       else:
 600          tabl_glob = merge(tabl_glob, tabres)
 601       tabl_glob.titr = 'G, K sur les differentes couronnes, Gmoy, Kmoy et ' \
 602                        'Gpmax fonctions du temps'
 603
 604       # 3.6. ----- traitement selon identification / prédiction
 605       d_para = {
 606          'INTERPOL' : ['LIN', 'LIN'],
 607          'NOM_PARA' : 'INST',
 608          'PROL_DROITE' : 'CONSTANT',
 609          'PROL_GAUCHE' : 'CONSTANT',
 610       }
 611
 612       # 3.6.1. --- identification
 613       if identification:
 614          KJ_CRIT = l_crit[iocc]
 615          # on verifie que KJ_CRIT soit compris dans l'intervalle [KMOY_min, KMOY_max]
 616          valkmoy = tabres.KMOY.values()
 617          if not (min(valkmoy) <= KJ_CRIT <= max(valkmoy)):
 618 #                               'constant utilisé).')
 619             UTMESS('A','RUPTURE0_1')
 620          if is_2D:
 621             # définition des fonctions pour faire les interpolations
 622             d_para.update({ 'NOM_RESU' : 'DELTALMAX' })
 623             # DeltaMax en fonction du temps
 624             fdL = t_fonction(tb_Gpmax.INST.values(), tb_Gpmax.DELTALMAX.values(), d_para)
 625             # Gpmax fonction du temps
 626             d_para.update({ 'NOM_RESU' : 'GPMAX' })
 627             fGp = t_fonction(tb_Gpmax.INST.values(), tb_Gpmax.GPMAX.values(), d_para)
 628
 629             d_para.update({ 'NOM_PARA' : 'KMOY',
 630                            'NOM_RESU' : 'INST', })
 631             valkmoy = tabres.KMOY.values()
 632             # temps en fonction de Kmoy
 633             finv = t_fonction(valkmoy, tabres.INST.values(), d_para)
 634
 635             # valeurs à mettre dans la table
 636             # temps correspondant à KJ_CRIT
 637             ti   = finv(KJ_CRIT)
 638             # GP correspondant au temps critique
 639             Gpi  = fGp(ti)
 640             d_ident = {
 641                'KJ_CRIT'   : KJ_CRIT,
 642                'INST'      : ti,
 643                'GPMAX'     : Gpi,
 644                'KGPMAX'    : fKj(Gpi, **dict_constantes),
 645                'DELTALMAX' : fdL(ti),
 646             }
 647             lv_ident.append(d_ident)
 648          else:
 649
 650             l_i_noeuds_sommets = range(0,len(l_noeuds_fissure),pas)
 651             t_noeud_Kcrit = Table(para=tabres.para)
 652
 653             # On determine quel noeud sommet maximise KMOY au cours du temps:
 654
 655             row_KMOY_max = tabres.KMOY.MAXI()
 656             noeud_KMOY_max = row_KMOY_max.NOEUD.values()[0]
 657
 658             # avec le noeud ou KJ_CRIT est atteint, on regarde GP a gauche et a droite.
 659             # le GP le plus grand correspond au GPmax
 660             # on garde la tranche ou GP est le plus grand
 661
 662             id_noeud_KMOY_max = list(l_noeuds_fissure).index(noeud_KMOY_max)
 663             if id_noeud_KMOY_max==0:
 664                # "Gpmax sur la 1ere tranche"
 665                nume_tranche_Gpmax = 1
 666             elif id_noeud_KMOY_max==(len(l_noeuds_fissure)-1):
 667                # "Gpmax sur la derniere tranche"
 668                nume_tranche_Gpmax = nb_tranches
 669             else:
 670                # "Gpmax sur une tranche intermediaire"
 671                Gpi_tot = Table(para=tb_Gpmax.para)
 672                Gpi_gauche = tb_Gpmax.NUME_TRANCHE==(id_noeud_KMOY_max/pas)
 673                Gpi_tot.append(Gpi_gauche.rows[0])
 674                Gpi_droite = tb_Gpmax.NUME_TRANCHE==(id_noeud_KMOY_max/pas+1)
 675                Gpi_tot.append(Gpi_droite.rows[0])
 676                Gpi_tab = Gpi_tot.GPMAX.MAXI()
 677                nume_tranche_Gpmax = Gpi_tab.NUME_TRANCHE.values()[0]
 678
 679             # tb_Gpmax_TrancheCrit est une table de la meme nature que la table 2D tb_Gpmax
 680             # i.e. valeurs sur la seule tranche qui nous interesse (ou on sait
 681             # que KJ_CRIT sera atteint)
 682
 683             tb_Gpmax_TrancheCrit = tb_Gpmax.NUME_TRANCHE==nume_tranche_Gpmax
 684
 685             # avec le noeud ou KJ_CRIT est atteint, on determine le temps
 686             # critique par interpolation
 687             tabres_NoeudCrit = tabres.NOEUD==noeud_KMOY_max
 688
 689             # la suite est idem 2D, seuls les noms des tables changent
 690
 691             # définition des fonctions pour faire les interpolations
 692             d_para.update({ 'NOM_RESU' : 'DELTALMAX' })
 693             # DeltaMax en fonction du temps
 694             fdL = t_fonction(tb_Gpmax_TrancheCrit.INST.values(),
 695                               tb_Gpmax_TrancheCrit.DELTALMAX.values(),
 696                               d_para)
 697
 698             # Gpmax fonction du temps
 699             d_para.update({ 'NOM_RESU' : 'GPMAX' })
 700             fGp = t_fonction(tb_Gpmax_TrancheCrit.INST.values(),
 701                               tb_Gpmax_TrancheCrit.GPMAX.values(),
 702                               d_para)
 703
 704
 705             d_para.update({ 'NOM_PARA' : 'KMOY',
 706                            'NOM_RESU' : 'INST', })
 707             valkmoy = tabres_NoeudCrit.KMOY.values()
 708             # temps en fonction de Kmoy
 709             finv = t_fonction(valkmoy, tabres_NoeudCrit.INST.values(), d_para)
 710
 711             # valeurs à mettre dans la table
 712             # temps correspondant a KJ_CRIT
 713             ti   = finv(KJ_CRIT)
 714             # GP correspondant au temps critique
 715             Gpi  = fGp(ti)
 716             # par rapport a 2D, on ajoute 'NUME_TRANCHE'
 717             d_ident = {
 718                'KJ_CRIT'      : KJ_CRIT,
 719                'INST'         : ti,
 720                'NUME_TRANCHE' : int(nume_tranche_Gpmax),
 721                'GPMAX'        : Gpi,
 722                'KGPMAX'       : fKj(Gpi, **dict_constantes),
 723                'DELTALMAX'    : fdL(ti),
 724             }
 725             lv_ident.append(d_ident)
 726
 727       # 3.6.2. --- prédiction
 728       else:
 729          pass
 730
 731    # 4. ----- Construction de la table résultat si demandée
 732    # 4.1. --- identification
 733    if identification:
 734       if is_2D:
 735          para_tab_ident=('KJ_CRIT', 'INST', 'GPMAX', 'KGPMAX', 'DELTALMAX')
 736       else:
 737          para_tab_ident=('KJ_CRIT', 'INST', 'NUME_TRANCHE', 'GPMAX', 'KGPMAX', 'DELTALMAX')
 738       tab_ident = Table(rows=lv_ident,
 739                         para=para_tab_ident,
 740                         typ = ('R')*len(para_tab_ident),
 741                         titr='Identification aux valeurs de tenacités critiques')
 742       dprod_result = tab_ident.dict_CREA_TABLE()
 743       if info >= 2:
 744          print tab_ident
 745
 746    # 4.2. --- prédiction
 747    else:
 748       # définition de la fonction GPcrit = f(TEMP)
 749       d_para.update({ 'NOM_PARA' : 'TEMP',
 750                       'NOM_RESU' : 'GP_CRIT', })
 751       fGpc = t_fonction(mccalc['TEMP'], mccalc['GP_CRIT'], d_para)
 752
 753       # en 3D, GPMAX et DELTALMAX ne sont pas dans tab_glob: les recuperer de tb_Gpmax
 754
 755       if is_2D:
 756          tab_pred = tabl_glob['NUME_ORDRE', 'INST', 'TEMP', 'DELTALMAX', 'GPMAX']
 757       else:
 758          tab_pred = tb_Gpmax['NUME_ORDRE', 'INST', 'NUME_TRANCHE', 'DELTALMAX', 'GPMAX']
 759          # on recupere TEMP de tabl_glob
 760          tab_temp_tranche=Table(para=['NUME_ORDRE', 'NUME_TRANCHE', 'TEMP'])
 761          # on fait une moyenne de la temperature sur les noeuds d'une tranche
 762          for ordre in l_numord:
 763             tabl_glob_ORDRE_i = tabl_glob.NUME_ORDRE==ordre
 764             temp_ORDRE_i = tabl_glob_ORDRE_i.TEMP.values()
 765             for i_tranche in range(nb_tranches):
 766                l_temp_noeuds_tranche_i = temp_ORDRE_i[i_tranche : i_tranche+pas+1]
 767                temp_tranche_i = moyenne(*l_temp_noeuds_tranche_i)
 768                d = {'NUME_ORDRE': ordre, 'NUME_TRANCHE': i_tranche+1, 'TEMP': temp_tranche_i}
 769                tab_temp_tranche.append(d)
 770          tab_pred = merge(tab_pred, tab_temp_tranche, ['NUME_ORDRE', 'NUME_TRANCHE'])
 771
 772       tab_pred.fromfunction('GP_CRIT', fGpc, 'TEMP')
 773       tab_pred.fromfunction('PREDICTION', crit, ('GP_CRIT', 'GPMAX'))
 774       tab_pred.titr = 'Comparaison Gpmax à Gpcrit(T)'
 775       dprod_result = tab_pred.dict_CREA_TABLE()
 776
 777    # 9. ----- création de la table_sdaster résultat
 778    dprod = tabl_glob.dict_CREA_TABLE()
 779    result = CREA_TABLE(**dprod)
 780    tabresult = CREA_TABLE(**dprod_result)
 781
 782
 783
 784 # -----------------------------------------------------------------------------
 785 def CallRCVALE(TEMP, para, MATER):
 786    """Fonction appelant RCVALE et retourne la valeur d'un paramètre.
 787    """
 788    valres, flag_ok = MATER.RCVALE('ELAS', 'TEMP', TEMP, para)
 789    assert list(flag_ok).count('OK') != 1, \
 790          'Erreur lors de la récupération des valeurs du matériau.'
 791    return valres
 792
 793 # -----------------------------------------------------------------------------
 794 def fKj(G, YOUNG, NU):
 795    """Calcul de Kj à partir de G (formule d'Irwin)
 796    """
 797    Kj=(abs(G * YOUNG / (1.0 - NU**2)))**0.5
 798    return Kj
 799
 800 # -----------------------------------------------------------------------------
 801 def fDL(ICOP, pascop):
 802    """DeltaL = numéro copeau * pas d'entaille
 803    """
 804    return ICOP * pascop
 805
 806 # -----------------------------------------------------------------------------
 807 def fGp_Etot(TOTALE, DELTAL, R, syme=False):
 808    """Gp(Etotale, K), deltal pris dans le context global.
 809       ICOP      : numéro du copeau,
 810       DELTAL    : liste des epaisseurs des copeaux
 811       R         : rayon en axisymetrique,
 812                   longueur de l'élément 1D situé sur le front d'entaille si modèle 3D.
 813       syme      : True s'il y a symétrie.
 814    """
 815    import types
 816    fact_syme = 1.
 817    if syme:
 818       fact_syme = 2.
 819    Gp_Etot = fact_syme * TOTALE / (DELTAL * R )
 820    return Gp_Etot
 821
 822 # -----------------------------------------------------------------------------
 823 def MaxRelatif(table, nom_para):
 824    """Extrait le dernier maxi du champ `nom_para` de la table.
 825    """
 826    l_val = getattr(table, nom_para).values()
 827    l_val.reverse()
 828    Vlast = l_val[0]
 829    for val in l_val:
 830       if val < Vlast:
 831          break
 832       Vlast = val
 833    return getattr(table, nom_para) == Vlast
 834
 835 # -----------------------------------------------------------------------------
 836 def crit(GP_CRIT, GPMAX):
 837    """Retourne 1 quand GP_CRIT > GPMAX
 838    """
 839    if GPMAX > GP_CRIT:
 840       return 1
 841    else:
 842       return 0
 843
 844 # -----------------------------------------------------------------------------
 845 def moyenne(*args):
 846    """Fonction moyenne
 847    """
 848    return sum(args)/len(args)
 849
 850 def moyenne_positive(*args):
 851    """Fonction moyenne
 852    """
 853    return sum([abs(a) for a in args])/len(args)
 854
 855 def mysum(*args):
 856    """Fonction sum.
 857       La fonction sum ne peut pas etre appelee sur une liste de parametre
 858       d'une table via fromfunction
 859    """
 860    return sum(args)
 861
 862 def sum_and_check(*args):
 863    """Fonction sum.
 864       Verifie si la somme est positive.
 865       Si la somme est negative, on la met egale a zero
 866    """
 867    somme = sum(args)
 868    if somme<0:
 869       somme=0
 870    return somme
 871
 872 # On recupere des infos sur le fond de fissure
 873 def getFondFissInfo(fondfiss):
 874    # >FONFISS .FOND      .NOEU        <
 875    # >FONFISS .FOND      .TYPE        <
 876    import aster
 877    l_noeuds_fissure = aster.getvectjev(fondfiss.nom.ljust(8)+'.FOND      .NOEU        ')
 878    type_mailles = aster.getvectjev(fondfiss.nom.ljust(8)+'.FOND      .TYPE        ')
 879    if (type_mailles[0].strip() == 'SEG3' ):
 880       pas = 2
 881    else:
 882       pas = 1
 883    return l_noeuds_fissure, pas
 884
 885 ########################################################################
 886 # determination de la distance min entre 2 points consécutifs de la ligne de coupe
 887
 888 def largeur_tranche(nom_maillage, l_noms_noeuds_fissure, pas, i_tranche):
 889    # >MA      .COORDO    .VALE        <
 890    from math import sqrt
 891    import aster
 892
 893    # tuple des noms des noeuds du maillage
 894    t_noms_noeuds_maillage = aster.getvectjev(nom_maillage.ljust(8)+'.NOMNOE')
 895    # on convertit en liste pour utiliser la methode index
 896    # qui est plus optimal qu'une boucle sur les indices du tuple
 897    l_noms_noeuds_maillage = list(t_noms_noeuds_maillage)
 898
 899    l_numeros_noeuds_fissure = []
 900    for i in range(0,len(l_noms_noeuds_fissure),pas):
 901       nom = l_noms_noeuds_fissure[i]
 902       index = l_noms_noeuds_maillage.index(nom)
 903       l_numeros_noeuds_fissure.append(index)
 904
 905    coor1=aster.getvectjev(nom_maillage.ljust(8)+'.COORDO    .VALE        ',
 906                         l_numeros_noeuds_fissure[i_tranche]*3,3)
 907    coor2=aster.getvectjev(nom_maillage.ljust(8)+'.COORDO    .VALE        ',
 908                         l_numeros_noeuds_fissure[i_tranche+1]*3,3)
 909
 910    d=sqrt( (coor1[0]-coor2[0])**2+(coor1[1]-coor2[1])**2+(coor1[2]-coor2[2])**2)
 911    return d
 912
 913 def mergeLineInTable(multiTable, lineTable, nb_noeuds):
 914    # on ajoute a la table multiTable les colonnes de lineTable
 915    # pour chaque nume_ordre autant de fois qu'il y a de nb_noeuds
 916    from Utilitai.Table      import Table, merge
 917
 918    l_ordre = lineTable.NUME_ORDRE
 919    l_para = lineTable.copy().para
 920    l_para.remove('NUME_ORDRE')
 921    for i, ordre in enumerate(l_ordre):
 922       multiTable_i = multiTable.NUME_ORDRE==ordre
 923       row_i = lineTable.rows[i]
 924       for para in l_para:
 925          valeur_i = row_i[para]
 926          multiTable_i[para] = [valeur_i] * nb_noeuds
 927       if i==0:
 928          newTable=multiTable_i
 929       else:
 930          newTable = merge(newTable, multiTable_i)
 931
 932    return newTable