Aster/Cata/cataSTA10/Macro/dyna_iss_vari_ops.py

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