Aster/Cata/cataSTA9/Macro/calc_spec_ops.py

   1 #@ MODIF calc_spec_ops Macro  DATE 21/10/2008   AUTEUR CORUS M.CORUS
   2
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2008  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 import copy
  22 import types
  23 from SD.sd_fonction import sd_fonction
  24
  25 # -----------------------------------------------------------------------------
  26 class FonctionError(Exception): pass
  27 class ParametreError(FonctionError):      pass  # probleme de NOM_PARA
  28 class InterpolationError(FonctionError):  pass
  29 class ProlongementError(FonctionError):   pass
  30
  31 # -----------------------------------------------------------------------------
  32
  33
  34 def calc_spec_ops(self,TAB_ECHANT,ECHANT,INTERSPE,TRANSFERT,TITRE,INFO,**args):
  35 #  ------------------------------------------------------------------
  36 #  Calcul d'une matrice interspectrale
  37 #  a partir de fonctions reelles
  38
  39    import aster
  40    from types import ListType, TupleType
  41    EnumTypes = (ListType, TupleType)
  42    from Accas               import _F
  43    from Utilitai.Utmess     import  UTMESS
  44    import Numeric
  45    import FFT
  46
  47    commande='CALC_SPEC'
  48
  49    ier = 0
  50    # La macro compte pour 1 dans la numerotation des commandes
  51    self.set_icmd(1)
  52
  53    # Le concept sortant (de type table_sdaster ou derive) est tab
  54    self.DeclareOut('tabout', self.sd)
  55
  56    # On importe les definitions des commandes a utiliser dans la macro
  57    # Le nom de la variable doit etre obligatoirement le nom de la commande
  58    CREA_TABLE    = self.get_cmd('CREA_TABLE')
  59    CALC_TABLE    = self.get_cmd('CALC_TABLE')
  60    DEFI_FONCTION = self.get_cmd('DEFI_FONCTION')
  61
  62 #--- Verifications sur les entrees --#
  63
  64    if (ECHANT==None and TAB_ECHANT==None) :
  65       raise FonctionError, 'Vous devez specifier des fonctions en entree'
  66
  67    if TAB_ECHANT==None : TAB_ECHANT=[]
  68    if ECHANT==None : ECHANT=[]
  69    if INTERSPE==None : INTERSPE=[]
  70    if TRANSFERT==None : TRANSFERT=[]
  71    if len(TAB_ECHANT)*len(ECHANT) !=0 :
  72       raise FonctionError, 'Vous pouvez specifier une table_fonction ou' + ' une liste de fonctions en entree, mais pas les deux'
  73    if len(TRANSFERT)*len(INTERSPE) !=0 :
  74       raise FonctionError, 'Vous ne pouvez specifier qu' +"'"+'un type de calcul par appel'
  75
  76
  77
  78 #-- Recuperation des entrees --#
  79
  80    l_f=[]
  81    l_t=[]
  82    l_G=[]
  83    l_H=[]
  84 #   for occ in TAB_ECHANT :
  85 #      l_t.append(('TAB_ECHANT',occ))
  86    if TAB_ECHANT:  #MC
  87       l_t = TAB_ECHANT.List_F()[0]
  88
  89    for occ in ECHANT :
  90       l_f.append(('ECHANT',occ))
  91    for occ in INTERSPE :
  92       l_G.append(('INTERSPE',occ))
  93    for occ in TRANSFERT :
  94       l_H.append(('TRANSFERT',occ))
  95
  96
  97 # Pour dimensionner les fenetres :
  98 # Cas ECHANT : on recupere simplement la premiere fonction
  99 # Cas_TAB_ECHANT : on recupere toutes les fonctions
 100
 101    if len(l_f) >0 :
 102       vale_sig=l_f[0][1]['FONCTION'].Valeurs();
 103       l_ech=len(vale_sig[0])
 104       dt=vale_sig[0][1]-vale_sig[0][0]
 105    else :
 106       #tab_ast=l_t[0][1]['NOM_TAB'];
 107       tab_ast=l_t['NOM_TAB']  #MC
 108       tab_py=tab_ast.EXTR_TABLE();
 109
 110       nom_fonc= tab_py['FONCTION'].values()['FONCTION']
 111       fonc_py = [sd_fonction(fonc) for fonc in nom_fonc]
 112       temp=fonc_py[0].VALE.get();
 113       dt=temp[1]-temp[0];
 114
 115       l_ech_temp=l_t['LONGUEUR_ECH'];
 116       recouvr_temp=l_t['RECOUVREMENT'];
 117       l_ech_t=[l_ech_temp[0]['DUREE'] , l_ech_temp[0]['POURCENT'],l_ech_temp[0]['NB_PTS'] ];
 118       recouvr_t=[recouvr_temp[0]['DUREE'] , recouvr_temp[0]['POURCENT'],recouvr_temp[0]['NB_PTS'] ];
 119       if l_ech_t.count(None)==3 : l_ech=len(temp)/2;
 120       if recouvr_t.count(None)==3 : recouvr=0;
 121       if l_ech_t.count(None)<2 :
 122          raise FonctionError, 'Vous ne pouvez utiliser qu'+"'"+ 'un mot clef pour definir la longueur des echantillons'
 123       if recouvr_t.count(None)<2 :
 124          raise FonctionError, 'Vous ne pouvez utiliser qu'+"'"+'un mot clef pour definir la longueur de recouvrement des echantillons'
 125       for i1 in range(3) :
 126           if l_ech_t[i1] !=None :
 127              if   i1 == 0 :
 128                 l_ech=int(Numeric.floor(l_ech_t[i1]/dt));
 129              elif i1 == 1 :
 130                 l_ech=int(Numeric.floor((len(temp)/2)*l_ech_t[i1]*0.01));
 131              elif i1 == 2 :
 132                 l_ech=int(Numeric.floor(l_ech_t[i1]))
 133       if l_ech > len(temp)/2 :
 134          raise FonctionError, 'Vous devez specifier une longueur d'+"'"+'echantillon inferieure a la longueur totale de l'+"'"+'acquisition'
 135       for i1 in range(3) :
 136           if recouvr_t[i1] !=None :
 137              if   i1 == 0 :
 138                 recouvr=int(Numeric.floor(recouvr_t[i1]/dt));
 139              elif i1 == 1 :
 140                 recouvr=int(Numeric.floor((l_ech)*recouvr_t[i1]*0.01));
 141              elif i1 == 2 :
 142                 recouvr=int(Numeric.floor(recouvr_t[i1]))
 143       if recouvr > l_ech :
 144          raise FonctionError, 'La longueur de recouvrement ne peut exceder la longueur '
 145
 146 #-- Recuperation des fenetres
 147
 148    for occ in l_G+l_H :
 149       if occ[1]['FENETRE'] == 'RECT' :
 150          fene=[1.]*l_ech
 151       elif occ[1]['FENETRE'] == 'HAMM' :
 152          fene=[0.54-0.46*Numeric.cos(2*Numeric.pi*i1/(l_ech-1)) for i1 in range(l_ech)]
 153       elif occ[1]['FENETRE'] == 'HANN' :
 154          fene=[0.5-0.5*Numeric.cos(2*Numeric.pi*i1/(l_ech-1)) for i1 in range(l_ech)]
 155       elif occ[1]['FENETRE'] == 'EXPO' :
 156          para=occ[1]['DEFI_FENE']
 157          if len(para) != 2 :
 158             raise FonctionError, 'Erreur de taille dans DEFI_FENE : ' + 'la fenetre exponentielle est definie par exactement deux valeurs'
 159          fene=[1.]*int(para[0]-1)+[Numeric.exp(para[1]*(i1-int(para[0]-1))*dt) for i1 in range(int(para[0]-1),l_ech)]
 160       elif occ[1]['FENETRE'] == 'PART' :
 161          fene=occ[1]['DEFI_FENE']
 162          if len(fene) != l_ech :
 163             raise FonctionError, 'Erreur de taille dans DEFI_FENE : ' + 'La fenetre doit etre definie avec le meme nombre de points que les echantillons'
 164
 165    if len(TRANSFERT)+len(INTERSPE) == 0 : #-- on ne rentre rien : interspectre par defaut - fenetre rectangulaire
 166       fene=[1.]*l_ech
 167       INTERSPE=1.;
 168
 169
 170 #--          Recuperation des signaux           --#
 171 #-- Verifications et transformations de Fourier --#
 172 #--         Entrees sous formes de table        --#
 173
 174    tmp=[];
 175    lt=[];
 176    frq=[];
 177    fft=[];
 178    df=[];
 179    num_ord=[];
 180    num_mes=[];
 181
 182
 183    if TAB_ECHANT : # Cas TAB_ECHANT
 184       num_mes_temp= tab_py['NUME_MES'].values()['NUME_MES']
 185       max_mes=Numeric.maximum.reduce(num_mes_temp);
 186       num_ord_temp= tab_py['NUME_ORDRE_I'].values()['NUME_ORDRE_I']
 187       long_fonc=[len(fonc_py[i1].VALE.get()) for i1 in range(len(fonc_py))]
 188
 189       N_fen=int(Numeric.floor((Numeric.minimum.reduce(long_fonc)/2-l_ech)/(l_ech-recouvr))+1)
 190
 191       sig=[];
 192       dt=[];
 193       for i1 in range(len(fonc_py)) :
 194          vale=fonc_py[i1].VALE.get();
 195          temp=(list(vale[0:int(len(vale)/2)]));
 196          sig.append(list(vale[int(len(vale)/2):]));
 197          test_pas=Numeric.subtract(temp[1:],temp[0:-1])
 198          crit=test_pas.tolist();
 199          crit.sort();
 200          dt.append(crit[-1]);
 201          if abs((crit[-1]-crit[0])/crit[-1]) > 1.e-5 :
 202             raise FonctionError, 'L'+"'"+'echantillonage doit etre fait a pas constant'
 203
 204       for j1 in range(N_fen) :
 205          for i1 in range(len(fonc_py)) :
 206             fft.append(FFT.fft(Numeric.multiply(sig[i1][j1*(l_ech-recouvr):(j1*(l_ech-recouvr)+l_ech)],fene)))
 207             if j1 == 0 : df.append(1./(dt[i1])/l_ech);
 208             num_mes.append(num_mes_temp[i1]+max_mes*j1);
 209             num_ord.append(num_ord_temp[i1]);
 210
 211       test_df=Numeric.subtract(df[1:],df[0:-1])
 212       test_df=test_df.tolist();
 213       test_df.sort();
 214       if abs(test_df[-1]) > 1.e-5 :
 215          raise FonctionError, 'Toutes les fonctions doivent etre definies ' + 'avec la meme frequence d'+"'"+'echantillonage'
 216
 217       frq = [df[-1]*i1 for i1 in range(l_ech)]
 218
 219
 220 #--          Recuperation des signaux           --#
 221 #-- Verifications et transformations de Fourier --#
 222 #--         Entrees sous formes de fonction     --#
 223
 224    if ECHANT:
 225       for occ in l_f :
 226          vale_sig=occ[1]['FONCTION'].Valeurs();
 227          #-- pour les tests ulterieurs --#
 228          lt.append(len(vale_sig[0]))
 229          if len(vale_sig[0]) != len(vale_sig[1]) :
 230             raise FonctionError, 'Les vecteurs associes au temps '+'et aux echantillons doivent etre de meme longueur'
 231          num_mes.append(occ[1]['NUME_MES'])
 232          num_ord.append(occ[1]['NUME_ORDRE_I'])
 233
 234          tmp.append(vale_sig[0])
 235          test_pas=Numeric.subtract(vale_sig[0][1:],vale_sig[0][0:-1])
 236          crit=test_pas.tolist();
 237          crit.sort();
 238          if abs((crit[-1]-crit[0])/crit[-1]) > 1.e-5 :
 239             raise FonctionError, 'L'+"'"+'echantillonage doit etre fait a pas constant'
 240          fft.append(FFT.fft(Numeric.multiply(vale_sig[1],fene)))
 241          df.append(1./(crit[-1])/len(vale_sig[0]));
 242
 243
 244       #-- Verification des longueurs --#
 245
 246       test_long=Numeric.subtract(lt[1:],lt[0:-1])
 247       test_long=test_long.tolist();
 248       test_long.sort();
 249       if (test_long[-1]-test_long[0]) != 0 :
 250          raise FonctionError, 'Toutes les fonctions doivent etre definies avec le meme nombre de points'
 251
 252
 253       test_df=Numeric.subtract(df[1:],df[0:-1])
 254       test_df=test_df.tolist();
 255       test_df.sort();
 256       if abs(test_df[-1]) > 1.e-5 :
 257          raise FonctionError, 'Toutes les fonctions doivent etre definies '+'avec la meme frequence d'+"'"+'echantillonage'
 258
 259       frq = [df[-1]*i1 for i1 in range(lt[-1])]
 260
 261
 262 #-- index des numeros d'ordre pour le moyennage
 263
 264    uu=[];
 265    vv=[];
 266    uu=uu+num_ord;
 267    vv=vv+num_ord;
 268    uu.sort();
 269    ind_ord=[];
 270    list_ord=[];
 271    while  len(uu) > 0 :
 272       list_ord.append(uu[0])
 273       tt=[];
 274       for i1 in range(uu.count(uu[0])) :
 275          tt.append(vv.index(uu[0]))
 276          vv[tt[-1]]=0
 277       ind_ord.append(tt)
 278       uu=uu[int(uu.count(uu[0])):]
 279
 280 #-- Calcul de la matrice inter spectrale
 281
 282    if len(INTERSPE) != 0 :
 283       dimh   = (len(list_ord)*(len(list_ord)+1))/2
 284       l_fc=[];
 285       nume_i1=[]
 286       nume_j1=[]
 287
 288       for i1 in range(len(list_ord)) :
 289          for j1 in range(i1,len(list_ord)) :
 290             #-- on ne calcule les spectres que pour des numeros de mesures correspondants
 291             #-- Ca n'a a priori pas de sens de calculer l'interspectre entre deux signaux acquis a des instants differents
 292             #-- Par contre, on peut moyenner deux interspectres obtenus a des instants differents, sous reserve
 293             #-- de stationnarite et d'ergodicite du signal
 294             mes_i1=[num_mes[k1] for k1 in ind_ord[i1]]
 295             mes_j1=[num_mes[k1] for k1 in ind_ord[j1]]
 296             ind_mes=[];
 297             #-- recuperation des indices des fft a prendre en compte pour l'interspectre
 298             for k1 in range(len(mes_i1)) :
 299                if mes_i1[k1] in mes_j1 :
 300                   ind_mes.append([ind_ord[i1][k1],ind_ord[j1][mes_j1.index(mes_i1[k1])]])
 301
 302             #-- Calcul des interspectres
 303             dsp=[0.j]*l_ech;
 304             if len(ind_mes) > 0 :
 305                for l1 in range(len(ind_mes)) :
 306                   dsp_t=Numeric.multiply(Numeric.conjugate(fft[ind_mes[l1][0]]),fft[ind_mes[l1][1]])
 307                   dsp_t=Numeric.divide(dsp_t,l_ech*len(ind_mes))
 308                   dsp=Numeric.add(dsp,dsp_t)
 309                dsp=dsp.tolist();
 310                dsp_r=[];
 311
 312                for k1 in range(int(Numeric.floor(l_ech/2))) :
 313                   dsp_r=dsp_r+[frq[k1],dsp[k1].real,dsp[k1].imag]
 314
 315                _fonc = DEFI_FONCTION(NOM_PARA='FREQ',VALE_C=dsp_r,);
 316                l_fc.append(_fonc.nom)
 317                nume_i1.append(list_ord[i1])
 318                nume_j1.append(list_ord[j1])
 319
 320       mcfact=[]
 321       mcfact.append(_F(PARA='NOM_CHAM'    ,LISTE_K='DSP'  ))
 322       mcfact.append(_F(PARA='OPTION'      ,LISTE_K='TOUT' ))
 323       mcfact.append(_F(PARA='DIMENSION'   ,LISTE_I=(dimh,) ))
 324       mcfact.append(_F(PARA='NUME_ORDRE_I',LISTE_I=nume_i1 ))
 325       mcfact.append(_F(PARA='NUME_ORDRE_J',LISTE_I=nume_j1 ))
 326       mcfact.append(_F(PARA='FONCTION_C'  ,LISTE_K=l_fc  ))
 327       self.DeclareOut('tab_inte',self.sd)
 328       tab_inte=CREA_TABLE(LISTE=mcfact,
 329                           TITRE='',
 330                           TYPE_TABLE='TABLE_FONCTION')
 331
 332
 333 #-- Calcul des transferts
 334
 335    if len(TRANSFERT) != 0 :
 336
 337       l_fc=[];
 338       nume_i1=[]
 339       nume_j1=[]
 340
 341       #-- test sur les entrees pour les references --#
 342       if type(l_H[0][1]['REFER'])==int :
 343          refer=[];
 344          refer.append(l_H[0][1]['REFER'])
 345       elif type(l_H[0][1]['REFER'])==tuple :
 346          refer=list(l_H[0][1]['REFER'])
 347
 348       ind_refer=[];
 349       dimh   = len(refer)*(len(list_ord)-len(refer))
 350       for k1 in range(len(refer)) :
 351          for l1 in range(len(list_ord)) :
 352             if refer[k1] == list_ord[l1] : ind_refer.append(l1);
 353
 354       #-- H1 : interspectre / autospectre
 355       #-- H2 : autospectre / interspectre
 356       #-- CO : coherence entre H1 et H2.
 357
 358       if l_H[0][1]['ESTIM']!='HV' :
 359          for i1 in range(len(refer)) :
 360             for j1 in range(len(list_ord)) :
 361                if refer[i1] != list_ord[j1] :
 362                   mes_i1=[num_mes[k1] for k1 in ind_ord[ind_refer[i1]]]  #-- mesures des efforts
 363                   mes_j1=[num_mes[k1] for k1 in ind_ord[j1]]  #-- mesures des reponses
 364                   ind_mes=[];
 365                   #-- recuperation des indices des mesures a predre en compte pour les spectres
 366                   for k1 in range(len(mes_i1)) :
 367                      if mes_i1[k1] in mes_j1 :
 368                         ind_mes.append([ind_ord[i1][k1],ind_ord[j1][mes_j1.index(mes_i1[k1])]])
 369
 370                   #-- Calcul des FRF
 371                   Guu=[0.j]*l_ech;
 372                   Gyy=[0.j]*l_ech;
 373                   Gyu=[0.j]*l_ech;
 374                   if len(ind_mes) > 0 :
 375                      for l1 in range(len(ind_mes)) :
 376                         Guu_t=Numeric.multiply(Numeric.conjugate(fft[ind_mes[l1][0]]),fft[ind_mes[l1][0]])
 377                         Guu=Numeric.add(Guu,Guu_t)
 378                         Gyu_t=Numeric.multiply(Numeric.conjugate(fft[ind_mes[l1][1]]),fft[ind_mes[l1][0]])
 379                         Gyu=Numeric.add(Gyu,Gyu_t)
 380                         Gyy_t=Numeric.multiply(Numeric.conjugate(fft[ind_mes[l1][1]]),fft[ind_mes[l1][1]])
 381                         Gyy=Numeric.add(Gyy,Gyy_t)
 382
 383                      if l_H[0][1]['ESTIM']=='H1' :
 384                         frf=Numeric.divide(Numeric.conjugate(Gyu),Guu);
 385                         nom_frf='FRF-H1';
 386                      elif l_H[0][1]['ESTIM']=='H2' :
 387                         frf=Numeric.divide(Gyy,Gyu);
 388                         nom_frf='FRF-H2';
 389                      elif l_H[0][1]['ESTIM']=='CO' :
 390                         H1=Numeric.divide(Numeric.conjugate(Gyu),Guu);
 391                         H2=Numeric.divide(Gyy,Gyu);
 392                         frf=Numeric.divide(H1,H2);
 393                         nom_frf='FRF-COH';
 394
 395                      frf=frf.tolist();
 396                      frf_r=[];
 397
 398                      for k1 in range(int(Numeric.floor(l_ech/2))) :
 399                         frf_r=frf_r+[frq[k1],frf[k1].real,frf[k1].imag]
 400
 401                      _fonc = DEFI_FONCTION(NOM_PARA='FREQ',VALE_C=frf_r,);
 402                      l_fc.append(_fonc.nom)
 403                      nume_i1.append(refer[i1])
 404                      nume_j1.append(list_ord[j1])
 405
 406       #-- On remplit la table_fonction avec tout ce qui va bien
 407
 408       mcfact=[]
 409       mcfact.append(_F(PARA='NOM_CHAM'    ,LISTE_K=nom_frf ))
 410       mcfact.append(_F(PARA='OPTION'      ,LISTE_K='TOUT'  ))
 411       mcfact.append(_F(PARA='DIMENSION'   ,LISTE_I=(dimh,) ))
 412       mcfact.append(_F(PARA='NUME_ORDRE_I',LISTE_I=nume_i1 ))
 413       mcfact.append(_F(PARA='NUME_ORDRE_J',LISTE_I=nume_j1 ))
 414       mcfact.append(_F(PARA='FONCTION_C'  ,LISTE_K=l_fc  ))
 415       self.DeclareOut('tab_inte',self.sd)
 416       tab_inte=CREA_TABLE(LISTE=mcfact,
 417                           TITRE='',
 418                           TYPE_TABLE='TABLE_FONCTION')
 419