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