Aster/Cata/cataSTA10/Macro/calc_spec_ops.py

   1 #@ MODIF calc_spec_ops Macro  DATE 11/05/2010   AUTEUR COURTOIS M.COURTOIS
   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 numpy
  45    import numpy.fft as 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       print "test : l_ech = ", l_ech
 106    else :
 107       #tab_ast=l_t[0][1]['NOM_TAB'];
 108       tab_ast=l_t['NOM_TAB']  #MC
 109       tab_py=tab_ast.EXTR_TABLE();
 110
 111       nom_fonc= tab_py['FONCTION'].values()['FONCTION']
 112       fonc_py = [sd_fonction(fonc) for fonc in nom_fonc]
 113       temp=fonc_py[0].VALE.get();
 114       dt=temp[1]-temp[0];
 115
 116       l_ech_temp=l_t['LONGUEUR_ECH'];
 117       recouvr_temp=l_t['RECOUVREMENT'];
 118
 119       l_ech_t=[l_ech_temp[0]['DUREE'] , l_ech_temp[0]['POURCENT'],l_ech_temp[0]['NB_PTS'] ];
 120       recouvr_t=[recouvr_temp[0]['DUREE'] , recouvr_temp[0]['POURCENT'],recouvr_temp[0]['NB_PTS'] ];
 121
 122       if l_ech_t.count(None)==3 : l_ech=len(temp)/2;
 123       if recouvr_t.count(None)==3 : recouvr=0;
 124       if l_ech_t.count(None)<2 :
 125          raise FonctionError, 'Vous ne pouvez utiliser qu'+"'"+ 'un mot clef pour definir la longueur des echantillons'
 126       if recouvr_t.count(None)<2 :
 127          raise FonctionError, 'Vous ne pouvez utiliser qu'+"'"+'un mot clef pour definir la longueur de recouvrement des echantillons'
 128       for i1 in range(3) :
 129           if l_ech_t[i1] !=None :
 130              if   i1 == 0 :
 131                 l_ech=int(numpy.floor(l_ech_t[i1]/dt));
 132              elif i1 == 1 :
 133                 l_ech=int(numpy.floor((len(temp)/2)*l_ech_t[i1]*0.01));
 134              elif i1 == 2 :
 135                 l_ech=int(numpy.floor(l_ech_t[i1]))
 136       if l_ech > len(temp)/2 :
 137          raise FonctionError, 'Vous devez specifier une longueur d'+"'"+'echantillon inferieure a la longueur totale de l'+"'"+'acquisition'
 138       for i1 in range(3) :
 139           if recouvr_t[i1] !=None :
 140              if   i1 == 0 :
 141                 recouvr=int(numpy.floor(recouvr_t[i1]/dt));
 142              elif i1 == 1 :
 143                 recouvr=int(numpy.floor((l_ech)*recouvr_t[i1]*0.01));
 144              elif i1 == 2 :
 145                 recouvr=int(numpy.floor(recouvr_t[i1]))
 146       if recouvr > l_ech :
 147          raise FonctionError, 'La longueur de recouvrement ne peut exceder la longueur '
 148       print "test2 : l_ech = ", l_ech
 149
 150 #-- Recuperation des fenetres
 151
 152
 153    for occ in l_G+l_H :
 154       if occ[1]['FENETRE'] == 'RECT' :
 155          fene=[1.]*l_ech
 156       elif occ[1]['FENETRE'] == 'HAMM' :
 157          fene=[0.54-0.46*numpy.cos(2*numpy.pi*i1/(l_ech-1)) for i1 in range(l_ech)]
 158       elif occ[1]['FENETRE'] == 'HANN' :
 159          fene=[0.5-0.5*numpy.cos(2*numpy.pi*i1/(l_ech-1)) for i1 in range(l_ech)]
 160       elif occ[1]['FENETRE'] == 'EXPO' :
 161          para=occ[1]['DEFI_FENE']
 162          if len(para) != 2 :
 163             raise FonctionError, 'Erreur de taille dans DEFI_FENE : ' + 'la fenetre exponentielle est definie par exactement deux valeurs'
 164          fene=[1.]*int(para[0])+[numpy.exp(para[1]*(i1-int(para[0]-1))*dt) for i1 in range(int(para[0]),l_ech)]
 165       elif occ[1]['FENETRE'] == 'PART' :
 166          fene=occ[1]['DEFI_FENE']
 167          if len(fene) != l_ech :
 168             raise FonctionError, 'Erreur de taille dans DEFI_FENE : ' + 'La fenetre doit etre definie avec le meme nombre de points que les echantillons'
 169
 170       # normalisation de la fenetre
 171       fene=numpy.divide(fene,numpy.sqrt(numpy.sum(numpy.multiply(fene,fene)))).tolist()
 172
 173    if len(TRANSFERT)+len(INTERSPE) == 0 : #-- on ne rentre rien : interspectre par defaut - fenetre rectangulaire
 174       fene=[1.]*l_ech
 175       INTERSPE=1.;
 176
 177
 178 #--          Recuperation des signaux           --#
 179 #-- Verifications et transformations de Fourier --#
 180 #--         Entrees sous formes de table        --#
 181
 182    tmp=[];
 183    lt=[];
 184    frq=[];
 185    fft=[];
 186    df=[];
 187    num_ord=[];
 188    num_mes=[];
 189
 190
 191    if TAB_ECHANT : # Cas TAB_ECHANT
 192       num_mes_temp= tab_py['NUME_MES'].values()['NUME_MES']
 193       max_mes=numpy.maximum.reduce(num_mes_temp);
 194       num_ord_temp= tab_py['NUME_ORDRE_I'].values()['NUME_ORDRE_I']
 195       long_fonc=[len(fonc_py[i1].VALE.get()) for i1 in range(len(fonc_py))]
 196
 197       N_fen=int(numpy.floor((numpy.minimum.reduce(long_fonc)/2-l_ech)/(l_ech-recouvr))+1)
 198
 199       sig=[];
 200       dt=[];
 201       for i1 in range(len(fonc_py)) :
 202          vale=fonc_py[i1].VALE.get();
 203          temp=(list(vale[0:int(len(vale)/2)]));
 204          sig.append(list(vale[int(len(vale)/2):]));
 205          test_pas=numpy.subtract(temp[1:],temp[0:-1])
 206          crit=test_pas.tolist();
 207          crit.sort();
 208          dt.append(crit[-1]);
 209          if abs((crit[-1]-crit[0])/crit[-1]) > 1.e-5 :
 210             raise FonctionError, 'L'+"'"+'echantillonage doit etre fait a pas constant'
 211
 212       for j1 in range(N_fen) :
 213          for i1 in range(len(fonc_py)) :
 214             fft.append(FFT.fft(numpy.multiply(sig[i1][j1*(l_ech-recouvr):(j1*(l_ech-recouvr)+l_ech)],fene)))
 215             if j1 == 0 : df.append(1./(dt[i1])/l_ech);
 216             num_mes.append(num_mes_temp[i1]+max_mes*j1);
 217             num_ord.append(num_ord_temp[i1]);
 218
 219       if len(df)>1 :
 220          test_df=numpy.subtract(df[1:],df[0:-1])
 221          test_df=test_df.tolist();
 222          test_df.sort();
 223          if abs(test_df[-1]) > 1.e-5 :
 224             raise FonctionError, 'Toutes les fonctions doivent etre definies ' + 'avec la meme frequence d'+"'"+'echantillonage'
 225
 226       frq = [df[-1]*i1 for i1 in range(l_ech)]
 227
 228
 229 #--          Recuperation des signaux           --#
 230 #-- Verifications et transformations de Fourier --#
 231 #--         Entrees sous formes de fonction     --#
 232
 233    if ECHANT:
 234       for occ in l_f :
 235          vale_sig=occ[1]['FONCTION'].Valeurs();
 236          #-- pour les tests ulterieurs --#
 237          lt.append(len(vale_sig[0]))
 238          if len(vale_sig[0]) != len(vale_sig[1]) :
 239             raise FonctionError, 'Les vecteurs associes au temps '+'et aux echantillons doivent etre de meme longueur'
 240          num_mes.append(occ[1]['NUME_MES'])
 241          num_ord.append(occ[1]['NUME_ORDRE_I'])
 242
 243          tmp.append(vale_sig[0])
 244          test_pas=numpy.subtract(vale_sig[0][1:],vale_sig[0][0:-1])
 245          crit=test_pas.tolist();
 246          crit.sort();
 247          if abs((crit[-1]-crit[0])/crit[-1]) > 1.e-5 :
 248             raise FonctionError, 'L'+"'"+'echantillonage doit etre fait a pas constant'
 249          print "vale_sig[1]= ", len(vale_sig[1]), vale_sig[1]
 250          print "  fene = ",len(fene), fene
 251          fft.append(FFT.fft(numpy.multiply(vale_sig[1],fene)))
 252          df.append(1./(crit[-1])/len(vale_sig[0]));
 253
 254
 255       #-- Verification des longueurs --#
 256
 257       test_long=numpy.subtract(lt[1:],lt[0:-1])
 258       test_long=test_long.tolist();
 259       test_long.sort();
 260       if (test_long[-1]-test_long[0]) != 0 :
 261          raise FonctionError, 'Toutes les fonctions doivent etre definies avec le meme nombre de points'
 262
 263       if len(df) > 1 :
 264          test_df=numpy.subtract(df[1:],df[0:-1])
 265          test_df=test_df.tolist();
 266          test_df.sort();
 267          if abs(test_df[-1]) > 1.e-5 :
 268              raise FonctionError, 'Toutes les fonctions doivent etre definies '+'avec la meme frequence d'+"'"+'echantillonage'
 269
 270       frq = [df[-1]*i1 for i1 in range(lt[-1])]
 271
 272
 273 #-- index des numeros d'ordre pour le moyennage
 274
 275    uu=[];
 276    vv=[];
 277    uu=uu+num_ord;
 278    vv=vv+num_ord;
 279    uu.sort();
 280    ind_ord=[];
 281    list_ord=[];
 282    while  len(uu) > 0 :
 283       list_ord.append(uu[0])
 284       tt=[];
 285       for i1 in range(uu.count(uu[0])) :
 286          tt.append(vv.index(uu[0]))
 287          vv[tt[-1]]=0
 288       ind_ord.append(tt)
 289       uu=uu[int(uu.count(uu[0])):]
 290
 291 #-- Calcul de la matrice inter spectrale
 292
 293    if len(INTERSPE) != 0 :
 294       nb_ord = len(list_ord)
 295       dimh   = (nb_ord*(nb_ord+1))/2
 296       l_fc=[];
 297       nume_i1=[]
 298       nume_j1=[]
 299
 300       for i1 in range(nb_ord) :
 301          for j1 in range(i1,nb_ord) :
 302             #-- on ne calcule les spectres que pour des numeros de mesures correspondants
 303             #-- Ca n'a a priori pas de sens de calculer l'interspectre entre deux signaux acquis a des instants differents
 304             #-- Par contre, on peut moyenner deux interspectres obtenus a des instants differents, sous reserve
 305             #-- de stationnarite et d'ergodicite du signal
 306             mes_i1=[num_mes[k1] for k1 in ind_ord[i1]]
 307             mes_j1=[num_mes[k1] for k1 in ind_ord[j1]]
 308             ind_mes=[];
 309             #-- recuperation des indices des fft a prendre en compte pour l'interspectre
 310             for k1 in range(len(mes_i1)) :
 311                if mes_i1[k1] in mes_j1 :
 312                   ind_mes.append([ind_ord[i1][k1],
 313                                   ind_ord[j1][mes_j1.index(mes_i1[k1])]])
 314
 315             #-- Calcul des interspectres
 316             dsp=[0.j]*l_ech;
 317             if len(ind_mes) > 0 :
 318                for l1 in range(len(ind_mes)) :
 319                   dsp_t=numpy.multiply(numpy.conjugate(fft[ind_mes[l1][0]]),fft[ind_mes[l1][1]])
 320                   dsp_t=numpy.divide(dsp_t,l_ech*len(ind_mes))
 321                   dsp=numpy.add(dsp,dsp_t)
 322                dsp=dsp.tolist();
 323                dsp_r=[];
 324
 325                for k1 in range(int(numpy.floor(l_ech/2))) :
 326                   dsp_r=dsp_r+[frq[k1],dsp[k1].real,dsp[k1].imag]
 327
 328                _fonc = DEFI_FONCTION(NOM_PARA='FREQ',VALE_C=dsp_r,);
 329                l_fc.append(_fonc.nom)
 330                nume_i1.append(list_ord[i1])
 331                nume_j1.append(list_ord[j1])
 332
 333       mcfact=[]
 334       mcfact.append(_F(PARA='NOM_CHAM'    ,LISTE_K='DSP'   ,NUME_LIGN=(1,) ))
 335       mcfact.append(_F(PARA='OPTION'      ,LISTE_K='TOUT'  ,NUME_LIGN=(1,) ))
 336       mcfact.append(_F(PARA='DIMENSION'   ,LISTE_I=(dimh,) ,NUME_LIGN=(1,) ))
 337       mcfact.append(_F(PARA='NUME_ORDRE_I',LISTE_I=nume_i1 ,NUME_LIGN=range(2,dimh+2) ))
 338       mcfact.append(_F(PARA='NUME_ORDRE_J',LISTE_I=nume_j1 ,NUME_LIGN=range(2,dimh+2) ))
 339       mcfact.append(_F(PARA='FONCTION_C'  ,LISTE_K=l_fc    ,NUME_LIGN=range(2,dimh+2)))
 340       self.DeclareOut('tab_inte',self.sd)
 341       tab_inte=CREA_TABLE(LISTE=mcfact,
 342                           TITRE='',
 343                           TYPE_TABLE='TABLE_FONCTION')
 344
 345
 346
 347
 348 #-- Calcul des transferts
 349
 350    if len(TRANSFERT) != 0 :
 351
 352       l_fc=[];
 353       nume_i1=[]
 354       nume_j1=[]
 355
 356       #-- test sur les entrees pour les references --#
 357       if type(l_H[0][1]['REFER'])==int :
 358          refer=[];
 359          refer.append(l_H[0][1]['REFER'])
 360       elif type(l_H[0][1]['REFER'])==tuple :
 361          refer=list(l_H[0][1]['REFER'])
 362
 363       ind_refer=[];
 364       dimh   = len(refer)*(len(list_ord)-len(refer))
 365       for k1 in range(len(refer)) :
 366          for l1 in range(len(list_ord)) :
 367             if refer[k1] == list_ord[l1] : ind_refer.append(l1);
 368
 369       #-- H1 : interspectre / autospectre
 370       #-- H2 : autospectre / interspectre
 371       #-- CO : coherence entre H1 et H2.
 372
 373       if l_H[0][1]['ESTIM']!='HV' :
 374          for i1 in range(len(refer)) :
 375             for j1 in range(len(list_ord)) :
 376                if refer[i1] != list_ord[j1] :
 377                   mes_i1=[num_mes[k1] for k1 in ind_ord[ind_refer[i1]]]  #-- mesures des efforts
 378                   mes_j1=[num_mes[k1] for k1 in ind_ord[j1]]  #-- mesures des reponses
 379
 380                   ind_mes=[];
 381                   #-- recuperation des indices des mesures a predre en compte pour les spectres
 382                   for k1 in range(len(mes_i1)) :
 383                      if mes_i1[k1] in mes_j1 :
 384                         ind_ord[j1][mes_j1.index(mes_i1[k1])]
 385                         ind_mes.append([ind_ord[ind_refer[i1]][k1],ind_ord[j1][mes_j1.index(mes_i1[k1])]])
 386
 387                   #-- Calcul des FRF
 388                   Guu=[0.j]*l_ech;
 389                   Gyy=[0.j]*l_ech;
 390                   Gyu=[0.j]*l_ech;
 391                   if len(ind_mes) > 0 :
 392                      for l1 in range(len(ind_mes)) :
 393                         Guu_t=numpy.multiply(numpy.conjugate(fft[ind_mes[l1][0]]),fft[ind_mes[l1][0]])
 394                         Guu=numpy.add(Guu,Guu_t)
 395                         Gyu_t=numpy.multiply(numpy.conjugate(fft[ind_mes[l1][1]]),fft[ind_mes[l1][0]])
 396                         Gyu=numpy.add(Gyu,Gyu_t)
 397                         Gyy_t=numpy.multiply(numpy.conjugate(fft[ind_mes[l1][1]]),fft[ind_mes[l1][1]])
 398                         Gyy=numpy.add(Gyy,Gyy_t)
 399
 400                      if l_H[0][1]['ESTIM']=='H1' :
 401                         frf=numpy.divide(numpy.conjugate(Gyu),Guu);
 402                         nom_frf='FRF-H1';
 403                      elif l_H[0][1]['ESTIM']=='H2' :
 404                         frf=numpy.divide(Gyy,Gyu);
 405                         nom_frf='FRF-H2';
 406                      elif l_H[0][1]['ESTIM']=='CO' :
 407                         H1=numpy.divide(numpy.conjugate(Gyu),Guu);
 408                         H2=numpy.divide(Gyy,Gyu);
 409                         frf=numpy.divide(H1,H2);
 410                         nom_frf='FRF-COH';
 411
 412                      frf=frf.tolist();
 413                      frf_r=[];
 414
 415                      for k1 in range(int(numpy.floor(l_ech/2))) :
 416                         frf_r=frf_r+[frq[k1],frf[k1].real,frf[k1].imag]
 417
 418                      _fonc = DEFI_FONCTION(NOM_PARA='FREQ',VALE_C=frf_r,);
 419                      l_fc.append(_fonc.nom)
 420                      nume_i1.append(refer[i1])
 421                      nume_j1.append(list_ord[j1])
 422
 423       #-- On remplit la table_fonction avec tout ce qui va bien
 424
 425       mcfact=[]
 426       mcfact.append(_F(PARA='NOM_CHAM'    ,LISTE_K=nom_frf ))
 427       mcfact.append(_F(PARA='OPTION'      ,LISTE_K='TOUT'  ))
 428       mcfact.append(_F(PARA='DIMENSION'   ,LISTE_I=(dimh,) ))
 429       mcfact.append(_F(PARA='NUME_ORDRE_I',LISTE_I=nume_i1 ))
 430       mcfact.append(_F(PARA='NUME_ORDRE_J',LISTE_I=nume_j1 ))
 431       mcfact.append(_F(PARA='FONCTION_C'  ,LISTE_K=l_fc  ))
 432       self.DeclareOut('tab_inte',self.sd)
 433       tab_inte=CREA_TABLE(LISTE=mcfact,
 434                           TITRE='',
 435                           TYPE_TABLE='TABLE_FONCTION')
 436