Aster/Cata/cataSTA10/Macro/observation_ops.py

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