Aster/Cata/cataSTA9/Macro/post_gp_ops.py

   1 #@ MODIF post_gp_ops Macro  DATE 31/10/2006   AUTEUR REZETTE C.REZETTE
   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    macro = 'POST_GP'
  31    ier=0
  32    from Accas import _F
  33    from Utilitai.Utmess     import UTMESS
  34    from Utilitai.Table      import Table, merge
  35    from Utilitai.t_fonction import t_fonction
  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
  45    # ----- Comptage, commandes + déclaration concept sortant
  46    self.set_icmd(1)
  47    self.DeclareOut('result', self.sd)
  48    self.DeclareOut('tabresult', self['TABL_RESU'])
  49    info = self['INFO']
  50
  51    # 0. ----- Type de calcul
  52    identification = self['IDENTIFICATION'] != None
  53    if identification:
  54       # 0.1. --- identification : on boule sur les valeurs de TEMP.
  55       #          Pour chaque couple (T, Kjc(T)), on évalue les Ki, Kmoy et
  56       #          les valeurs de Gpmax, DeltaLmax, inst.max correspondantes.
  57       mccalc = self['IDENTIFICATION']
  58       l_crit = mccalc['KJ_CRIT']
  59       lv_ident = []
  60       l_temp = mccalc['TEMP']
  61    else:
  62       # 0.2. --- prédiction : on ne fait qu'une itération.
  63       #          Il faut un RESU_THER (sinon on utilise la température du
  64       #          premier Gpcrit et cà n'a pas trop d'intéret).
  65       #          A chaque instant, on regarde à quelle température est le
  66       #          fond d'entaille et on compare Gpmax à cet instant au Gpcrit.
  67       mccalc = self['PREDICTION']
  68       l_crit = mccalc['GP_CRIT']
  69       lv_pred = []
  70       l_temp = mccalc['TEMP'][0]
  71
  72    if not type(l_temp) in EnumTypes:
  73       l_temp = [l_temp,]
  74    if not type(l_crit) in EnumTypes:
  75       l_crit = [l_crit,]
  76
  77    # 1. ----- calcul de G-theta
  78    nbcour = len(self['THETA_2D'])
  79    l_tab = []
  80    for occ in self['THETA_2D']:
  81       dMC = occ.cree_dict_valeurs(occ.mc_liste)
  82
  83       __theta = CALC_THETA(MODELE=self['MODELE'],
  84                            DIRECTION=self['DIRECTION'],
  85                            THETA_2D=_F(GROUP_NO=dMC['GROUP_NO'],
  86                                        MODULE=1.0,
  87                                        R_INF=dMC['R_INF'],
  88                                        R_SUP=dMC['R_SUP']),)
  89
  90       __gtheta = CALC_G(THETA=_F(THETA=__theta),
  91                         EXCIT=self['EXCIT'].List_F(),
  92                         RESULTAT=self['RESULTAT'],
  93                         TOUT_ORDRE='OUI',
  94                         SYME_CHAR=self['SYME_CHAR'],
  95                         COMP_ELAS=self['COMP_ELAS'].List_F(),)
  96
  97       tab = __gtheta.EXTR_TABLE()
  98
  99       # une Table par couronne
 100       l_tab.append(tab)
 101
 102    # 2. ----- Calcul de l'energie élastique en exploitant les groupes de
 103    #          mailles fournis par la procedure de maillage
 104    l_copo = [grma.strip() for grma in self['GROUP_MA']]
 105    nbcop = len(l_copo)
 106    l_charg = [charg['CHARGE'] for charg in self['EXCIT']]
 107
 108    __ener = POST_ELEM(MODELE=self['MODELE'],
 109                         RESULTAT=self['RESULTAT'],
 110                         CHARGE=l_charg,
 111                         TOUT_ORDRE='OUI',
 112                         ENER_ELAS=_F(GROUP_MA=l_copo),
 113                         TITRE='Energie élastique',)
 114
 115    t_enel = __ener.EXTR_TABLE()
 116
 117    # 2.1. ----- Indice de chaque copeau et deltaL
 118    d_icop = dict(zip(l_copo, range(1, nbcop + 1)))
 119
 120    l_lieu = [grma.strip() for grma in t_enel.LIEU.values()]
 121    l_icop = [d_icop[grma] for grma in l_lieu]
 122    t_enel['ICOP'] = l_icop
 123    t_enel.fromfunction('DELTAL', fDL, 'ICOP', { 'pascop' : self['PAS_ENTAILLE'] })
 124
 125    # 2.2. ----- Calcul de Gp fonction de Ener.Totale et de deltaL
 126    t_enel.fromfunction('GP', fGp_Etot, ('TOTALE', 'ICOP'),
 127          { 'pascop' : self['PAS_ENTAILLE'],
 128            'syme'   : self['SYME_CHAR'] != 'SANS',
 129            'R'      : self['RAYON_AXIS'] })
 130
 131    # 2.3. ----- Tableau de Gp = f(icop) pour chaque instant
 132    if info >= 2:
 133       tGp_t_icop = t_enel['INST', 'DELTAL', 'GP']
 134       tGp_t_icop.titr = "Gp à chaque instant en fonction de la distance au " \
 135                         "fond d'entaille"
 136       tGp_t_icop.ImprTabCroise()
 137
 138    # 2.4. ----- Table Gpmax
 139    ttmp = t_enel['NUME_ORDRE', 'INST', 'ICOP', 'DELTAL', 'GP']
 140    l_numord = list(Set(ttmp.NUME_ORDRE.values()))
 141    l_numord.sort()
 142    for j in l_numord:
 143       tj = ttmp.NUME_ORDRE == j
 144       if self['CRIT_MAXI_GP'] == 'ABSOLU':
 145          t = tj.GP.MAXI()
 146       else:
 147          t = MaxRelatif(tj, 'GP')
 148       if j == 1:
 149          tb_Gpmax = t
 150       else:
 151          tb_Gpmax = tb_Gpmax | t
 152    tb_Gpmax.Renomme('GP', 'GPMAX')
 153    tb_Gpmax.Renomme('ICOP', 'ICOPMAX')
 154    tb_Gpmax.Renomme('DELTAL', 'DELTALMAX')
 155    tb_Gpmax.titr = 'Gpmax à chaque instant'
 156    if info >= 2:
 157       print tb_Gpmax
 158
 159    # 2.5. ----- extraction de la température en fond d'entaille
 160    if self['RESU_THER']:
 161       grno_fond = self['THETA_2D'][0]['GROUP_NO']
 162       __relev = POST_RELEVE_T(ACTION=_F(RESULTAT=self['RESU_THER'],
 163                                         OPERATION='EXTRACTION',
 164                                         INTITULE='Temperature',
 165                                         NOM_CHAM='TEMP',
 166                                         TOUT_ORDRE='OUI',
 167                                         NOM_CMP='TEMP',
 168                                         GROUP_NO=grno_fond,),)
 169       t_relev = __relev.EXTR_TABLE()['NUME_ORDRE', 'TEMP']
 170
 171
 172    # 3. ----- boucle sur les mots-clés facteurs
 173    #          opérations dépendant de la température
 174    MATER = self['MATER']
 175    flag_mat = True
 176
 177    for iocc, TEMP in enumerate(l_temp):
 178       # 3.0. ----- Temperature fonction du temps : si on n'a pas de RESU_THER,
 179       #            on prend la température d'identification.
 180       if not self['RESU_THER']:
 181          l_rows = [{'NUME_ORDRE' : i, 'TEMP' : TEMP} for i in l_numord]
 182          t_relev = Table(rows=l_rows, para=('NUME_ORDRE', 'TEMP'), typ=('R', 'R'))
 183          flag_mat = True
 184
 185       # 3.1. ----- extrait du matériau E(TEMP) et NU(TEMP) (si nécessaire)
 186       if flag_mat:
 187          t_relev.fromfunction('YOUNG', CallRCVALE, 'TEMP',
 188                { 'para' : 'E', 'MATER' : MATER })
 189          t_relev.fromfunction('NU', CallRCVALE, 'TEMP',
 190                { 'para' : 'NU', 'MATER' : MATER })
 191          #tb_Gpmax = merge(tb_Gpmax, t_relev, 'NUME_ORDRE')
 192          flag_mat = False
 193
 194       # 3.2. ----- paramètres
 195       dict_constantes = {
 196          'YOUNG' : CallRCVALE(TEMP, 'E', MATER),
 197          'NU'    : CallRCVALE(TEMP, 'NU', MATER),
 198          'R'     : self['RAYON_AXIS'],
 199       }
 200
 201       # 3.3. ----- calcul de Kj(G)
 202       l_tabi = []
 203       for k, tab in enumerate(l_tab):
 204          # fusion avec TEMP, E et nu.
 205          tab = merge(tab, t_relev, 'NUME_ORDRE')
 206
 207          # calcul de Kj(G) = K_i
 208          new_para = 'K_%d' % (k + 1)
 209          tab.fromfunction(new_para, fKj, ('G', 'YOUNG', 'NU'),
 210                           { 'R' : self['RAYON_AXIS'] })
 211
 212          # renomme G en G_i
 213          tab.Renomme('G', 'G_%d' % (k + 1))
 214          l_tabi.append(tab)
 215
 216       # 3.4 ----- Table des Gi, Ki sur les differentes couronnes + Kmoyen
 217       tabK_G = l_tabi[0]['NUME_ORDRE']
 218       for tab in l_tabi:
 219          tabK_G = merge(tabK_G, tab, 'NUME_ORDRE')
 220       tabK_G.titr = 'G et K sur les differentes couronnes + moyennes'
 221       tabK_G.fromfunction('GMOY', moyenne, ['G_%d' % (k + 1) for k in range(nbcour)])
 222       tabK_G.fromfunction('KMOY', moyenne, ['K_%d' % (k + 1) for k in range(nbcour)])
 223
 224       # 3.5. ----- Contribution à la table globale
 225       tabres = merge(tabK_G, tb_Gpmax, 'NUME_ORDRE')
 226       tabres['OCCURRENCE'] = [iocc + 1] * len(l_numord)
 227       if info >= 2:
 228          print tabres
 229
 230       # 3.5.1. --- Table globale
 231       if iocc == 0:
 232          tabl_glob = tabres
 233       else:
 234          tabl_glob = merge(tabl_glob, tabres)
 235       tabl_glob.titr = 'G, K sur les differentes couronnes, Gmoy, Kmoy et ' \
 236                        'Gpmax fonctions du temps'
 237
 238       # 3.6. ----- traitement selon identification / prédiction
 239       d_para = {
 240          'INTERPOL' : ['LIN', 'LIN'],
 241          'NOM_PARA' : 'INST',
 242          'PROL_DROITE' : 'CONSTANT',
 243          'PROL_GAUCHE' : 'CONSTANT',
 244       }
 245       # Gpmax fonction du temps
 246       d_para.update({ 'NOM_RESU' : 'GPMAX' })
 247       fGp = t_fonction(tabres.INST.values(), tabres.GPMAX.values(), d_para)
 248
 249       # 3.6.1. --- identification
 250       if identification:
 251          KJ_CRIT = l_crit[iocc]
 252          # définition des fonctions pour faire les interpolations
 253          d_para.update({ 'NOM_RESU' : 'DELTALMAX' })
 254          fdL = t_fonction(tabres.INST.values(), tabres.DELTALMAX.values(), d_para)
 255
 256          d_para.update({ 'NOM_PARA' : 'KMOY',
 257                          'NOM_RESU' : 'INST', })
 258          valkmoy = tabres.KMOY.values()
 259          finv = t_fonction(valkmoy, tabres.INST.values(), d_para)
 260
 261          if not (min(valkmoy) <= KJ_CRIT <= max(valkmoy)):
 262             UTMESS('A', macro, 'Interpolation hors du domaine (prolongement ' \
 263                                'constant utilisé).')
 264          # valeurs à mettre dans la table
 265          ti   = finv(KJ_CRIT)
 266          Gpi  = fGp(ti)
 267          d_ident = {
 268             'KJ_CRIT'   : KJ_CRIT,
 269             'INST'      : ti,
 270             'GPMAX'     : Gpi,
 271             'KGPMAX'    : fKj(Gpi, **dict_constantes),
 272             'DELTALMAX' : fdL(ti),
 273          }
 274          lv_ident.append(d_ident)
 275
 276       # 3.6.2. --- prédiction
 277       else:
 278          pass
 279
 280    # 4. ----- Construction de la table résultat si demandée
 281    # 4.1. --- identification
 282    if identification:
 283       tab_ident = Table(rows=lv_ident,
 284                         para=('KJ_CRIT', 'INST', 'GPMAX', 'KGPMAX', 'DELTALMAX'),
 285                         typ= ('R',       'R',    'R',     'R',      'R'),
 286                         titr='Identification aux valeurs de tenacités critiques')
 287       dprod_result = tab_ident.dict_CREA_TABLE()
 288       if info >= 2:
 289          print tab_ident
 290
 291    # 4.2. --- prédiction
 292    else:
 293       # définition de la fonction GPcrit = f(TEMP)
 294       d_para.update({ 'NOM_PARA' : 'TEMP',
 295                       'NOM_RESU' : 'GP_CRIT', })
 296       fGpc = t_fonction(mccalc['TEMP'], mccalc['GP_CRIT'], d_para)
 297
 298       tab_pred = tabl_glob['NUME_ORDRE', 'INST', 'TEMP', 'DELTALMAX', 'GPMAX']
 299       tab_pred.fromfunction('GP_CRIT', fGpc, 'TEMP')
 300       tab_pred.fromfunction('PREDICTION', crit, ('GP_CRIT', 'GPMAX'))
 301       tab_pred.titr = 'Comparaison Gpmax à Gpcrit(T)'
 302       dprod_result = tab_pred.dict_CREA_TABLE()
 303
 304    # 9. ----- création de la table_sdaster résultat
 305    dprod = tabl_glob.dict_CREA_TABLE()
 306    result = CREA_TABLE(**dprod)
 307    tabresult = CREA_TABLE(**dprod_result)
 308
 309
 310
 311 # -----------------------------------------------------------------------------
 312 def CallRCVALE(TEMP, para, MATER):
 313    """Fonction appelant RCVALE et retourne la valeur d'un paramètre.
 314    """
 315    valres, flag_ok = MATER.RCVALE('ELAS', 'TEMP', TEMP, para)
 316    assert list(flag_ok).count('OK') != 1, \
 317          'Erreur lors de la récupération des valeurs du matériau.'
 318    return valres
 319
 320 # -----------------------------------------------------------------------------
 321 def fKj(G, YOUNG, NU, R):
 322    """Calcul de Kj à partir de G (formule d'Irwin)
 323    """
 324    return (G / R * YOUNG / (1.0 - NU**2))**0.5
 325
 326 # -----------------------------------------------------------------------------
 327 def fDL(ICOP, pascop):
 328    """DeltaL = numéro copeau * pas d'entaille
 329    """
 330    return ICOP * pascop
 331
 332 # -----------------------------------------------------------------------------
 333 def fGp_Etot(TOTALE, ICOP, pascop, R, syme=False):
 334    """Gp(Etotale, K), deltal pris dans le context global.
 335       ICOP   : numéro du copeau,
 336       pascop : pas d'entaille.
 337       syme   : True s'il y a symétrie.
 338    """
 339    fact_axis = 1.
 340    if syme:
 341       fact_axis = 2.
 342    return fact_axis * TOTALE / (fDL(ICOP, pascop) * R)
 343
 344 # -----------------------------------------------------------------------------
 345 def MaxRelatif(table, nom_para):
 346    """Extrait le dernier maxi du champ `nom_para` de la table.
 347    """
 348    l_val = getattr(table, nom_para).values()
 349    l_val.reverse()
 350    Vlast = l_val[0]
 351    for val in l_val:
 352       if val < Vlast:
 353          break
 354       Vlast = val
 355    return getattr(table, nom_para) == Vlast
 356
 357 # -----------------------------------------------------------------------------
 358 def crit(GP_CRIT, GPMAX):
 359    """Retourne 1 quand GP_CRIT > GPMAX
 360    """
 361    if GPMAX > GP_CRIT:
 362       return 1
 363    else:
 364       return 0
 365
 366 # -----------------------------------------------------------------------------
 367 def moyenne(*args):
 368    """Fonction moyenne
 369    """
 370    return sum(args)/len(args)
 371