Télécharger la plaquette Cast3M


fr en

Télécharger proch3.eso

Retour à la liste

Numérotation des lignes : 

proch3

  1. C PROCH3    SOURCE    CB215821  20/11/25    13:37:13     10792          
  2. C     SUBROUTINE PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE,
  3. C     INSYM, MTAB3, I )
  4. ************************************************************************
  5. *
  6. *                             P R O C H 3
  7. *                             -----------
  8. *
  9. * FONCTION:
  10. * ---------
  11. *
  12. *     RECHERCHE DES NBMOD MODES PROPRES LES PLUS PROCHES DE FREQ
  13. *     ET CONSTRUIT L'OBJET 'SOLUTION' CORRESPONDANT.
  14. *
  15. * MODE D'APPEL:
  16. * -------------
  17. *
  18. *     CALL   PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE, INSYM )
  19. *
  20. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE
  21. * -----------
  22. *
  23. *     FREQ    REEL DP   (E)  FREQUENCE A APPROCHER PAR UNE FREQUENCE
  24. *                            PROPRE.
  25. *     NBMOD   ENTIER    (E)  ORDRE DE MULTIPLICITE DE LA FREQ PROPRE
  26. *     IPRIGI  ENTIER    (E)  POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT
  27. *                            LA MATRICE DE RIGIDITE.
  28. *     IPMASS  ENTIER    (E)  POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT
  29. *                            LA MATRICE MASSE.
  30. *     IPMODE  ENTIER    (S)  POINTEUR DE L'OBJET 'SOLUTION' CONTENANT
  31. *                            LE MODE PROPRE TROUVE.
  32. *     LIMAGE  BOOLEEN   (E)  VRAI SI ON SOUHAITE DE FREQ. NEGATIVES.
  33. *
  34. * LEXIQUE: (ORDRE ALPHABETIQUE)
  35. * --------
  36. *
  37. *     IPKW2M  ENTIER    POINTEUR DE LA 'RIGIDITE' "DECALEE"  K - W2.M
  38. *     IPLVER  LISTCHPO  LISTE DE CHPO A ITERER. ILS FORMENT UN SOUS-
  39. *                       ESPACE QUI CONVERGE VERS LE SOUS-ESPACE PROPRE
  40. *                       DOMINANT DE DIMENSION NBMOD. SORTIE: PARTIE REELLE
  41. *                       DES CHPOINTS SOLUTIONS
  42. *     IPLVEI  LISTCHPO  LISTE DE CHPO.SORTIE: PARTIE IMAGINAIRE
  43. *                       DES CHPOINTS SOLUTIONS
  44. *     IPVECP  ENTIER    POINTEUR DU 'CHPOINT' QUI CONTIENT DES NOMBRES
  45. *                       ALEATOIRES, PUIS UN VECTEUR PROPRE.
  46. *     OMEGA2  REEL DP   PULSATION PROPRE TROUVEE AU CARRE.
  47. *                       ( NE SERT A RIEN ! )
  48. *     PROPRE  REEL DP   TABLEAU CONTENANT DES CARACTERISTIQUES DU
  49. *                       MODE PROPRE CALCULE:
  50. *                            PROPRE(1) = FREQ. PROPRE ,
  51. *                            PROPRE(2) = MASSE GEN.
  52. *                            PROPRE(3)ET(4) ET(5) DEPL.GEN. SELON X,Y,Z
  53. *     W2      REEL DP   PULSATION AU CARRE A APPROCHER.
  54. *
  55. * MODE DE FONCTIONNEMENT:
  56. * -----------------------
  57. *
  58. *     LE CALCUL D'UN VECTEUR PROPRE SE FAIT PAR LA METHODE DE
  59. *     L'ITERATION D'UN SOUS-ESPACE, AVEC DECALAGE INITIAL ("SHIFTING")
  60. *     DE W2
  61. *
  62. * AUTEUR, DATE DE CREATION:
  63. * -------------------------
  64. *
  65. *     PASCAL MANIGOT     16 OCTOBRE 1984
  66. *
  67. * MODIF: Benoit PRABEL   MARS 2009
  68. *
  69. * LANGAGE:
  70. * --------
  71. *
  72. *     ESOPE
  73. *
  74. ************************************************************************
  75.  
  76.       SUBROUTINE PROCH3 ( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE
  77.      & ,INSYM, MTAB3, I)
  78.  
  79.       IMPLICIT INTEGER(I-N)
  80.       IMPLICIT REAL*8 (A-H,O-Z)
  81.  
  82. -INC CCREEL
  83. -INC PPARAM
  84. -INC CCOPTIO
  85. -INC SMLCHPO
  86. -INC SMLREEL
  87. -INC SMCHPOI
  88.  
  89. *     -- CONSTANTES --
  90.  
  91.       PARAMETER (LPROPR = 10)
  92.       PARAMETER (DEUXPI = (2.D0*XPI))
  93.  
  94.  
  95. *     -- ARGUMENTS --
  96.  
  97.       REAL*8  FREQ
  98.       INTEGER NBMOD, IPRIGI, IPMASS, INF0, IPMODE
  99.       LOGICAL LIMAGE
  100.       CHARACTER*(LOCOMP) MOTCLE
  101.  
  102.  
  103. *     -- VARIABLES LOCALES --
  104.  
  105.       POINTEUR IPLVER.MLCHPO, IPLVE1.MLCHPO, IPLVEI.MLCHPO
  106.       POINTEUR IPLVAR.MLREEL, IPLVA1.MLREEL, IPLVAI.MLREEL
  107.       INTEGER  NBVEC
  108.       REAL*8   OMEGA2, PROPRE(LPROPR)
  109.  
  110. ******************************************************************
  111. *                   INITIALISATION                               *
  112. ******************************************************************
  113.  
  114.  
  115. *     -- CREATION DE (K-W2M) --
  116.  
  117.       W2 = (FREQ * DEUXPI) ** 2
  118.       IF( LIMAGE ) THEN
  119.          W2 = SIGN(W2,FREQ)
  120.       ENDIF
  121.       CALL DECALE ( IPRIGI, IPMASS, W2, IPKW2M )
  122.       IF ( IERR .NE. 0 ) RETURN
  123.  
  124.  
  125. *     -- INITIALISATION DES ITERATIONS: CREATION D'UNE LISTE
  126. *        DE  NBVEC 'CHPOINT' ALEATOIRE --
  127.  
  128.  
  129. *     bp (03.2009) : si on souhaite eventuellement fournir 1 CHPOINT
  130.       call LIROBJ('LISTCHPO',MLCHP1,0,IRET1)
  131.       if(IRET1 .ne. 0) then
  132.        segact,MLCHP1
  133.        NCHPO1 = MLCHP1.ICHPOI(/1)
  134.        IPVECP = 0
  135. *       call COPIE2(MLCHP1.ICHPOI(1),IPVECP)
  136.        ICHP1 = MLCHP1.ICHPOI(1)
  137.        call COPIE2(ICHP1,IPVECP)
  138. *      A Faire : vérifier l'accord entre chpoint et matrices...?
  139.       else
  140.        CALL ALEAT1 (IPKW2M, IPVECP )
  141.       endif
  142.        IF ( IERR .NE. 0 ) RETURN
  143.        N1 = 1
  144.        SEGINI ,IPLVER
  145.        IPLVER.ICHPOI(1) = IPVECP
  146. **     on cherce le nombre de ddl sans relations
  147.        nddl = 0
  148.        ndlx = 0
  149.        mchpoi = ipvecp
  150.        segact mchpoi
  151.        nsous = ipchp(/1)
  152.         do 81 isous = 1, nsous
  153.           msoupo = ipchp(isous)
  154.           segact msoupo
  155.           mpoval = ipoval
  156.           segact mpoval
  157.           nn = vpocha(/1)
  158.           nc1 = vpocha(/2)
  159.           nddl = nddl + nn*nc1
  160.         do 95 inc = 1,nc1
  161.         if (nocomp(inc).eq. 'LX  ') ndlx = ndlx + nn
  162.  95     continue
  163.         segdes msoupo, mpoval
  164.  81     continue
  165.         segdes mchpoi
  166.  
  167.       nddl = nddl - ndlx - (ndlx/2)
  168.       nbvec = max(2*nbmod, 8)
  169.       if (nbvec .gt. nddl) nbvec = nddl
  170.  
  171. *     on continue le tirage au sort (nous aurions pu passer par Lanzos!)
  172.       DO 100 IB100 = 1, NBVEC-1
  173.          if ( (IRET1 .ne. 0) .and. (IB100 .lt. NCHPO1) ) then
  174.            IPVECP = 0
  175. *           call COPIE2(MLCHP1.ICHPOI(1 + IB100),IPVECP)
  176.            ICHP1 = MLCHP1.ICHPOI(1 + IB100)
  177.            call COPIE2(ICHP1,IPVECP)
  178.          else
  179.            CALL ALEAT1 (IPKW2M, IPVECP )
  180.          endif
  181.          IF ( IERR .NE. 0 ) RETURN
  182.          IPLVER.ICHPOI(**) = IPVECP
  183.  100  CONTINUE
  184.       SEGDES ,IPLVER
  185.       if(IRET1 .ne. 0)  SEGDES,MLCHP1
  186.  
  187. * cette orthonormalisation n'est pas indispensable, dans la plupart
  188. * des cas,
  189. * mais, on l'a gardée pour tester la singularité de la masse
  190. * et mettre à jour nbvec si necessaire. On n'a besoin de ce test
  191. * qu'au debut, car si l'on a eliminé les eventuels vecteurs du noyau
  192. * de la masse ils ne peuvent pas ressurgir lors du processus iteratif
  193. * car ils correspondent à des frequences infinies
  194. *      IF ( INSYM .EQ. 0) THEN
  195.          CALL GRAAMO ( IPLVER, IPMASS )
  196.          IF ( IERR .NE. 0 ) RETURN
  197.           segact IPLVER
  198.           nbvec = IPLVER.ICHPOI(/1)
  199.           segdes iplver
  200.           if (nbmod .gt. nbvec) nbmod = nbvec
  201. *       END IF
  202.  
  203.  
  204. ******************************************************************
  205. *                        RESOLUTION                              *
  206. ******************************************************************
  207.  
  208. *     -- RESOLUTION PAR ITERATION DU SOUS-ESPACE IPLVER
  209. *     ET NORMALISATION DES MODES PROPRES --
  210.  
  211. *on differentie NBMOD=dim du sous espace initialement recherché et
  212. *NBMOD2 qui peut etre modifié par sespac si dernier vecteur complexe
  213.       NBMOD2 = NBMOD
  214.  
  215.       IF (INSYM .EQ. 0) THEN
  216.  
  217.          CALL SESPA  ( IPLVAR, IPLVER, NBMOD, IPKW2M, IPMASS )
  218.          IF ( IERR .NE. 0 ) RETURN
  219.          segact iplver
  220.          nblgn = IPLVER.ICHPOI(/1)
  221.          ICOMP = 0
  222.          mlchpo = iplver
  223.          segact mlchpo*mod
  224.          do 101 ibmod = 1,nbmod
  225.                ix = ichpoi(ibmod)
  226.                CALL MOTS1( IPLMOT, MOTCLE )
  227.                call norma1 (ix,iplmot,motcle,ixx)
  228.                call dtchpo(ix)
  229.                ichpoi(ibmod) = ixx
  230.  101     continue
  231.          segdes mlchpo
  232.  
  233.         ELSE
  234.  
  235.            CALL SESPAC(IPLVAR,IPLVAI,IPLVER,NBMOD2,IPKW2M,IPMASS)
  236.            ICOMP = 0
  237. *
  238.         END IF
  239.  
  240.  
  241. ******************************************************************
  242. *      CREATION DE L'OBJET SOLUTION                              *
  243. ******************************************************************
  244.  
  245. *     -- TRANSLATION DES FREQUENCES PROPRES -- 
  246.  
  247.       SEGACT ,IPLVAR*MOD
  248. *      IF (ICOMP .EQ. 0) THEN
  249. *---- Cas d'un appel a sespa (matrices symetriques => mode Reel)
  250.       IF (INSYM .EQ. 0) THEN
  251.          DO 200 IB200 =1, NBMOD
  252.             FREQPR = IPLVAR.PROG(IB200)
  253.             CALL W2FREQ ( FREQPR, W2, OMEGA2, FREQP1, LIMAGE )
  254.             IF ( IERR .NE. 0 ) RETURN
  255.             IPLVAR.PROG(IB200) = FREQP1
  256.  200     CONTINUE
  257.  
  258. *---- Cas d'un appel a sespaC (matrices nonsymetriques => mode Reel ou Complexe)
  259.       ELSE
  260.          SEGACT IPLVAI*MOD
  261. *         DO 400 IB400 = 1, NBMOD2
  262.          DO 400 IB400 = 1, NBMOD2
  263.             FREQPR = IPLVAR.PROG(IB400)
  264.             FREQPI = IPLVAI.PROG(IB400)
  265.             CALL W2FRQC ( FREQPR, FREQPI, W2,XRW2,XIW2, FREQP1, FREQP2)
  266.             IF ( IERR .NE. 0 ) RETURN
  267. *            IPLVAR.PROG(IB400) = FREQP1
  268. *            IPLVAI.PROG(IB400) = FREQP2
  269. *           on se débarasse des erreurs d'arrondis de W2FRQC
  270.             if (FREQPI .eq. 0.) then
  271.               if ( (FREQPR + W2) .lt. 0.) then
  272.                 IPLVAR.PROG(IB400)= 0.D0
  273.                 IPLVAI.PROG(IB400)= FREQP2
  274.               else
  275.                 IPLVAR.PROG(IB400)= FREQP1
  276.                 IPLVAI.PROG(IB400)= 0.D0
  277.               endif
  278.             else
  279.               IPLVAR.PROG(IB400)= FREQP1
  280.               IPLVAI.PROG(IB400)= FREQP2
  281.             endif
  282. *           WRITE(6,*) IB400,' FREQUENCE =',FREQP1,'+i',FREQP2,
  283. *     &     '  soit',(IPLVAR.PROG(IB400)),'+i',(IPLVAI.PROG(IB400))
  284.  400     CONTINUE
  285.       END IF
  286.  
  287. *     -- ON GARDE LES NBMOD2 PREMIERS ELEMENTS PROPRES --
  288. *             -- ET ON DETRUIT LES AUTRES --
  289.  
  290. *     -- LES VALEURS PROPRES --
  291.  
  292. *      JG = NBMOD
  293.       JG = NBMOD2
  294.       SEGADJ IPLVAR
  295.       SEGDES IPLVAR
  296.  
  297.       IF (INSYM .NE. 0) THEN
  298. *         JG = NBMOD
  299.          JG = NBMOD2
  300.          SEGADJ IPLVAI
  301.          SEGDES IPLVAI
  302.       ENDIF
  303.  
  304. *     -- LES MODES PROPRES --
  305.       SEGACT IPLVER
  306. *     si le dernier mode est complexe, alors il nous faut le (nbmod + 1)^ieme
  307. *     qui correspond a la partie imaginaire indissociable de la partie reelle,
  308. *     ce qui justifie l'introduction de nbmod2
  309.       DO 550 IB550 = (NBMOD2 + 1), NBVEC
  310.          IPCHPO = IPLVER.ICHPOI(IB550)
  311.          CALL DTCHPO ( IPCHPO )
  312.  550  CONTINUE
  313.       N1 = NBMOD2
  314.       SEGADJ IPLVER
  315.       SEGDES IPLVER
  316.  
  317.  
  318. *     --  CLASSEMENT DES MODES PROPRES --
  319.  
  320.       IF (INSYM .eq. 0) THEN
  321. c        on n'a retenu que le nombre de mode demandes,
  322. *        on trie par lambda croissant de maniere a avoir le bon numero
  323.          CALL ORDVEC ( IPLVAR, IPLVER,.false. )
  324.       ENDIF
  325. * rem: mode complexes deja classés par module de la val propre croissant
  326.  
  327.  
  328. *     -- AFFICHAGE --
  329.  
  330.       IF ( IIMPI .EQ. 2 ) THEN
  331.          WRITE ( IOIMP, 1000 ) FREQ
  332.  1000    FORMAT( /8X, 'SOLUTION POUR LA FREQUENCE:', E12.6,
  333.      1        /8X, '---------------------------',  /   )
  334.       ENDIF
  335.  
  336.  
  337. *     -- CONSTRUCTION DE LA SOLUTION --
  338.  
  339.       IF (INSYM .EQ. 0) THEN
  340.          CALL CRSOLU (W2,IPLVAR,IPLVER,NBMOD,IPKW2M,IPMASS,IPMODE)
  341.       ELSE
  342.          CALL CRSOL1(W2,IPLVAR,IPLVAI,IPLVER,NBMOD,NBMOD2,
  343.      &   IPKW2M,IPMASS,MTAB3,I,INF0)
  344.       END IF
  345.  
  346.  
  347. *     -- SUPPRESSION DES OBJETS DE TRAVAIL --
  348.  
  349.       CALL DTRIGI ( IPKW2M )
  350.       IF ( IERR .NE. 0 ) RETURN
  351.       CALL DTLCHP ( IPLVER )
  352.       IF ( IERR .NE. 0 ) RETURN
  353.       CALL DTLREE ( IPLVAR )
  354.       IF ( IERR .NE. 0 ) RETURN
  355.       IF (INSYM .NE. 0) CALL DTLREE ( IPLVAI )
  356.  
  357.       RETURN
  358.       END
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales