Aster/Cata/cataSTA9/Macro/dyna_iss_vari_ops.py

   1 #@ MODIF dyna_iss_vari_ops Macro  DATE 21/04/2008   AUTEUR ZENTNER I.ZENTNER
   2 # -*- coding: iso-8859-1 -*-
   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 from Accas import _F
  21 import string
  22
  23 def dyna_iss_vari_ops(self, NOM_CMP, PRECISION, INTERF,MATR_COHE, FREQ_INIT,UNITE_RESU_FORC,
  24                        NB_FREQ, PAS, UNITE_RESU_IMPE, TYPE, MATR_GENE ,INFO,
  25                          **args):
  26    """
  27       Macro DYNA_ISS_VARI
  28    """
  29    ier=0
  30    import Numeric as Num
  31    import LinearAlgebra as LinAl
  32    import MLab
  33    import os
  34    import aster
  35    diag = MLab.diag
  36    max = MLab.max
  37    min = MLab.min
  38    sum = Num.sum
  39    abs = Num.absolute
  40    conj = Num.conjugate
  41    from Utilitai.Table import Table
  42    from Utilitai.Utmess import  UTMESS
  43
  44    def get_group_coord(group):
  45       """Retourne les coordonnees des noeuds du groupe 'group'
  46       """
  47       l_ind = Num.array(coll_grno.get('%-8s' % group, [])) - 1
  48       return Num.take(t_coordo, l_ind)
  49
  50
  51    # On importe les definitions des commandes a utiliser dans la macro
  52
  53    COMB_MATR_ASSE = self.get_cmd('COMB_MATR_ASSE')
  54    LIRE_IMPE_MISS = self.get_cmd('LIRE_IMPE_MISS')
  55    LIRE_FORC_MISS = self.get_cmd('LIRE_FORC_MISS')
  56    COMB_MATR_ASSE = self.get_cmd('COMB_MATR_ASSE')
  57
  58    CREA_CHAMP = self.get_cmd('CREA_CHAMP')
  59    DYNA_LINE_HARM = self.get_cmd('DYNA_LINE_HARM')
  60    DETRUIRE= self.get_cmd('DETRUIRE')
  61
  62    DEFI_FONCTION  = self.get_cmd('DEFI_FONCTION')
  63    CREA_TABLE     = self.get_cmd('CREA_TABLE')
  64
  65    # Comptage commandes + declaration concept sortant
  66    self.set_icmd(1)
  67    self.DeclareOut('tab_out', self.sd)
  68    macro='DYNA_ISS_VARI'
  69 #--------------------------------------------------------
  70    dgene = MATR_GENE[0].cree_dict_valeurs(MATR_GENE[0].mc_liste)
  71    if dgene['MATR_AMOR'] != None:
  72      aster.affiche('MESSAGE',' MATR_AMOR existe')
  73      __ma_amort = MATR_GENE['MATR_AMOR']
  74    else:
  75      __ma_amort=COMB_MATR_ASSE(CALC_AMOR_GENE=_F(MASS_GENE = MATR_GENE['MATR_MASS'] ,
  76                                         RIGI_GENE = MATR_GENE['MATR_RIGI'] ,
  77                                         AMOR_REDUIT= (  0.0,),
  78                                          ),
  79                                   );
  80      aster.affiche('MESSAGE',' MATR_AMOR pas donnee, on prend AMOR_REDUIT=0.0,')
  81 #   dint = INTERF[0].cree_dict_valeurs(INTERF[0].mc_liste)
  82 #   dcoh = MATR_COHE[0].cree_dict_valeurs(MATR_COHE[0].mc_liste)
  83
  84    from SD.sd_maillage import sd_maillage
  85    from SD.sd_base_modale import sd_base_modale
  86    from SD.sd_resultat import sd_resultat
  87    from SD.sd_cham_gene import sd_cham_gene
  88    # MAILLAGE
  89    nom_bamo = MATR_GENE['MATR_RIGI'].REFA.get()[0]
  90    nume_ddl = aster.getvectjev(nom_bamo[0:8] + '           .REFD        ' )[3]
  91    nom_mail = aster.getvectjev( nume_ddl[0:19] + '.REFN        ' )[0]
  92    num_mail = sd_maillage(nom_mail)
  93    # MODELE, DDLGENE
  94    nom_ddlgene = MATR_GENE['MATR_RIGI'].REFA.get()[1]
  95    nom_modele = aster.getvectjev( nume_ddl[0:19] + '.LILI        ' )[1]
  96    nume_resu = self.jdc.sds_dict[string.strip(nom_bamo)]
  97    nume_ddlgene = self.jdc.sds_dict[string.strip(nom_ddlgene)]
  98    nume_modele = self.jdc.sds_dict[string.strip(nom_modele[0:8])]
  99
 100    #TEST base modale
 101    nom_bamo1 = MATR_GENE['MATR_MASS'].REFA.get()[0]
 102    nom_bamo2 = MATR_GENE['MATR_RIGI'].REFA.get()[0]
 103    if string.strip(nom_bamo) != string.strip(nom_bamo1) or string.strip(nom_bamo) != string.strip(nom_bamo2) or string.strip(nom_bamo1) != string.strip(nom_bamo2):
 104       UTMESS('F','ALGORITH5_42')
 105
 106
 107    nbnot, nbl, nbma, nbsm, nbsmx, dime = num_mail.DIME.get()
 108
 109    # coordonnees des noeuds
 110    l_coordo = num_mail.COORDO.VALE.get()
 111    t_coordo = Num.array(l_coordo)
 112    t_coordo.shape = nbnot, 3
 113    # groupes de noeuds
 114    coll_grno = num_mail.GROUPENO.get()
 115    GROUP_NO_INTER=INTERF['GROUP_NO_INTERF']
 116    noe_interf = get_group_coord(GROUP_NO_INTER)
 117    #  print noe_interf
 118    nbno, nbval = noe_interf.shape
 119    if INFO==2:
 120       aster.affiche('MESSAGE','NBNO INTERFACE : '+str(nbno))
 121   # MODES
 122    nbval, nbmodt,nbmodd,nbmods = nume_resu.UTIL.get()
 123
 124
 125    nbmodt2 = MATR_GENE['MATR_RIGI'].DESC.get()[1]
 126    if nbmodt2 != nbmodt:
 127        UTMESS('F','ALGORITH5_42')
 128
 129    if INFO==2:
 130       texte = 'NOMBRE DE MODES: '+str(nbmodt)+'   MODES DYNAMIQUES: '+str(nbmodd)+'   MODES STATIQUES: '+str(nbmods)
 131       aster.affiche('MESSAGE',texte)
 132       aster.affiche('MESSAGE','COMPOSANTE '+NOM_CMP)
 133    SPEC = Num.zeros((NB_FREQ,nbmodt,nbmodt), Num.Float)+1j
 134 #
 135 #---------------------------------------------------------------------
 136   # BOUCLE SUR LES FREQUENCES
 137    VITE_ONDE = MATR_COHE['VITE_ONDE']
 138    alpha = MATR_COHE['PARA_ALPHA']
 139    abscisse = [None]*NB_FREQ
 140
 141    for k in range(0,NB_FREQ):
 142       freqk=FREQ_INIT+PAS*k
 143       aster.affiche('MESSAGE','FREQUENCE DE CALCUL: '+str(freqk))
 144
 145       # Matrice de coherence
 146       XX=noe_interf[:,0]
 147       YY=noe_interf[:,1]
 148
 149       XN=Num.repeat(XX,nbno)
 150       YN=Num.repeat(YY,nbno)
 151       XR=Num.reshape(XN,(nbno,nbno))
 152       YR=Num.reshape(YN,(nbno,nbno))
 153       XRT=Num.transpose(XR)
 154       YRT=Num.transpose(YR)
 155       DX=XR-XRT
 156       DY=YR-YRT
 157       DIST=DX**2+DY**2
 158       COHE=Num.exp(-(DIST*(alpha*freqk/VITE_ONDE)**2.))
 159
 160       # On desactive temporairement les FPE qui pourraient etre generees (a tord!) par blas
 161       aster.matfpe(-1)
 162       eig, vec =LinAl.eigenvectors(COHE)
 163       aster.matfpe(1)
 164       eig=eig.real
 165       vec=vec.real
 166       # on rearrange selon un ordre decroissant
 167       eig = Num.where(eig < 1.E-10, 0.0, eig)
 168       order = (Num.argsort(eig)[::-1])
 169       eig = Num.take(eig, order)
 170       vec = Num.take(vec, order, 0)
 171
 172       #-----------------------
 173       # Nombre de modes POD a retenir
 174       etot=sum(diag(COHE))
 175       ener=0.0
 176       nbme=0
 177
 178       if INFO==2:
 179          aster.affiche('MESSAGE','ETOT :'+str(etot))
 180       while nbme < nbno:
 181          ener= eig[nbme]+ener
 182          prec=ener/etot
 183          nbme=nbme+1
 184          if INFO==2:
 185             aster.affiche('MESSAGE','VALEUR PROPRE  '+str(nbme)+' : '+str(eig[nbme-1]))
 186          if prec > PRECISION :
 187             break
 188
 189       aster.affiche('MESSAGE','NOMBRE DE MODES POD RETENUS : '+str(nbme))
 190       aster.affiche('MESSAGE','PRECISION (ENERGIE RETENUE) : '+str(prec))
 191
 192       PVEC=Num.zeros((nbme,nbno), Num.Float)
 193       for k1 in range(0,nbme):
 194          PVEC[k1, 0:nbno]=Num.sqrt(eig[k1])*vec[k1]
 195       # CALCUL DE FS variable-------------------------------
 196       XO=Num.zeros((nbme,nbmods), Num.Float)
 197       if NOM_CMP=='DX':
 198          COMP = 1
 199       elif NOM_CMP=='DY':
 200          COMP = 2
 201       elif NOM_CMP=='DZ':
 202          COMP = 3
 203
 204    #---------MODES interface
 205       # ----- boucle sur les modes statiques
 206       for mods in range(0,nbmods):
 207          nmo = nbmodd+mods+1
 208          __CHAM=CREA_CHAMP( TYPE_CHAM='NOEU_DEPL_R',
 209                 OPERATION='EXTR',
 210                 NUME_ORDRE=nmo,
 211                 RESULTAT = nume_resu  ,
 212                 NOM_CHAM = 'DEPL'
 213                       );
 214          MCMP =__CHAM.EXTR_COMP(NOM_CMP,[GROUP_NO_INTER]).valeurs
 215
 216          NNO =__CHAM.EXTR_COMP(NOM_CMP,[GROUP_NO_INTER], topo=1).noeud
 217
 218
 219          som=sum(MCMP)
 220          max1=max(MCMP)
 221          min1=min(MCMP)
 222          maxm=max([abs(max1),abs(min1)])
 223       #CALCUL DE XO
 224 #  on recupere la composante COMP (dx,dy,dz) des modes et on projete
 225          #  CAS 1: MODES DE CORPS RIGIDE
 226          if INTERF['MODE_INTERF'] =='CORP_RIGI':
 227             for modp in range(0,nbme):
 228                #modes de translation
 229                if mods+1 <=3:
 230                   if abs(som)<10.E-6:
 231                      XO[modp,mods]=0.0
 232                   else :
 233                      fact=1./som
 234                      XO[modp,mods]=fact*Num.innerproduct(MCMP,PVEC[modp])
 235                #modes de rotation
 236                else:
 237                   if maxm<10.E-6:
 238                      if som<10.E-6:
 239                         XO[modp,mods]=0.0
 240                      else :
 241                         UTMESS('F','ALGORITH6_86')
 242                   else :
 243                      fact = 1./(nbno)
 244                      XO[modp,mods]=1./(maxm**2.)*fact*Num.innerproduct(MCMP,PVEC[modp])
 245
 246          # CAS 2: MODES EF
 247          if INTERF['MODE_INTERF'] =='TOUT':
 248             for modp in range(0,nbme):
 249                if abs(som)<10.E-6:
 250                   if maxm<10.E-6:
 251                      XO[modp,mods]=0.0
 252                   else:
 253                      UTMESS('F','UTILITAI5_89')
 254                else:
 255                   fact=1./som
 256                   XO[modp,mods]=fact*Num.innerproduct(MCMP,PVEC[modp])
 257
 258          DETRUIRE(CONCEPT=_F(NOM=(__CHAM)),INFO=1)
 259
 260    #----Impedances etc.-----------------------------------------------------------------
 261
 262       if k>0:
 263          DETRUIRE(CONCEPT=_F(NOM=(__impe,__fosi,__rito)),INFO=1)
 264
 265       __impe = LIRE_IMPE_MISS(BASE=nume_resu,
 266                            TYPE=TYPE,
 267                            NUME_DDL_GENE=nume_ddlgene,
 268                            UNITE_RESU_IMPE= UNITE_RESU_IMPE,
 269                            FREQ_EXTR=freqk,
 270                            );
 271       __rito=COMB_MATR_ASSE(COMB_C=(
 272                                 _F(MATR_ASSE=__impe,
 273                                  COEF_C=1.0+0.j,),
 274                                 _F(MATR_ASSE=MATR_GENE['MATR_RIGI'],
 275                                  COEF_C=1.0+0.j,),
 276                                  ),
 277                                  SANS_CMP='LAGR',
 278                                  );
 279       __fosi = LIRE_FORC_MISS(BASE=nume_resu,
 280                            NUME_DDL_GENE=nume_ddlgene,
 281                            NOM_CMP=NOM_CMP,
 282                            NOM_CHAM='DEPL',
 283                            UNITE_RESU_FORC = UNITE_RESU_FORC,
 284                            FREQ_EXTR=freqk,);
 285       # impedance
 286       MIMPE=__impe.EXTR_MATR_GENE()
 287       #  extraction de la partie modes interface
 288       KRS = MIMPE[nbmodd:nbmodt,nbmodd:nbmodt]
 289
 290       # force sismique pour verif
 291 #      FS0=__fosi.EXTR_VECT_GENE_C()
 292 #      FSE=FS0[nbmodd:nbmodt][:]
 293       SP=Num.zeros((nbmodt,nbmodt),Num.Float)
 294       for k1 in range(0,nbme):
 295          #  calcul de la force sismique mode POD par mode POD
 296          FS = Num.matrixmultiply(KRS,XO[k1])
 297          Fzero=Num.zeros((1,nbmodd),Num.Float)
 298          FS2=Num.concatenate((Fzero,Num.reshape(FS,(1,nbmods))),1)
 299       #  Calcul harmonique
 300          __fosi.RECU_VECT_GENE_C(FS2[0])
 301          __dyge = DYNA_LINE_HARM(MODELE= nume_modele,
 302                           MATR_MASS = MATR_GENE['MATR_MASS'],
 303                           MATR_RIGI = __rito,
 304                           FREQ = freqk,
 305                           MATR_AMOR = __ma_amort,
 306                           EXCIT =_F ( VECT_ASSE = __fosi,
 307                                       COEF_MULT= 1.0,
 308                                   ),
 309                         );
 310          #  recuperer le vecteur modal depl calcule par dyge
 311          desc = __dyge.DESC.get()
 312          assert desc[0].strip() == 'DEPL', 'Champ DEPL non trouvé'
 313          nomcham = __dyge.TACH.get()[1][0].strip()
 314          cham = sd_cham_gene(nomcham)
 315          RS = Num.array(cham.VALE.get())
 316          SP=SP+RS*conj(RS[:,Num.NewAxis])
 317          DETRUIRE(CONCEPT=_F(NOM=(__dyge)),INFO=1)
 318
 319
 320       SPEC[k]=SP
 321
 322       abscisse[k]= freqk
 323 ##---------------------------------------------------------------------
 324 #  Ecriture des tables
 325 #---------------------------------------------------------------------
 326 #   ------ CREATION DE L OBJET TABLE
 327    tab = Table()
 328    tab.append({'NOM_CHAM' : 'DSP', 'OPTION' : 'TOUT',  'DIMENSION' : nbmodt})
 329    foncc=Num.array([None]*NB_FREQ*3)
 330    for k2 in range(nbmodt):
 331       for k1 in range(k2+1):
 332          ks=0
 333          for k in range(NB_FREQ) :
 334             foncc[ks]=abscisse[k]
 335             foncc[ks+1]= SPEC[k][k1,k2].real
 336             foncc[ks+2]= SPEC[k][k1,k2].imag
 337             ks=ks+3
 338          _f = DEFI_FONCTION(NOM_PARA='FREQ',
 339                          NOM_RESU='SPEC',
 340                          VALE_C  = foncc.tolist() )
 341
 342       # Ajout d'une ligne dans la Table
 343          tab.append({'NUME_ORDRE_I' : k1+1, 'NUME_ORDRE_J' : k2+1, 'FONCTION_C' : _f.nom})
 344
 345
 346    # Creation du concept en sortie
 347    tab_out = CREA_TABLE(TYPE_TABLE='TABLE_FONCTION',
 348                         **tab.dict_CREA_TABLE())
 349    return ier