Aster/Cata/cataSTA8/Macro/post_dyna_alea_ops.py

   1 #@ MODIF post_dyna_alea_ops Macro  DATE 10/10/2006   AUTEUR MCOURTOI M.COURTOIS
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2006  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20
  21 def post_dyna_alea_ops(self,INTE_SPEC,NUME_VITE_FLUI,TOUT_ORDRE,NUME_ORDRE_I,
  22                        NOEUD_I,OPTION,MOMENT,TITRE,INFO,**args):
  23    import aster
  24    from types import ListType, TupleType
  25    EnumTypes = (ListType, TupleType)
  26    from Accas               import _F
  27    from Utilitai.Utmess     import UTMESS
  28    from Utilitai.t_fonction import t_fonction
  29    from Utilitai.Table      import Table
  30    import Numeric
  31    import math
  32    from math import pi,sqrt
  33
  34    commande='POST_DYNA_ALEA'
  35
  36    ier = 0
  37    # La macro compte pour 1 dans la numerotation des commandes
  38    self.set_icmd(1)
  39
  40    # Le concept sortant (de type table_sdaster ou dérivé) est tab
  41    self.DeclareOut('tabout', self.sd)
  42
  43    # On importe les definitions des commandes a utiliser dans la macro
  44    # Le nom de la variable doit etre obligatoirement le nom de la commande
  45    CREA_TABLE    = self.get_cmd('CREA_TABLE')
  46    CALC_TABLE    = self.get_cmd('CALC_TABLE')
  47    IMPR_TABLE    = self.get_cmd('IMPR_TABLE')
  48    RECU_FONCTION = self.get_cmd('RECU_FONCTION')
  49    IMPR_FONCTION = self.get_cmd('IMPR_FONCTION')
  50
  51    intespec=INTE_SPEC.EXTR_TABLE()
  52
  53 #  ------------------------------------------------------------------
  54 #  Liste des moments spectraux
  55 #  repérer le type de l'interspectre et son nom
  56 #                1- concept interspectre
  57 #                2- table de table d interspectre
  58
  59    if 'NUME_VITE_FLUI' in intespec.para :
  60       if TOUT_ORDRE!=None :
  61          jnuor=intespec['NUME_VITE_FLUI'].values()['NUME_VITE_FLUI']
  62          jvite=dict([(i,0) for i in jnuor]).keys()
  63       else :
  64         jvite=[NUME_VITE_FLUI,]
  65    else :
  66       jvite  =[None]
  67
  68 #  ------------------------------------------------------------------
  69 #  Repérer les couples d'indices selectionnés
  70 #  vérification de l'égalité du nombre d indices en i et j
  71
  72    if NUME_ORDRE_I!=None :
  73      l_ind_i=NUME_ORDRE_I
  74      l_ind_j=args['NUME_ORDRE_J']
  75      if type(l_ind_i) not in EnumTypes : l_ind_i=[l_ind_i]
  76      if type(l_ind_j) not in EnumTypes : l_ind_j=[l_ind_j]
  77      if len(l_ind_i)!=len(l_ind_j) :
  78         txt  = "il faut autant d indices en I et J"
  79         UTMESS('F',commande, txt)
  80      listpara=['NUME_ORDRE_I','NUME_ORDRE_J']
  81      listtype=['I','I']
  82      dicotabl={'NUME_ORDRE_I'  : l_ind_i  ,\
  83                'NUME_ORDRE_J'  : l_ind_j  , }
  84    elif NOEUD_I!=None :
  85      l_ind_i=NOEUD_I
  86      l_ind_j=args['NOEUD_J']
  87      l_cmp_i=args['NOM_CMP_I']
  88      l_cmp_j=args['NOM_CMP_J']
  89      if type(l_cmp_i) not in EnumTypes : l_cmp_i=[l_cmp_i]
  90      if type(l_cmp_j) not in EnumTypes : l_cmp_j=[l_cmp_j]
  91      if type(l_ind_i) not in EnumTypes : l_ind_i=[l_ind_i]
  92      if type(l_ind_j) not in EnumTypes : l_ind_j=[l_ind_j]
  93      if len(l_ind_i)!=len(l_ind_j) :
  94         txt  = "il faut autant d indices en I et J"
  95         UTMESS('F',commande, txt)
  96      if len(l_cmp_i)!=len(l_cmp_j) :
  97         txt  = "il faut autant de composantes en I et J"
  98         UTMESS('F',commande, txt)
  99      if len(l_ind_i)!=len(l_cmp_i) :
 100         txt  = "il faut autant de composantes que de noeuds"
 101         UTMESS('F',commande, txt)
 102      listpara=['NOEUD_I','NOEUD_J','NOM_CMP_I','NOM_CMP_J']
 103      listtype=['K8','K8','K8','K8',]
 104      dicotabl={'NOEUD_I'  : l_ind_i,\
 105                'NOEUD_J'  : l_ind_j,\
 106                'NOM_CMP_I': l_cmp_i,\
 107                'NOM_CMP_J': l_cmp_j }
 108 #  ------------------------------------------------------------------
 109 #  Cas de tous les indices centraux
 110
 111    elif OPTION!=None :
 112       if 'NUME_ORDRE_I' in intespec.para :
 113          inuor=intespec['NUME_ORDRE_I'].values()['NUME_ORDRE_I']
 114          imode=dict([(i,0) for i in inuor]).keys()
 115          l_ind_i=imode
 116          l_ind_j=imode
 117          listpara=['NUME_ORDRE_I','NUME_ORDRE_J']
 118          listtype=['I','I']
 119          dicotabl={'NUME_ORDRE_I'  : l_ind_i  ,\
 120                    'NUME_ORDRE_J'  : l_ind_j  , }
 121       else :
 122          if 'NUME_VITE_FLUI' in intespec.para :
 123             intespec=intespec.NUME_VITE_FLUI==jvite[0]
 124          l_ind_i=intespec['NOEUD_I'].values()['NOEUD_I']
 125          l_ind_j=intespec['NOEUD_J'].values()['NOEUD_J']
 126          if len(l_ind_i)!=len(l_ind_j) :
 127             txt  = "il faut autant d indices en I et J"
 128             UTMESS('F',commande, txt)
 129          l_cmp_i=intespec['NOM_CMP_I'].values()['NOM_CMP_I']
 130          l_cmp_j=intespec['NOM_CMP_J'].values()['NOM_CMP_J']
 131          if (len(l_ind_i)!=len(l_cmp_i) or len(l_ind_j)!=len(l_cmp_j)) :
 132             txt  = "il faut autant de composantes que de noeuds"
 133             UTMESS('F',commande, txt)
 134          l_l=zip(zip(l_ind_i,l_cmp_i),zip(l_ind_j,l_cmp_j))
 135          l_ind_i=[]
 136          l_ind_j=[]
 137          l_cmp_i=[]
 138          l_cmp_j=[]
 139          for ai,aj in l_l :
 140              if ai==aj :
 141                 l_ind_i.append(ai[0])
 142                 l_ind_j.append(aj[0])
 143                 l_cmp_i.append(ai[1])
 144                 l_cmp_j.append(aj[1])
 145          listpara=['NOEUD_I','NOEUD_J','NOM_CMP_I','NOM_CMP_J']
 146          listtype=['K8','K8','K8','K8',]
 147          dicotabl={'NOEUD_I'  : l_ind_i*len(jvite)  ,\
 148                    'NOEUD_J'  : l_ind_j*len(jvite)  ,\
 149                    'NOM_CMP_I': l_cmp_i*len(jvite)  ,\
 150                    'NOM_CMP_J': l_cmp_j*len(jvite) }
 151
 152    if jvite[0]!=None :
 153       listpara.append('NUME_VITE_FLUI')
 154       listtype.append('I')
 155       dicotabl['NUME_VITE_FLUI']=[]
 156 #  ------------------------------------------------------------------
 157 #  Liste des moments spectraux
 158
 159    l_moments=[0,1,2,3,4]
 160    if MOMENT!=None :
 161       l_moments=l_moments+list(MOMENT)
 162       l_moments=dict([(i,0) for i in l_moments]).keys()
 163
 164 #  ------------------------------------------------------------------
 165 #  Boucle sur les tables
 166
 167    l_ind=zip(l_ind_i,l_ind_j)
 168    for vite in jvite :
 169      if INFO==2 :
 170         texte='POUR LA MATRICE INTERSPECTRALE '+INTE_SPEC.nom+'\n'
 171         aster.affiche('MESSAGE',texte)
 172      for ind in l_ind :
 173         mcfact=[]
 174         if vite!=None :
 175           dicotabl['NUME_VITE_FLUI'].append(vite)
 176           mcfact.append(_F(NOM_PARA='NUME_VITE_FLUI',VALE_I=vite))
 177         if 'NOEUD_I' in listpara :
 178           mcfact.append(_F(NOM_PARA='NOEUD_I',VALE_K=ind[0]))
 179           mcfact.append(_F(NOM_PARA='NOEUD_I',VALE_K=ind[1]))
 180           if INFO==2 :
 181              aster.affiche('MESSAGE','INDICES :'+ind[0]+' - '+ind[1]+'\n')
 182         else :
 183           mcfact.append(_F(NOM_PARA='NUME_ORDRE_I',VALE_I=ind[0]))
 184           mcfact.append(_F(NOM_PARA='NUME_ORDRE_J',VALE_I=ind[1]))
 185         if INFO==2 :
 186              aster.affiche('MESSAGE','INDICES :'+str(ind[0])+' - '\
 187                                                 +str(ind[1])+'\n')
 188         __fon1=RECU_FONCTION(TABLE        = INTE_SPEC,
 189                              NOM_PARA_TABL= 'FONCTION_C',
 190                              FILTRE       = mcfact, )
 191         val  = __fon1.Valeurs()
 192         fvalx= Numeric.array(val[0])
 193         fvaly= Numeric.array(val[1])
 194         frez = fvalx[0]
 195
 196 #--- moments spectraux
 197
 198         val_mom={}
 199         for i_mom in l_moments :
 200             trapz     = Numeric.zeros(len(fvaly),Numeric.Float)
 201             trapz[0]  = 0.
 202             valy      = fvaly*(2*pi*fvalx)**i_mom
 203             trapz[1:] = (valy[1:]+valy[:-1])/2*(fvalx[1:]-fvalx[:-1])
 204             prim_y    = Numeric.cumsum(trapz)
 205             val_mom[i_mom] = prim_y[-1]
 206         for i_mom in l_moments :
 207           chmo='LAMBDA_'+str(i_mom).zfill(2)
 208           if dicotabl.has_key(chmo) : dicotabl[chmo].append(val_mom[i_mom])
 209           else :
 210                  dicotabl[chmo]=[val_mom[i_mom],]
 211                  listpara.append(chmo)
 212                  listtype.append('R')
 213
 214 #--- fonctions statistiques
 215
 216         pstat  = {'ECART'           :0.,\
 217                   'NB_PASS_ZERO_P_S':0.,\
 218                   'NB_EXTREMA_P_S'  :0.,\
 219                   'FACT_IRRE'       :0.,\
 220                   'FREQ_APPAR'      :0.,}
 221         if (NUME_VITE_FLUI or frez>=0.) :
 222 #--- cas NUME_VITE_FLUI, seule la partie positive du spectre est utilisée
 223 #--- Il faut donc doubler lambda  pour calculer le bon écart type
 224             pstat['ECART'] = sqrt(val_mom[0]*2.)
 225         else :
 226             pstat['ECART'] = sqrt(val_mom[0])
 227         if abs(val_mom[2])>=1e-20 :
 228               pstat['NB_EXTREMA_P_S'] = 1./pi*sqrt(val_mom[4]/val_mom[2])
 229         if abs(val_mom[0])>=1e-20 :
 230            pstat['NB_PASS_ZERO_P_S'] = 1./pi*sqrt(val_mom[2]/val_mom[0])
 231            pstat['FREQ_APPAR'] = 0.5*pstat['NB_PASS_ZERO_P_S']
 232            if abs(val_mom[4])>=1e-20 :
 233               pstat['FACT_IRRE'] = sqrt( val_mom[2]*val_mom[2]/val_mom[0]/val_mom[4])
 234
 235         for key in pstat.keys():
 236           if dicotabl.has_key(key) : dicotabl[key].append(pstat[key])
 237           else :
 238                  dicotabl[key]=[pstat[key],]
 239                  listpara.append(key)
 240                  listtype.append('R')
 241
 242 #--- construction de la table produite
 243
 244    mcfact=[]
 245    for i in range(len(listpara)) :
 246       if listtype[i]=='R':
 247          mcfact.append(_F(PARA=listpara[i] ,LISTE_R=dicotabl[listpara[i]] ))
 248       if listtype[i]=='K8':
 249          mcfact.append(_F(PARA=listpara[i] ,LISTE_K=dicotabl[listpara[i]] ))
 250       if listtype[i]=='I':
 251          mcfact.append(_F(PARA=listpara[i] ,LISTE_I=dicotabl[listpara[i]] ))
 252    tabout = CREA_TABLE(LISTE=mcfact,TITRE = 'POST_DYNA_ALEA concept : '+self.sd.nom)
 253
 254    return ier