Aster/Cata/cataSTA10/Macro/observation_ops.py

   1 #@ MODIF observation_ops Macro  DATE 13/01/2011   AUTEUR PELLET J.PELLET
   2
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2011  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
  22 # RESPONSABLE BODEL C.BODEL
  23
  24 def observation_ops(self,
  25                     PROJECTION   = None,
  26                     MODELE_1     = None,
  27                     MODELE_2     = None,
  28                     RESULTAT     = None,
  29                     MATR_A       = None,
  30                     MATR_B       = None,
  31                     MODI_REPERE  = None,
  32                     NOM_CHAM     = None,
  33                     FILTRE       = None,
  34                     EPSI_MOYENNE = None,
  35                     **args):
  36
  37     """
  38      Ecriture de la macro MACRO_OBSERVATION
  39     """
  40     ier=0
  41
  42
  43     # La macro compte pour 1 dans la numerotation des commandes
  44     self.set_icmd(1)
  45
  46     # on transforme le mc MODI_REPERE pour ne pas le confondre avec l'operateur
  47     # du meme nom
  48     MODIF_REPERE = MODI_REPERE
  49
  50     # importation de commandes
  51     import aster
  52     from Accas import _F
  53     from Utilitai.Utmess     import UTMESS
  54     from Cata.cata import mode_meca, dyna_harmo, evol_elas,dyna_trans
  55     MODI_REPERE = self.get_cmd('MODI_REPERE')
  56     PROJ_CHAMP  = self.get_cmd('PROJ_CHAMP')
  57     CREA_CHAMP  = self.get_cmd('CREA_CHAMP')
  58     CREA_RESU   = self.get_cmd('CREA_RESU')
  59     POST_RELEVE_T   = self.get_cmd('POST_RELEVE_T')
  60     DETRUIRE = self.get_cmd('DETRUIRE')
  61
  62     # dans **args, on range les options de PROJ_CHAMP facultatives, et dont on
  63     # ne sert pas par la suite
  64     mcfact = args
  65
  66     if not isinstance(NOM_CHAM,tuple):
  67         NOM_CHAM = [NOM_CHAM]
  68
  69     TYPE_CHAM = None
  70     RESU = None
  71
  72     if isinstance( RESULTAT, evol_elas): TYPE_RESU='EVOL_ELAS'
  73     if isinstance( RESULTAT, dyna_trans): TYPE_RESU='DYNA_TRANS'
  74     if isinstance( RESULTAT, dyna_harmo): TYPE_RESU='DYNA_HARMO'
  75     if isinstance( RESULTAT, mode_meca): TYPE_RESU='MODE_MECA'
  76
  77     self.DeclareOut( 'RESU', self.sd)
  78
  79     # recuperation du maillage associe au modele numerique
  80     _maillag = aster.getvectjev( MODELE_1.nom.ljust(8) + '.MODELE    .LGRF' )
  81     maillage = _maillag[0].strip()
  82     mayanum = self.get_concept(maillage)
  83
  84     # modele numerique 2D ou 3D
  85     typmod= mayanum.DIME.get()
  86     typmod = typmod[5]
  87
  88     # recuperation du maillage associe au modele experimental
  89     _maillag = aster.getvectjev( MODELE_2.nom.ljust(8) + '.MODELE    .LGRF' )
  90     maillage = _maillag[0].strip()
  91     mayaexp = self.get_concept(maillage)
  92
  93     # cham_mater et cara_elem pour le resultat a projeter
  94     iret,ibid,nom_cara_elem = aster.dismoi('F','CARA_ELEM',RESULTAT.nom,'RESULTAT')
  95     if len(nom_cara_elem) > 0    :
  96         assert nom_cara_elem.strip() != "#PLUSIEURS" , nom_cara_elem
  97         if nom_cara_elem.strip() == "#AUCUN" :
  98             cara_elem = None
  99         else :
 100             cara_elem = self.get_concept(nom_cara_elem.strip())
 101     else:
 102         cara_elem = None
 103
 104     iret,ibid,nom_cham_mater = aster.dismoi('F','CHAM_MATER',RESULTAT.nom,'RESULTAT')
 105     if len(nom_cham_mater) > 0 :
 106         assert nom_cham_mater.strip() != "#PLUSIEURS" , nom_cham_mater
 107         if nom_cham_mater.strip() == "#AUCUN" :
 108             cham_mater = None
 109         else :
 110             cham_mater = self.get_concept(nom_cham_mater.strip())
 111     else:
 112         cham_mater = None
 113
 114     # afreq pour les frequences propres
 115     if isinstance( RESULTAT, mode_meca):
 116         # frequences propres
 117         from Cata.cata import RECU_TABLE
 118         __freq  = RECU_TABLE(CO=RESULTAT,
 119                              NOM_PARA='FREQ',);
 120         afreq  = __freq.EXTR_TABLE().Array('NUME_ORDRE','FREQ')
 121         # noms des matrices
 122         if MATR_A !=None or MATR_B !=None:
 123             # recherche du nume_ddl associe
 124             iret,ibid,nom_nume_ddl = aster.dismoi('F','NOM_NUME_DDL',MATR_A.nom,'MATR_ASSE')
 125             NUME_DDL = self.get_concept(nom_nume_ddl)
 126             # coherence avec le nom associe a MODELE_2 :
 127             iret,ibid,nom_modele = aster.dismoi('F','NOM_MODELE',nom_nume_ddl,'NUME_DDL')
 128             if nom_modele.strip() != MODELE_2.nom.strip():
 129                 UTMESS('F','CALCESSAI0_10')
 130         else:
 131             UTMESS('A','CALCESSAI0_9')
 132             NUME_DDL = None
 133
 134     else:
 135         afreq = None
 136         NUME_DDL = None
 137
 138     indice = range(len(RESULTAT.LIST_VARI_ACCES()['NUME_ORDRE']))
 139
 140 #***********************************************
 141 #  PHASE DE CALCUL DE LA DEFORMATION MOYENNE AUX NOEUDS
 142 #  CHAMP CALCULE SUR LE MODELE NUMERIQUE
 143 #***********************************************
 144
 145     resu_epsi = None
 146     if EPSI_MOYENNE != None :
 147       for nomcham in NOM_CHAM:
 148         if nomcham == 'EPSI_NOEU':
 149           if isinstance( RESULTAT, dyna_harmo):
 150             TYPE_CHAM  = 'NOEU_EPSI_C'
 151           else:
 152             TYPE_CHAM  = 'NOEU_EPSI_R'
 153
 154       if TYPE_CHAM == None:
 155             UTMESS('F', 'UTILITAI8_24', valk = ['NOEU_EPSI', nomcham])
 156       else:
 157         num_ordr = RESULTAT.LIST_VARI_ACCES()['NUME_ORDRE']
 158
 159         if isinstance( RESULTAT, evol_elas):
 160             list_inst = RESULTAT.LIST_VARI_ACCES()['INST']
 161         if isinstance( RESULTAT, dyna_trans):
 162             list_inst = RESULTAT.LIST_VARI_ACCES()['INST']
 163         if isinstance( RESULTAT, dyna_harmo):
 164             list_freq = RESULTAT.LIST_VARI_ACCES()['FREQ']
 165
 166         liste = []
 167
 168       # il faut calculer le champ complet
 169         if typmod == 2:
 170             nom_cmps = ['EPXX','EPYY','EPZZ','EPXY',]
 171         else:
 172             nom_cmps = ['EPXX','EPYY','EPZZ','EPXY','EPXZ','EPYZ',]
 173
 174         argsi = {'ACTION' : [], }
 175         lnoeuds = { }
 176         nb_mcfact = 0
 177         seuil = []
 178         masque = []
 179         for epsi_moye in EPSI_MOYENNE :
 180             mcfactr = { }
 181             mcfacti = { }
 182             l_noeud = None
 183             val_masque = []
 184             seuil_lu = epsi_moye['SEUIL_VARI']
 185             if type(seuil_lu) == tuple :
 186                 val_seuil = seuil_lu[0]
 187             else:
 188                 val_seuil = seuil_lu
 189             seuil.append(val_seuil)
 190             masque_lu = epsi_moye['MASQUE']
 191             if type(masque_lu) != tuple :
 192                 val_masque.append(masque_lu)
 193             else:
 194                 val_masque = masque_lu
 195             masque.append(val_masque)
 196             for typ in ['NOEUD','GROUP_NO','MAILLE','GROUP_MA']:
 197               if epsi_moye[typ] != None:
 198                 l_noeud = find_no(mayanum, {typ : epsi_moye[typ]})
 199                 nb_mcfact = nb_mcfact + 1
 200                 for i in range(len(l_noeud)):
 201                     l_noeud[i]=l_noeud[i].strip()
 202                 lnoeuds[str(nb_mcfact)]=l_noeud
 203
 204             if l_noeud == None:
 205                 UTMESS('F','MODELISA3_13',valk=['EPSI_MOYENNE'])
 206
 207             if TYPE_CHAM[-1:] == 'C':
 208               mcfactr = { 'NOM_CMP'   : nom_cmps,
 209                           'OPERATION' : 'EXTRACTION',
 210                           'INTITULE' : str('R'+str(nb_mcfact)),
 211                           'FORMAT_C' : 'REEL',
 212                           'NOEUD' : l_noeud,
 213                           'NOM_CHAM'  : 'EPSI_NOEU',
 214                           'RESULTAT'  : RESULTAT,
 215                           'NUME_ORDRE'  : num_ordr,
 216                         }
 217               argsi['ACTION'].append(mcfactr)
 218               mcfacti = { 'NOM_CMP'   : nom_cmps,
 219                           'OPERATION' : 'EXTRACTION',
 220                           'INTITULE' : str('I'+str(nb_mcfact)),
 221                           'FORMAT_C' : 'IMAG',
 222                           'NOEUD' : l_noeud,
 223                           'NOM_CHAM'  : 'EPSI_NOEU',
 224                           'RESULTAT'  : RESULTAT,
 225                           'NUME_ORDRE'  : num_ordr,
 226                          }
 227               argsi['ACTION'].append(mcfacti)
 228             else:
 229               mcfactr = { 'NOM_CMP'   : nom_cmps,
 230                           'OPERATION' : 'EXTRACTION',
 231                           'INTITULE' : str(nb_mcfact),
 232                           'NOEUD' : l_noeud,
 233                           'NOM_CHAM'  : 'EPSI_NOEU',
 234                           'RESULTAT'  : RESULTAT,
 235                           'NUME_ORDRE'  : num_ordr,
 236                          }
 237               argsi['ACTION'].append(mcfactr)
 238
 239         _tepsi=POST_RELEVE_T(
 240                    **argsi)
 241
 242         table=_tepsi.EXTR_TABLE()
 243
 244         DETRUIRE( CONCEPT= _F( NOM = _tepsi ), INFO=1)
 245
 246         mcfact2 = { }
 247         __chame = [None]*len(indice)
 248         for ind in indice:
 249          argsa = {'AFFE' : [], }
 250          for mcfacta in range(nb_mcfact):
 251           l_noeud_mcfact = lnoeuds[str(mcfacta+1)]
 252           l_vmoye=[]
 253           l_cmp_vari=[]
 254           seuil_mc=seuil[mcfacta]
 255           masque_mc=masque[mcfacta]
 256           for cmp in nom_cmps:
 257            lur = 0
 258            lui = 0
 259            l_valr= []
 260            l_vali= []
 261            l_valc= []
 262            l_val= []
 263            for row in table.rows:
 264             if TYPE_CHAM[-1:] == 'C':
 265              if row['INTITULE'].strip() == str('R'+str(mcfacta+1))   \
 266                       and row['NUME_ORDRE'] == num_ordr[ind] :
 267                 l_valr.append(row[cmp])
 268                 lur = 1
 269              elif row['INTITULE'].strip() == str('I'+str(mcfacta+1))   \
 270                       and row['NUME_ORDRE'] == num_ordr[ind] :
 271                 l_vali.append(row[cmp])
 272                 lui = 1
 273
 274             else:
 275              if row['INTITULE'].strip() == str(mcfacta+1)   \
 276                       and row['NUME_ORDRE'] == num_ordr[ind]:
 277                   l_val.append(row[cmp])
 278
 279            if TYPE_CHAM[-1:] == 'C':
 280              if lur and lui :
 281                 if len(l_valr) != len(l_vali):
 282                   UTMESS('F','POSTRELE_59')
 283                 for i in range(len(l_valr)):
 284                   l_valc.append(complex(l_valr[i],l_vali[i]))
 285
 286                 lur = 0
 287                 lui = 0
 288              else:
 289                 UTMESS('F','POSTRELE_59')
 290
 291           # on regarde a la fois la partie reelle et la partie imag pour les complexes
 292              vmoyer = sum(l_valr)/len(l_valr)
 293              vmoyei = sum(l_vali)/len(l_vali)
 294              vmoye = sum(l_valc)/len(l_valc)
 295              vminr = min(l_valr)
 296              vmini = min(l_vali)
 297              vmaxr = max(l_valr)
 298              vmaxi = max(l_vali)
 299              if vmoyer > 0:
 300                if (vmaxr > vmoyer*(1.+seuil_mc)) or (vminr < vmoyer*(1-seuil_mc)):
 301                  l_cmp_vari.append(cmp)
 302              if vmoyei > 0:
 303                if (vmaxi > vmoyei*(1.+seuil_mc)) or (vmini < vmoyei*(1-seuil_mc)):
 304                  l_cmp_vari.append(cmp)
 305              if vmoyer < 0:
 306                if (vminr > vmoyer*(1.-seuil_mc)) or (vmaxr < vmoyer*(1+seuil_mc)):
 307                  l_cmp_vari.append(cmp)
 308              if vmoyei < 0:
 309                if (vmini > vmoyei*(1.-seuil_mc)) or (vmaxi < vmoyei*(1+seuil_mc)):
 310                  l_cmp_vari.append(cmp)
 311            else:
 312              vmoye = sum(l_val)/len(l_val)
 313              vmin = min(l_val)
 314              vmax = max(l_val)
 315              if vmoye > 0:
 316                if (vmax > vmoye*(1.+seuil_mc)) or (vmin < vmoye*(1-seuil_mc)):
 317                  l_cmp_vari.append(cmp)
 318              if vmoye < 0:
 319                if (vmin > vmoye*(1.-seuil_mc)) or (vmax < vmoye*(1+seuil_mc)):
 320                  l_cmp_vari.append(cmp)
 321
 322            l_vmoye.append(vmoye)
 323
 324           if len(l_cmp_vari) > 0:
 325            for cmp in nom_cmps:
 326             vu = 0
 327             for cmp_vari in l_cmp_vari:
 328              if cmp_vari not in masque_mc:
 329               if cmp == cmp_vari and not vu:
 330                 if EPSI_MOYENNE[mcfacta]['MAILLE'] != None:
 331                   entite = str('MAILLE : '+str(EPSI_MOYENNE[mcfacta]['MAILLE']))
 332                 if EPSI_MOYENNE[mcfacta]['GROUP_MA'] != None:
 333                   entite = str('GROUP_MA : '+str(EPSI_MOYENNE[mcfacta]['GROUP_MA']))
 334                 UTMESS('A','OBSERVATION_8',vali=[num_ordr[ind]],valr=[seuil_mc],valk=[entite,cmp])
 335                 vu = 1
 336
 337
 338           if TYPE_CHAM[-1:] == 'C':
 339             mcfactc = { 'NOM_CMP'   : nom_cmps,
 340                          'NOEUD' : l_noeud_mcfact,
 341                         'VALE_C'  : l_vmoye,
 342                        }
 343           else:
 344             mcfactc = { 'NOM_CMP'   : nom_cmps,
 345                          'NOEUD' : l_noeud_mcfact,
 346                         'VALE'  : l_vmoye,
 347                        }
 348
 349           argsa['AFFE'].append(mcfactc)
 350
 351          __chame[ind] = CREA_CHAMP( OPERATION  = 'AFFE',
 352                           MODELE = MODELE_1,
 353                           PROL_ZERO = 'OUI',
 354                           TYPE_CHAM  = TYPE_CHAM,
 355                           OPTION   = 'EPSI_NOEU',
 356                           **argsa
 357                           );
 358
 359          if isinstance( RESULTAT, mode_meca):
 360                   mcfact2 = {'CHAM_GD'    : __chame[ind],
 361                        'MODELE'     : MODELE_1,
 362                        'NUME_MODE'    : int(afreq[ind,0]),
 363                        'FREQ'    : afreq[ind,1],
 364                        }
 365
 366          if isinstance( RESULTAT, evol_elas):
 367                     mcfact2 = {'CHAM_GD'    : __chame[ind],
 368                        'MODELE'     : MODELE_1,
 369                        'INST'    : list_inst[ind],
 370                        }
 371
 372          if isinstance( RESULTAT, dyna_trans):
 373                     mcfact2 = {'CHAM_GD'    : __chame[ind],
 374                        'MODELE'     : MODELE_1,
 375                        'INST'    : list_inst[ind],
 376                        }
 377
 378          if isinstance( RESULTAT, dyna_harmo):
 379                 mcfact2 = {'CHAM_GD'    : __chame[ind],
 380                            'MODELE'     : MODELE_1,
 381                            'FREQ'    : list_freq[ind],
 382                           }
 383
 384          if cham_mater is not None:
 385                     mcfact2['CHAM_MATER'] = cham_mater
 386          if cara_elem is not None:
 387                     mcfact2['CARA_ELEM'] = cara_elem
 388
 389          liste.append(mcfact2)
 390
 391         resu_epsi = 'EPSI_NOEU'
 392         RESU = CREA_RESU(
 393                               OPERATION = 'AFFE',
 394                               TYPE_RESU = TYPE_RESU,
 395                               NOM_CHAM  = 'EPSI_NOEU',
 396                               AFFE      = liste,
 397                              );
 398
 399
 400 #***********************************************
 401 #  BOUCLE SUR LES NOM_CHAM
 402 #***********************************************
 403
 404     for nomcham in NOM_CHAM:
 405
 406      if nomcham == 'DEPL' or nomcham == 'VITE' or nomcham == 'ACCE':
 407              if isinstance( RESULTAT, dyna_harmo):
 408                   TYPE_CHAM = 'NOEU_DEPL_C'
 409              else:
 410                   TYPE_CHAM = 'NOEU_DEPL_R'
 411      elif nomcham == 'EPSI_NOEU':
 412              if isinstance( RESULTAT, dyna_harmo):
 413                TYPE_CHAM  = 'NOEU_EPSI_C'
 414              else:
 415                TYPE_CHAM  = 'NOEU_EPSI_R'
 416      else:
 417                UTMESS('F', 'ELEMENTS4_48', valk = [nomcham])
 418
 419 #***********************************************
 420 #  PHASE DE PROJECTION
 421 #***********************************************
 422
 423      if PROJECTION == 'OUI':
 424       if resu_epsi and nomcham == 'EPSI_NOEU':
 425         __proj = PROJ_CHAMP(METHODE='COLLOCATION',
 426                           RESULTAT = RESU,
 427                           MODELE_1 = MODELE_1,
 428                           MODELE_2 = MODELE_2,
 429                           NOM_CHAM = nomcham,
 430                           **mcfact
 431                          )
 432       else:
 433         __proj = PROJ_CHAMP(METHODE = 'COLLOCATION',
 434                           RESULTAT = RESULTAT,
 435                           MODELE_1 = MODELE_1,
 436                           MODELE_2 = MODELE_2,
 437                           NOM_CHAM = nomcham,
 438                           **mcfact
 439                          )
 440       modele = MODELE_2
 441      else:
 442       if resu_epsi and nomcham == 'EPSI_NOEU':
 443         __proj = RESU
 444       else:
 445         __proj = RESULTAT
 446       modele = MODELE_1
 447
 448
 449 #***********************************************
 450 #  PHASE DE CHANGEMENT DE REPERE
 451 #***********************************************
 452 # Le changement de repere se fait dans les routines exterieures crea_normale et crea_repere
 453
 454 # On range dans le mcfact MODI_REPERE une liste de modifications. ex :
 455 #     MODI_REPERE = ( _F( GROUP_NO  = toto,
 456 #                         REPERE    = 'NORMALE' ),
 457 #                         CONDITION = (1.0,0.0,0.0),
 458 #                         NOM_PARA  = 'X')
 459 #                     _F( NOEUD   = ('a','z','e'...),
 460 #                         REPERE  = 'CYLINDRIQUE',
 461 #                         ORIGINE = (0.,0.,0.),
 462 #                         AXE_Z   = (0.,1.,0.), ),
 463 #                     _F( GROUP_NO = titi,
 464 #                         REPERE   = 'UTILISATEUR',
 465 #                         ANGL_NAUT = (alpha, beta, gamma), )
 466 #                    )
 467
 468
 469      if MODIF_REPERE != None :
 470         num_ordr = __proj.LIST_VARI_ACCES()['NUME_ORDRE']
 471
 472         for modif_rep in MODIF_REPERE :
 473           modi_rep = modif_rep.val
 474           type_cham = modif_rep['TYPE_CHAM']
 475           nom_cmp = modif_rep['NOM_CMP']
 476
 477           if type_cham == 'TENS_2D':
 478                 nomchamx = 'EPSI_NOEU'
 479           elif type_cham == 'TENS_3D':
 480                 nomchamx = 'EPSI_NOEU'
 481           else:
 482                 nomchamx = 'DEPL'
 483
 484           if nomcham == nomchamx:
 485             mcfact1 = { 'NOM_CMP'   : nom_cmp,
 486                         'TYPE_CHAM' : type_cham,
 487                         'NOM_CHAM'  : nomcham }
 488
 489             mcfact2 = { }
 490
 491             if modi_rep['REPERE'] == 'DIR_JAUGE' :
 492                 vect_x = None
 493                 vect_y = None
 494                 if modi_rep.has_key('VECT_X'):
 495                     vect_x = modi_rep['VECT_X']
 496                 if modi_rep.has_key('VECT_Y'):
 497                     vect_y = modi_rep['VECT_Y']
 498
 499                 # il faut des mailles pour les tenseurs d'ordre 2
 500                 for typ in ['MAILLE','GROUP_MA',]:
 501                   if modi_rep.has_key(typ) :
 502                     if PROJECTION == 'OUI':
 503                       maya = mayaexp
 504                     else:
 505                       maya = mayanum
 506                     list_ma = find_ma(maya, {typ : modi_rep[typ]})
 507
 508                 mcfact1.update({ 'MAILLE'    : list_ma })
 509                 angl_naut = crea_repere_xy(vect_x, vect_y)
 510
 511                 mcfact2.update({ 'REPERE'    : 'UTILISATEUR',
 512                                  'ANGL_NAUT' : angl_naut })
 513
 514                 argsm = {'MODI_CHAM'   : mcfact1,
 515                         'DEFI_REPERE' : mcfact2 }
 516
 517                 __bidon = MODI_REPERE( RESULTAT    = __proj,
 518                                       NUME_ORDRE  = num_ordr,
 519                                       **argsm)
 520                 DETRUIRE( CONCEPT= _F( NOM = __proj ), INFO=1)
 521                 __proj = __bidon
 522
 523             if modi_rep['REPERE'] == 'NORMALE' :
 524                 # Cas ou l'utilisateur choisit de creer les reperes locaux
 525                 # selon la normale. On fait un changement de repere local
 526                 # par noeud
 527                 for option in ['VECT_X','VECT_Y','CONDITION_X','CONDITION_Y'] :
 528                     if modi_rep.has_key(option):
 529                         vect = { option : modi_rep[option] }
 530                 if len(vect) != 1 :
 531                     UTMESS('E','UTILITAI7_9')
 532
 533                 chnorm = crea_normale(self, MODELE_1, MODELE_2, NUME_DDL,
 534                                       cham_mater, cara_elem)
 535                 modele = MODELE_2
 536                 chnormx = chnorm.EXTR_COMP('DX',[],1)
 537                 ind_noeuds = chnormx.noeud
 538                 nom_allno = [mayaexp.NOMNOE.get()[k-1] for k in ind_noeuds]
 539
 540                 # on met les noeuds conernes sous forme de liste et on va
 541                 # chercher les noeuds des mailles pour 'MAILLE' et 'GROUP_MA'
 542                 for typ in ['NOEUD','GROUP_NO','MAILLE','GROUP_MA']:
 543                     if modi_rep.has_key(typ) :
 544                         list_no_exp = find_no(mayaexp, {typ : modi_rep[typ]})
 545
 546                 # boucle sur les noeuds pour modifier les reperes.
 547                 __bid = [None]*(len(list_no_exp) + 1)
 548                 __bid[0] = __proj
 549                 k = 0
 550                 for nomnoe in list_no_exp:
 551                     ind_no = nom_allno.index(nomnoe)
 552                     angl_naut = crea_repere(chnorm, ind_no, vect)
 553
 554                     mcfact1.update({ 'NOEUD'     : nomnoe })
 555                     mcfact2.update({ 'REPERE'    : 'UTILISATEUR',
 556                                      'ANGL_NAUT' : angl_naut})
 557                     argsm = {'MODI_CHAM'   : mcfact1,
 558                             'DEFI_REPERE' : mcfact2 }
 559                     __bid[k+1] = MODI_REPERE( RESULTAT    = __bid[k],
 560                                               TOUT_ORDRE  = 'OUI',
 561                                               CRITERE     = 'RELATIF',
 562                                               **argsm)
 563                     k = k + 1
 564
 565                 __proj = __bid[-1:][0]
 566
 567
 568             if modi_rep['REPERE'] == 'UTILISATEUR' or modi_rep['REPERE'] == 'CYLINDRIQUE':
 569
 570               if type_cham == 'TENS_2D' or type_cham == 'TENS_3D' :
 571                   for typ in ['MAILLE','GROUP_MA',]:
 572                     if modi_rep.has_key(typ) :
 573                         mcfact1.update({typ : modi_rep[typ]})
 574               else:
 575                 for typ in ['NOEUD','GROUP_NO','MAILLE','GROUP_MA']:
 576                     if modi_rep.has_key(typ) :
 577                         mcfact1.update({typ : modi_rep[typ]})
 578
 579               if modi_rep['REPERE'] == 'CYLINDRIQUE' :
 580                     origine = modi_rep['ORIGINE']
 581                     axe_z   = modi_rep['AXE_Z']
 582                     mcfact2.update({ 'REPERE'  : 'CYLINDRIQUE',
 583                                      'ORIGINE' : origine,
 584                                      'AXE_Z'   : axe_z })
 585
 586               elif modi_rep['REPERE'] == 'UTILISATEUR' :
 587                     angl_naut = modi_rep['ANGL_NAUT']
 588                     mcfact2.update({ 'REPERE'    : 'UTILISATEUR',
 589                                      'ANGL_NAUT' : angl_naut })
 590
 591               argsm = {'MODI_CHAM'   : mcfact1,
 592                       'DEFI_REPERE' : mcfact2 }
 593
 594               __bidon = MODI_REPERE( RESULTAT    = __proj,
 595                                      CRITERE     = 'RELATIF',
 596                                      **argsm)
 597               DETRUIRE( CONCEPT= _F( NOM = __proj ), INFO=1)
 598               __proj = __bidon
 599
 600
 601 #*************************************************
 602 # Phase de selection des DDL de mesure
 603 #*************************************************
 604
 605      resu_filtre = None
 606      if FILTRE != None:
 607          num_ordr = __proj.LIST_VARI_ACCES()['NUME_ORDRE']
 608
 609          __chamf = [None]*len(indice)
 610          if isinstance( RESULTAT, evol_elas):
 611             list_inst = __proj.LIST_VARI_ACCES()['INST']
 612          if isinstance( RESULTAT, dyna_trans):
 613             list_inst = __proj.LIST_VARI_ACCES()['INST']
 614          if isinstance( RESULTAT, dyna_harmo):
 615             list_freq = __proj.LIST_VARI_ACCES()['FREQ']
 616
 617
 618          liste = []
 619
 620          for ind in indice:
 621              mcfact2 = { }
 622              filtres = []
 623              __chamex = CREA_CHAMP(TYPE_CHAM  = TYPE_CHAM,
 624                                   OPERATION  = 'EXTR',
 625                                   RESULTAT   = __proj,
 626                                   NOM_CHAM   = nomcham,
 627                                   NUME_ORDRE = num_ordr[ind],);
 628
 629              for mcfiltre in FILTRE :
 630               filtre = mcfiltre.val
 631               try:
 632                 nomchamx = filtre['NOM_CHAM']
 633               except KeyError:
 634                 nomchamx = None
 635
 636               if nomchamx == None or nomchamx == nomcham:
 637
 638                 mcfact1 = {}
 639
 640                 atraiter = None
 641                 if filtre['DDL_ACTIF'][0][0] == 'E' and nomcham == 'EPSI_NOEU':
 642                    atraiter = nomcham
 643                 elif filtre['DDL_ACTIF'][0][0] == 'D' and nomcham == 'DEPL':
 644                    atraiter = nomcham
 645                 elif filtre['DDL_ACTIF'][0][0] == 'D' and nomcham == 'VITE':
 646                    atraiter = nomcham
 647                 elif filtre['DDL_ACTIF'][0][0] == 'D' and nomcham == 'ACCE':
 648                    atraiter = nomcham
 649
 650                 if atraiter:
 651                  for typ in ['NOEUD', 'GROUP_NO', 'MAILLE', 'GROUP_MA']:
 652                     if filtre.has_key(typ) :
 653                         mcfact1.update({typ : filtre[typ]})
 654                  mcfact1.update({'NOM_CMP' : filtre['DDL_ACTIF'],
 655                                 'CHAM_GD' : __chamex })
 656                  filtres.append(mcfact1)
 657
 658              if len(filtres) > 0:
 659               if nomcham == 'DEPL' or nomcham == 'VITE' or nomcham == 'ACCE':
 660                 __chamf[ind] = CREA_CHAMP(TYPE_CHAM = TYPE_CHAM,
 661                                        OPERATION = 'ASSE',
 662                                        MODELE    = modele,
 663                                        ASSE      = filtres
 664                                        );
 665
 666               elif nomcham == 'EPSI_NOEU':
 667                 __chamf[ind] = CREA_CHAMP(TYPE_CHAM = TYPE_CHAM,
 668                                        OPERATION = 'ASSE',
 669                                        PROL_ZERO = 'OUI',
 670                                        MODELE    = modele,
 671                                        ASSE      = filtres,
 672                                        );
 673               else:
 674                 valk = []
 675                 valk.append(nomcham)
 676                 valk.append('DEPL VITE ACCE')
 677                 valk.append('EPSI_NOEU')
 678                 UTMESS('F','OBSERVATION_6',valk)
 679
 680               argsr = {}
 681               if isinstance( RESULTAT, mode_meca):
 682                 mcfact2 = {'CHAM_GD' : __chamf[ind],
 683                        'MODELE'      : modele,
 684                        'NUME_MODE'   : int(afreq[ind,0]),
 685                        'FREQ'        : afreq[ind,1],
 686                        }
 687             ## on recopie les matrices associees au MODELE_2 dans le resultat final
 688             ## NB : ce n'est peut-etre pas propre, car ces matrices ne
 689             ## veulent plus rien dire si on a filtre des DDL !!!!!
 690                 argsr = { 'MATR_A' : MATR_A,
 691                          'MATR_B' : MATR_B
 692                         }
 693
 694               if isinstance( RESULTAT, evol_elas):
 695                 mcfact2 = {'CHAM_GD'    : __chamf[ind],
 696                        'MODELE'     : modele,
 697                        'INST'    : list_inst[ind],
 698                        }
 699
 700               if isinstance( RESULTAT, dyna_trans):
 701                 mcfact2 = {'CHAM_GD'    : __chamf[ind],
 702                        'MODELE'     : modele,
 703                        'INST'    : list_inst[ind],
 704                        }
 705
 706               if isinstance( RESULTAT, dyna_harmo):
 707                 mcfact2 = {'CHAM_GD'    : __chamf[ind],
 708                        'MODELE'     : MODELE_2,
 709                        'FREQ'    : list_freq[ind],
 710                        }
 711
 712               if cham_mater is not None:
 713                 mcfact2['CHAM_MATER'] = cham_mater
 714               if cara_elem is not None:
 715                 mcfact2['CARA_ELEM'] = cara_elem
 716
 717               liste.append(mcfact2)
 718
 719              DETRUIRE( CONCEPT= _F( NOM = __chamex ), INFO=1)
 720
 721          if len(filtres) > 0 and len(liste) > 0:
 722             resu_filtre = nomcham
 723             if RESU:
 724              RESU = CREA_RESU( reuse = RESU,
 725                               OPERATION = 'AFFE',
 726                               TYPE_RESU = TYPE_RESU,
 727                               NOM_CHAM  = nomcham,
 728                               AFFE      = liste,
 729                               **argsr
 730                              );
 731             else:
 732              RESU = CREA_RESU(
 733                               OPERATION = 'AFFE',
 734                               TYPE_RESU = TYPE_RESU,
 735                               NOM_CHAM  = nomcham,
 736                               AFFE      = liste,
 737                               **argsr
 738                              )
 739
 740
 741 #*************************************************
 742 # Recopie de __proj dans RESU si celle-ci
 743 # n'a pas encore ete faite via FILTRE
 744 #*************************************************
 745
 746      if resu_filtre == None:
 747         RESU = PROJ_CHAMP(METHODE = 'COLLOCATION',
 748                           RESULTAT = __proj,
 749                           MODELE_1 = modele,
 750                           MODELE_2 = modele,
 751                           NUME_DDL = NUME_DDL,
 752                           NOM_CHAM = nomcham,
 753                           **mcfact
 754                          )
 755
 756
 757     return ier
 758
 759
 760
 761
 762 #**********************************************
 763 # RECUPERATION DES NORMALES
 764 #**********************************************
 765
 766 def crea_normale(self, modele_1, modele_2,
 767                  nume_ddl, cham_mater=None, cara_elem=None):
 768     """Cree un champ de vecteurs normaux sur le maillage experimental, par
 769        projection du champ de normales cree sur le maillage numerique
 770        les mailles doivent etre des elements de <peau> (facettes)
 771     """
 772     PROJ_CHAMP  = self.get_cmd('PROJ_CHAMP')
 773     CREA_CHAMP  = self.get_cmd('CREA_CHAMP')
 774     CREA_RESU   = self.get_cmd('CREA_RESU')
 775     DEFI_GROUP  = self.get_cmd('DEFI_GROUP')
 776     import aster
 777     from Accas import _F
 778     # recherche du maillage associe au modele numerique
 779     nom_modele_num = modele_1.nom
 780     _maillag = aster.getvectjev( nom_modele_num.ljust(8) + '.MODELE    .LGRF' )
 781     maillage = _maillag[0].strip()
 782     mayanum = self.get_concept(maillage)
 783
 784
 785     DEFI_GROUP( reuse = mayanum,
 786                 MAILLAGE      = mayanum,
 787                 CREA_GROUP_MA = _F( NOM  = '&&TOUMAI',
 788                                     TOUT = 'OUI' )
 789                );
 790
 791     __norm1 = CREA_CHAMP( MODELE    = modele_1,
 792                           OPERATION = 'NORMALE',
 793                           TYPE_CHAM = 'NOEU_GEOM_R',
 794                           GROUP_MA  = '&&TOUMAI',
 795                          );
 796
 797     DEFI_GROUP( reuse = mayanum,
 798                 MAILLAGE      = mayanum,
 799                 DETR_GROUP_MA = _F( NOM  = '&&TOUMAI' )
 800                );
 801
 802
 803     __norm2 = CREA_CHAMP( OPERATION = 'ASSE',
 804                           TYPE_CHAM = 'NOEU_DEPL_R',
 805                           MODELE    = modele_1,
 806                           ASSE      = _F( TOUT='OUI',
 807                                           CHAM_GD=__norm1,
 808                                           NOM_CMP=('X','Y','Z'),
 809                                           NOM_CMP_RESU=('DX','DY','DZ')
 810                                          )
 811                          );
 812
 813     affe_dct = {'CHAM_GD' : __norm2,
 814                 'INST' : 1,
 815                 'MODELE' : modele_1}
 816     if cham_mater is not None:
 817         affe_dct["CHAM_MATER"] = cham_mater
 818     if cara_elem is not None:
 819         affe_dct["CARA_ELEM"] = cara_elem
 820
 821     __norm3 = CREA_RESU( OPERATION = 'AFFE',
 822                          TYPE_RESU = 'EVOL_ELAS',
 823                          NOM_CHAM  = 'DEPL',
 824                          AFFE      = _F(**affe_dct)
 825                        );
 826
 827
 828     __norm4 = PROJ_CHAMP( METHODE='COLLOCATION',
 829                           RESULTAT   = __norm3,
 830                           MODELE_1   = modele_1,
 831                           MODELE_2   = modele_2,
 832                           NOM_CHAM   = 'DEPL',
 833                           TOUT_ORDRE = 'OUI',
 834                           NUME_DDL   = nume_ddl,
 835                           );
 836
 837     # __norm5 : toutes les normales au maillage au niveau des capteurs
 838     __norm5 = CREA_CHAMP( RESULTAT   = __norm4,
 839                           OPERATION  = 'EXTR',
 840                           NUME_ORDRE = 1,
 841                           NOM_CHAM   = 'DEPL',
 842                           TYPE_CHAM  = 'NOEU_DEPL_R',
 843                           );
 844
 845
 846     return __norm5
 847
 848 #**********************************************************************
 849 # Calcul des angles nautiques pour le repere local associe a la normale
 850 #**********************************************************************
 851
 852 def crea_repere(chnorm, ind_no, vect):
 853
 854     """Creation d'un repere orthonormal a partir du vecteur normale et
 855        d'une equation supplementaire donnee par l'utilisateur sous forme
 856        de trois parametres et du vecteur de base concerne.
 857     """
 858
 859     import numpy
 860
 861     nom_para = vect.keys()[0] # nom_para = 'VECT_X' ou 'VECT_Y'
 862
 863     # 1) pour tous les noeuds du maillage experimental, recuperer la normale
 864     #    calculee a partir du maillage numerique
 865     chnormx = chnorm.EXTR_COMP('DX',[],1)
 866     chnormy = chnorm.EXTR_COMP('DY',[],1)
 867     chnormz = chnorm.EXTR_COMP('DZ',[],1)
 868
 869     normale = [chnormx.valeurs[ind_no],
 870                chnormy.valeurs[ind_no],
 871                chnormz.valeurs[ind_no]]
 872
 873     # 2.1) soit l'utilisateur a donne un deuxieme vecteur explicitement
 874     # (option VECT_X Ou VECT_Y). Dans ce cas la, le 3e est le produit
 875     # vectoriel des deux premiers.
 876     if nom_para == 'VECT_X' :
 877         vect1 = numpy.array(list(vect[nom_para])) # vect x du reploc
 878         vect2 = numpy.cross(normale,vect1)
 879         reploc = numpy.array([vect1.tolist(), vect2.tolist(), normale])
 880         reploc = numpy.transpose(reploc)
 881
 882     elif nom_para == 'VECT_Y' :
 883         vect2 = numpy.array(list(vect[nom_para])) # vect y du reploc
 884         vect1 = numpy.cross(vect2, normale)
 885         reploc = numpy.array([vect1.tolist(), vect2.tolist(), normale])
 886         reploc = numpy.transpose(reploc)
 887
 888     # 2.2) TODO : plutot que de donner explicitement un vecteur du repere
 889     # local avec VECT_X/Y, on devrait aussi pouvoir donner une condition
 890     # sous forme d'une equation sur un de ces vecteurs. Par exemple,
 891     # CONDITION_X = (0.,1.,0.) signifierait que le vecteur X1 verifie
 892     #x(X1) + y(X1) + z(X1) = 0
 893     elif nom_para == 'CONDITION_X':
 894         pass
 895     elif nom_para == 'CONDITION_Y':
 896         pass
 897
 898     # 3) Calcul de l'angle nautique associe au repere local
 899     angl_naut = anglnaut(reploc)
 900
 901     return angl_naut
 902
 903 #**********************************************************************
 904 # Calcul des angles nautiques pour le repere associe a VECT_X et VECT_Y
 905 #**********************************************************************
 906
 907 def crea_repere_xy(vect_x, vect_y):
 908
 909     """Calcul des angles nautiques a partir des directions vect_x et vect_y.
 910        Si vect_x != None et vect_y != None alors on impose le premier vecteur de base
 911            colineaire a vect_x et le deuxieme vecteur dans le plan (vect_x,vect_y)
 912        Si vect_x != None et vect_y == None alors on impose le premier vecteur de base
 913            colineaire a vect_x
 914        Si vect_x == None et vect_y != None alors on impose le deuxieme vecteur de base
 915            colineaire a vect_y
 916        Si vect_x == None et vect_y == None alors on ne fait rien
 917     """
 918
 919     import numpy
 920     from Utilitai.Utmess import UTMESS
 921
 922     if vect_x == None and vect_y == None:
 923       angl_naut = (0.,0.,0.)
 924     else:
 925       if vect_x and vect_y:
 926         vx = numpy.array(list(vect_x))
 927         vy = numpy.array(list(vect_y))
 928         vect1 = vx
 929         vect3 = numpy.cross(vx,vy)
 930         vect2 = numpy.cross(vect3,vx)
 931
 932       elif vect_x:
 933         vx = numpy.array(list(vect_x))
 934         vy1 = numpy.cross((1.,0.,0.),vx)
 935         vy2 = numpy.cross((0.,1.,0.),vx)
 936         vy3 = numpy.cross((0.,0.,1.),vx)
 937         n1 = norm(vy1)
 938         n2 = norm(vy2)
 939         n3 = norm(vy3)
 940         nmax = max(n1,n2,n3)
 941         if nmax == n1:
 942             vy = vy1
 943         elif nmax == n2:
 944             vy = vy2
 945         elif nmax == n3:
 946             vy = vy3
 947         else:
 948             UTMESS('F','UTILITAI_7')
 949         vect3 = numpy.cross(vx,vy)
 950         vect1 = vx
 951         vect2 = numpy.cross(vect3,vect1)
 952
 953       elif vect_y:
 954         vy = numpy.array(list(vect_y))
 955         vx1 = numpy.cross((1.,0.,0.),vy)
 956         vx2 = numpy.cross((0.,1.,0.),vy)
 957         vx3 = numpy.cross((0.,0.,1.),vy)
 958         n1 = norm(vx1)
 959         n2 = norm(vx2)
 960         n3 = norm(vx3)
 961         nmax = max(n1,n2,n3)
 962         if nmax == n1:
 963             vx = vx1
 964         elif nmax == n2:
 965             vx = vx2
 966         elif nmax == n3:
 967             vx = vx3
 968         else:
 969             UTMESS('F','UTILITAI_7')
 970         vect3 = numpy.cross(vx,vy)
 971         vect2 = vy
 972         vect1 = numpy.cross(vect2, vect3)
 973
 974
 975       norm12=numpy.dot(vect1,vect1)
 976       norm22=numpy.dot(vect2,vect2)
 977       norm32=numpy.dot(vect3,vect3)
 978       if norm12 == 0 or norm22 == 0 or norm32 == 0:
 979           UTMESS('F','UTILITAI_7')
 980       else:
 981           reploc = numpy.array([vect1.tolist(),vect2.tolist(),vect3.tolist()])
 982           reploc = numpy.transpose(reploc)
 983           angl_naut = anglnaut(reploc)
 984
 985     return angl_naut
 986
 987 #*****************************************************************************
 988 # Aller chercher une liste de noeuds pour un mot cle 'NOEUD', 'GROUP_NO'
 989 # 'MAILLE' ou 'GROUP_MA'
 990 #*****************************************************************************
 991
 992 def find_no(maya,mcsimp):
 993     """ Si on demande une liste de noeuds, c'est simple, on retourne les noeuds
 994         Si on demande une liste de groupes de noeuds, on va chercher les noeuds
 995           dans ces groupes, en faisant attention a ne pas etre redondant
 996         Si on demande un liste de mailles, on va chercher dans le .CONNEX
 997           du maillage les indices, puis les noms des noeuds concernes
 998         etc...
 999     """
1000
1001     import numpy
1002
1003     list_no = []
1004     if mcsimp.has_key('GROUP_NO') and type(mcsimp['GROUP_NO']) != tuple :
1005         mcsimp['GROUP_NO'] = [mcsimp['GROUP_NO']]
1006     if mcsimp.has_key('MAILLE') and type(mcsimp['MAILLE']) != tuple :
1007         mcsimp['MAILLE'] = [mcsimp['MAILLE']]
1008     if mcsimp.has_key('GROUP_MA') and type(mcsimp['GROUP_MA']) != tuple :
1009         mcsimp['GROUP_MA'] = [mcsimp['GROUP_MA']]
1010
1011     if mcsimp.has_key('NOEUD') :
1012         list_no = list(mcsimp['NOEUD'])
1013     elif mcsimp.has_key('GROUP_NO') :
1014         for group in mcsimp['GROUP_NO'] :
1015             list_ind_no = list(numpy.array(maya.GROUPENO.get()[group.ljust(8)])-1)
1016             for ind_no in list_ind_no :
1017                 nomnoe = maya.NOMNOE.get()[ind_no]
1018                 if nomnoe not in list_no :
1019                     list_no.append(nomnoe)
1020     elif mcsimp.has_key('MAILLE') :
1021         for mail in mcsimp['MAILLE'] :
1022             for index in range(len(maya.NOMMAI.get())):
1023                 if maya.NOMMAI.get()[index].strip() == mail:
1024                     nu_ma = index
1025             for ind_no in maya.CONNEX.get()[nu_ma+1]:
1026                     nomnoe = maya.NOMNOE.get()[ind_no-1]
1027                     if nomnoe not in list_no:
1028                         list_no.append(nomnoe)
1029     elif mcsimp.has_key('GROUP_MA') :
1030         for group in mcsimp['GROUP_MA'] :
1031             list_nu_ma = list(numpy.array(maya.GROUPEMA.get()
1032                                             [group.ljust(8)]) - 1)
1033             for nu_ma in list_nu_ma:
1034                 for ind_no in maya.CONNEX.get()[nu_ma+1]:
1035                     nomnoe = maya.NOMNOE.get()[ind_no-1]
1036                     if nomnoe not in list_no:
1037                         list_no.append(nomnoe)
1038
1039     return list_no
1040
1041 #*****************************************************************************
1042 # Aller chercher une liste de mailles pour un mot cle 'MAILLE' ou 'GROUP_MA'
1043 #*****************************************************************************
1044
1045 def find_ma(maya,mcsimp):
1046     """ Si mot cle MAILLE, on retourne la liste des mailles
1047         Si mot cle GROUP_MA, on va chercher les mailles dans ces groupes
1048     """
1049
1050     import numpy
1051
1052     list_ma = []
1053     if mcsimp.has_key('GROUP_MA') and type(mcsimp['GROUP_MA']) != tuple :
1054         mcsimp['GROUP_MA'] = [mcsimp['GROUP_MA']]
1055     if mcsimp.has_key('MAILLE') and type(mcsimp['MAILLE']) != tuple :
1056         mcsimp['MAILLE'] = [mcsimp['MAILLE']]
1057
1058     if mcsimp.has_key('MAILLE') :
1059         for mail in mcsimp['MAILLE'] :
1060             list_ma.append(mail)
1061     elif mcsimp.has_key('GROUP_MA') :
1062         for group in mcsimp['GROUP_MA'] :
1063             list_ind_ma = list(numpy.array(maya.GROUPEMA.get()[group.ljust(8)])-1)
1064             for ind_ma in list_ind_ma :
1065                 nommail = maya.NOMMAI.get()[ind_ma]
1066                 if nommail not in list_ma :
1067                     list_ma.append(nommail)
1068
1069     return list_ma
1070
1071
1072 #************************************************************************************
1073 # Quelques utilitaires de calculs d'angles nautiques (viennent de zmat004a.comm
1074 #************************************************************************************
1075
1076 def norm(x):
1077     """Calcul de la norme euclidienne d'un vecteur"""
1078     import numpy
1079     tmp = numpy.sqrt(numpy.dot(x,x))
1080     return tmp
1081
1082 def anglnaut(P):
1083
1084
1085     """Calcule les angles nautiques correspondant a un repere local
1086        NB : seuls les deux premiers vecteurs de P (les images respectives
1087        de X et Y) sont utiles pour le calcul des angles
1088     """
1089
1090     import copy
1091     import numpy
1092     # expression des coordonnees globales des 3 vecteurs de base locale
1093     y = numpy.array([0.,1.,0.])
1094
1095     xg = P[:,0]
1096     yg = P[:,1]
1097     zg = P[:,2]
1098
1099     # calcul des angles nautiques
1100     x1=copy.copy(xg)
1101     # x1 projection de xg sur la plan xy, non norme
1102     x1[2]=0.
1103     # produit scalaire X xg
1104     normx = norm(x1)
1105     if normx == 0.: # on impose alpha = 0 pour lever l'indetermination
1106         COSA=1.
1107         SINA=0.
1108     else:
1109         COSA=x1[0]/normx
1110         #produit vectoriel X xg
1111         SINA=x1[1]/normx
1112     ar=numpy.arctan2(SINA,COSA)
1113     alpha=ar*180/numpy.pi
1114
1115     COSB=norm(x1)
1116     SINB=-xg[2]
1117     beta=numpy.arctan2(SINB,COSB)*180/numpy.pi
1118
1119     P2=numpy.zeros((3,3))
1120     P2[0,0]=numpy.cos(ar)
1121     P2[1,0]=numpy.sin(ar)
1122     P2[1,1]=numpy.cos(ar)
1123     P2[0,1]=-numpy.sin(ar)
1124     y1=numpy.dot(P2,y)
1125     y1n=y1/norm(y1)
1126
1127     # calcul de gamma
1128     COSG=numpy.dot(y1n,yg)
1129     SING=numpy.dot(xg,numpy.cross(y1n,yg))
1130     gamma=numpy.arctan2(SING,COSG)*180/numpy.pi
1131
1132     return alpha,beta,gamma
1133
1134
1135 ##  NB : Equations de passage : un vecteur de coordonnees globales (X,Y,Z) a pour
1136 ##  coordonnees locales (X1,Y1,Z1) avec
1137 ##  _                  _  _                   _  _                   _  _ _     _  _
1138 ## | 1     0      0     || cos(B) 0    -sin(B) ||  cos(A)  sin(A)   0 || X |   | X1 |
1139 ## | 0   cos(G)  sin(G) ||   0    1      0     || -sin(A)  cos(A)   0 || Y | = | Y1 |
1140 ## |_0  -sin(G)  cos(G)_||_sin(B) 0     cos(B)_||_   0       0      1_||_Z_|   |_Z1_|
1141 ##
1142 ##  A (alpha), B(beta), gamma (G) sont les angle nautiques que l'on donne habituellemet
1143 ##  dans les MODI_REPERE. Les equations a resoudre sont les suivantes :
1144 ##       cos(A)cos(B)                      = reploc[0][0]
1145 ##      -cos(G)sin(A) + sin(G)cos(A)sin(B) = reploc[0][1]
1146 ##       sin(A)sin(G) + cos(A)sin(B)cos(G) = reploc[0][2]
1147 ##
1148 ##       sin(A)cos(B)                      = reploc[1][0]
1149 ##       cos(A)cos(G) + sin(A)sin(B)sin(G) = reploc[1][1]
1150 ##      -cos(A)sin(G) + sin(A)sin(B)cos(G) = reploc[1][2]
1151 ##
1152 ##                                 -sin(B) = reploc[2][0]
1153 ##                            cos(B)sin(G) = reploc[2][1]
1154 ##                            cos(B)cos(G) = reploc[2][2]
1155
1156
1157