Télécharger qzbasc.eso

Retour à la liste

Numérotation des lignes :

  1. C QZBASC SOURCE PV 20/03/30 21:23:18 10567
  2. SUBROUTINE QZBASC (ALFR, ALFI, BETA, MATZ, MELEME, MCOMP,
  3. & IPBC,ERR,NWANTED)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. *
  7. ************************************************************************
  8. * CREATION D'UNE BASE DE MODES PROPRES COMPLEXES POUR VIBC
  9. * ______________________________________________________________________
  10. *
  11. * AUTEUR : Nicolas BENECH, 13 avril 1995
  12. * MODIF : BP, 2016-01-15 on reecrit tout
  13. * ______________________________________________________________________
  14. *
  15. * MODULE(S) APPELANT(S) : VIBRAC
  16. *
  17. * MODULE(S) APPELE(S) : CRTABL, ACCTAB, ECCTAB, CREPO1, MUCHPO
  18. * ______________________________________________________________________
  19. *
  20. * EN ENTREE :
  21. * -ALFR : partie reelle des valeurs propres
  22. * -ALFI : partie imaginaire des valeurs propres
  23. * -BETA : denominateur (reel) des valeurs propres
  24. * (denominateur nul --> valeur propre infinie)
  25. * -MATZ : vecteurs propres complexes
  26. * -MELEME : maillage support des chpoints
  27. * -MCOMP : liste des composantes des chpoints
  28. * ______________________________________________________________________
  29. *
  30. * EN SORTIE :
  31. * -IPBC : base de modes complexes
  32. * ______________________________________________________________________
  33. *
  34. * REMARQUE : la variable SEUIL definie au début des données générales
  35. * permet d'identifier les valeurs propres infinies, et
  36. * d'annuler le mode correspondant.
  37. * Sa valeur a ete fixee a 100 x EPSLON(1.0D0) pour le moment
  38. ************************************************************************
  39. *
  40. -INC CCREEL
  41. *-
  42.  
  43. -INC PPARAM
  44. -INC CCOPTIO
  45. -INC SMELEME
  46. -INC SMLCHPO
  47. -INC SMLMOTS
  48. -INC SMCHPOI
  49. -INC SMRIGID
  50. -INC SMTABLE
  51. -INC SMLREEL
  52. -INC SMLENTI
  53. -INC SMCOORD
  54. * segment de travail = MTRAV modifie pour NNCHPO chpoints
  55. SEGMENT MTRAV
  56. CHARACTER*4 INCO(NNIN)
  57. REAL*8 BB(NNCHPO,NNIN,NNNOE)
  58. INTEGER IBIN(NNIN,NNNOE),IGEO(NNNOE),NHAR(NNIN)
  59. ENDSEGMENT
  60. *
  61. * ******** INCO(NNIN) CONTIENT LE NOMS DES NNIN INCONNUES DIFFERENTES.
  62. *
  63. * ******** BB(I,J) EST LA VALEUR DE LA IEME INCONNUE DU CHAMP POUR
  64. * ******** LE JEME NOEUD DU TABLEAU IGEO.
  65. *
  66. * ******** IBIN(I,J)=1 OU 0. 1 INDIQUE QUE LA I EME INCONNUE DU CHAMP
  67. * ******** EXISTE POUR LE J EME NOEUD DU TABLEAU IGEO.
  68. *
  69. * ******** IGEO(I) EST LE NUMERO A METTRE DANS UN OBJET MELEME POUR
  70. * ******** REFERENCER LE IEME NOEUD
  71. *
  72. * ******** NHAR(I) EST LE NUMERO D'HARMONIQUE SI CALCUL AXI OU
  73. * ******** SIGNIFIE CONTRAINTE PLANE,DEFORMATION PLANE OU DEF PLAN GEN
  74. *
  75. SEGMENT ICPR(nbpts)
  76. *
  77. REAL*8 XVAL,SEUIL
  78. INTEGER I, J, K, IC, NUMAFF,ERR
  79. LOGICAL MODANN, MODREL, AFFICH,CONV
  80. CHARACTER*4 NOMDDL
  81. *
  82. POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL
  83. POINTEUR MATZ.XMATRI
  84. POINTEUR MCOMP.MLMOTS
  85. *
  86.  
  87. ************************************************************************
  88. * DONNEES GENERALES *
  89. ************************************************************************
  90.  
  91. * Ecriture des messages pour verification
  92.  
  93. AFFICH = IIMPI.GE.21
  94. c AFFICH = .TRUE.
  95. *
  96. IF (AFFICH)
  97. & WRITE (IOIMP,*) 'QZBASC: Extraction des donnees generales...'
  98. *
  99. * seuil pour le denominateur d'une valeur propre
  100.  
  101. * SEUIL = (EPSLON(1.0D0)*100)
  102. SEUIL = 1.D-99
  103.  
  104. * nombre de modes calcules
  105. NBMOD1=ALFR.PROG(/1)
  106.  
  107.  
  108. **** si on souhaite un nombre de modes inferieur, on devra alors trier
  109.  
  110. IF(NWANTED.GT.0.AND.NWANTED.LT.NBMOD1) THEN
  111. c tri selon 1/module**2 : il faut le calculer explicitement
  112. SEGINI,MLREEL=BETA
  113. I=1
  114. 88 IF(I.GT.NBMOD1) GOTO 22
  115. PROG(I) = (PROG(I)**2) / (ALFI.PROG(I)**2 + ALFR.PROG(I)**2)
  116. c bp : si c'est un complexe, on force l egalite du module
  117. c comme le tri d'ORDON1 est stable, les complexes i et i+1
  118. c resteront a des indices successifs
  119. IF(ABS(ALFI.PROG(I)).NE. 0.0D0) THEN
  120. I=I+1
  121. PROG(I)=PROG(I-1)
  122. ENDIF
  123. I=I+1
  124. GOTO 88
  125. 22 CONTINUE
  126. IORDRE=1
  127. CALL ORDON1(MLREEL,.FALSE.,.TRUE.,IORDRE)
  128. segsup,MLREEL
  129. MLENTI=IORDRE
  130. segact,MLENTI
  131. ELSE
  132. IORDRE=0
  133. NWANTED=NBMOD1
  134. ENDIF
  135.  
  136.  
  137. **** preparation des CHPOINTs deformees *******************************
  138.  
  139. * nombre de ddls (=dimension de MELEME = dimension de MCOMP)
  140. SEGACT,MELEME,MCOMP
  141. NDDL=MELEME.NUM(/2)
  142.  
  143. * on cree ICPR
  144. SEGINI,ICPR
  145. * on dimensionne au maxi MTRAV
  146. NNIN=NDDL
  147. NNNOE=NDDL
  148. NNCHPO=NBMOD1
  149. SEGINI,MTRAV
  150. NNIN=0
  151. NNNOE=0
  152.  
  153. *---- traitement de chaque ddl ----------------------------------
  154. DO 100 J=1,NDDL
  155.  
  156. c - NOEUD
  157. IP=NUM(1,J)
  158. JNOE=ICPR(IP)
  159. c nouveau noeud #IP de numero local #JNOE
  160. IF(JNOE.EQ.0) THEN
  161. NNNOE=NNNOE+1
  162. JNOE=NNNOE
  163. ICPR(IP)=JNOE
  164. IGEO(JNOE)=IP
  165. ENDIF
  166.  
  167. c - COMPOSANTE
  168. NOMDDL=MCOMP.MOTS(J)
  169. IF(NNIN.EQ.0) GOTO 111
  170. DO 110 JIN=1,NNIN
  171. IF(NOMDDL.EQ.INCO(JIN)) GOTO 112
  172. c NOMDDL trouvee dans INCO(JIN)
  173. 110 CONTINUE
  174.  
  175. 111 CONTINUE
  176. c NOMDDL pas trouvee dans INCO(JIN) -> on l'ajoute a la fin
  177. NNIN=NNIN+1
  178. JIN=NNIN
  179. INCO(JIN)=NOMDDL
  180. c bp : pour l'instant on laisse NHAR=NIFOUR, mais il faudrait
  181. c recuperer IRIGEL(5,:) des rigidites d'entree ou autre ...?
  182. NHAR(JIN)=NIFOUR
  183.  
  184. 112 CONTINUE
  185. c NOMDDL trouvee dans INCO(JIN) et noeud #IP dans IGEO(JNOE)
  186.  
  187. c Remplissage de IBIN et BB
  188. IBIN(JIN,JNOE)=J
  189. DO I=1,NBMOD1
  190. BB(I,JIN,JNOE)=MATZ.RE(J,I,1)
  191. ENDDO
  192.  
  193. 100 CONTINUE
  194. *---- fin de la boucle sur les ddl ----------------------------------
  195.  
  196. SEGSUP,MCOMP,ICPR
  197. SEGDES,MELEME
  198. SEGADJ,MTRAV
  199. c write(*,*) 'IGEO=',(IGEO(iou),iou=1,NNNOE)
  200. c write(*,*) 'INCO=',(INCO(iou),iou=1,NNIN)
  201. c on va generer N1 chpoints en 1 passage dans une copie de CRECHP
  202. N1=NBMOD1
  203. SEGINI,MLCHPO
  204. CALL CRECH4(MTRAV,MLCHPO)
  205. SEGSUP,MTRAV
  206.  
  207.  
  208. *******************************************
  209. * Creation de la table BASE_DE_MODES *
  210. *******************************************
  211. *
  212. IF (AFFICH) WRITE(IOIMP,*) 'Creation de la table BASE_DE_MODES...'
  213. *
  214. CALL CRTABL(IPTAB2)
  215. CALL ECCTAB(IPTAB2,'MOT',0,0.0D0,'SOUSTYPE',.TRUE.,0,
  216. & 'MOT',0,0.0D0,'BASE_DE_MODES',.TRUE.,0)
  217. CALL ECCTAB(IPTAB2,'MOT',0,0.0D0,'MAILLAGE',.TRUE.,0,
  218. & 'MAILLAGE',0,0.0D0,' ',.TRUE.,MELEME)
  219.  
  220.  
  221. ************************************************************************
  222. * BOUCLE SUR LES MODES *
  223. ************************************************************************
  224. *
  225. III=1
  226. 80 IF (III .GT. NWANTED) GOTO 20
  227.  
  228. * selection des NWANTED plus petits modes
  229. IF(IORDRE.NE.0) THEN
  230. I = LECT(III)
  231. ELSE
  232. I = III
  233. ENDIF
  234. *
  235. MODREL = (ABS(ALFI.PROG(I)).EQ. 0.0D0)
  236. IF (AFFICH) THEN
  237. IF (MODREL) THEN
  238. WRITE (IOIMP,*) III,': MODE ',I,' REEL',ALFI.PROG(I),
  239. & '-i',ALFR.PROG(I),' / 2pi',BETA.PROG(I)
  240. ELSE
  241. WRITE (IOIMP,*) III,': MODE ',I,' COMPLEXE',ALFI.PROG(I),
  242. & '-i',ALFR.PROG(I),' / 2pi',BETA.PROG(I)
  243. ENDIF
  244. ENDIF
  245. *
  246. * ---- frequence infinie ?
  247. *
  248. MODANN = (ABS(BETA.PROG(I)).LT.SEUIL)
  249.  
  250. * ---- oui
  251. IF (MODANN) THEN
  252. 90 WRITE (IOIMP,*) 'Attention !!! Mode ',I,
  253. & ' annule : frequence infinie.'
  254. CALL CRTABL(MTAB1)
  255. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'SOUSTYPE', .TRUE.,0,
  256. & 'MOT',0,0.0D0,'MODE_ANNULE',.TRUE.,0)
  257. CALL ECCTAB(IPTAB2,'ENTIER',I,0.0D0,' ',.TRUE.,0,
  258. & 'TABLE', 0, 0.0D0,' ',.TRUE.,MTAB1)
  259. IF (.NOT. MODREL) THEN
  260. MODREL = .TRUE.
  261. III = III + 1
  262. GO TO 90
  263. ENDIF
  264.  
  265. * ---- non
  266. ELSE
  267. *
  268. *------- Deformees reelle (1) et imaginaire (2) du mode i
  269.  
  270. * Recup des CHPOINTs :
  271. MCHPO1=ICHPOI(I)
  272.  
  273. * Cas reel : il faut un chpoint nul -> on crée un chpoint vide
  274. IF (MODREL) THEN
  275. NAT=1
  276. NSOUPO=0
  277. SEGINI,MCHPO2
  278. MCHPO2.IFOPOI=IFOMOD
  279. MCHPO2.MOCHDE='CHPOINT CREE PAR QZBASC'
  280. * Cas complexe : on fait la Recup
  281. ELSE
  282. MCHPO2=ICHPOI(I+1)
  283. ENDIF
  284. c write(*,*) '>>> mode',I,'phiR=',MCHPO1,'phiI=',MCHPO2
  285. *
  286. *
  287. ************************************************
  288. * Creation de la table MODE *
  289. ************************************************
  290. *
  291. *----- valeur propre (reelle ou complexe)
  292. *
  293. * Attention : on enregistre f = w/2pi
  294. * avec iw = landa donc w=-i landa !!!
  295. *
  296. IF (AFFICH) WRITE (*,*) 'Construction de la table MODE ...'
  297. *
  298. CALL CRTABL(MTAB1)
  299. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'SOUSTYPE',.TRUE.,0,
  300. & 'MOT',0,0.0D0,'MODE_COMPLEXE',.TRUE.,0)
  301. *
  302. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'NUMERO_MODE',.TRUE.,0,
  303. & 'ENTIER',I,0.0D0,' ',.TRUE.,0)
  304. CALL CREPO1(0.0D0,0.0D0,0.0D0,IPOIN)
  305. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  306. & 'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  307. XVAL=(-1*ALFR.PROG(I))/(2.*XPI*BETA.PROG(I))
  308. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'FREQUENCE_IMAGINAIRE',.TRUE.,0
  309. $ ,'FLOTTANT',0,XVAL,' ',.TRUE.,0)
  310. XVAL=(ALFI.PROG(I))/(2.*XPI*BETA.PROG(I))
  311. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'FREQUENCE_REELLE',.TRUE.,0,
  312. & 'FLOTTANT',0,XVAL,' ',.TRUE.,0)
  313. *
  314. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'DEFORMEE_MODALE_REELLE',
  315. & .TRUE.,0,'CHPOINT',0,0.0D0,' ',.TRUE.,MCHPO1)
  316. CALL ECCTAB(MTAB1,'MOT',0,0.0D0,'DEFORMEE_MODALE_IMAGINAIRE',
  317. & .TRUE.,0,'CHPOINT',0,0.0D0,' ',.TRUE.,MCHPO2)
  318. *
  319. CALL ECCTAB(IPTAB2,'ENTIER',III,0.0D0,' ',.TRUE.,0,'TABLE',
  320. & 0,0.0D0,' ',.TRUE.,MTAB1)
  321. *
  322. IF (MODREL) GOTO 70
  323. *
  324. *----- Valeur propre complexe conjuguee
  325. *
  326. *
  327. SEGINI, MTAB2=MTAB1
  328. SEGDES, MTAB2
  329. III = III + 1
  330. IF (AFFICH) WRITE (*,*) 'Mode conjugue, no : ',I
  331. CALL ECCTAB(MTAB2,'MOT',0,0.0D0,'NUMERO_MODE',.TRUE.,0,
  332. & 'ENTIER',I,0.0D0,' ',.TRUE.,0)
  333. CALL ECCTAB(MTAB2,'MOT',0,0.0D0,'FREQUENCE_REELLE',.TRUE.,0,
  334. & 'FLOTTANT',0,-1.*XVAL,' ',.TRUE.,0)
  335. XVAL=-1.D0
  336. CALL MUCHPO(MCHPO2,XVAL,IRET,1)
  337. CALL ECCTAB(MTAB2,'MOT',0,0.0D0,'DEFORMEE_MODALE_REELLE',
  338. & .TRUE.,0,'CHPOINT',0,0.0D0,' ',.TRUE.,MCHPO1)
  339. CALL ECCTAB(MTAB2,'MOT',0,0.0D0,'DEFORMEE_MODALE_IMAGINAIRE',
  340. & .TRUE.,0,'CHPOINT',0,0.0D0,' ',.TRUE.,IRET)
  341. *
  342. CALL ECCTAB(IPTAB2,'ENTIER',III,0.0D0,' ',.TRUE.,0,'TABLE',
  343. & 0,0.0D0,' ',.TRUE.,MTAB2)
  344. *
  345. IF (AFFICH) WRITE (*,*) 'Construction du conjugue ok'
  346. 70 CONTINUE
  347. ENDIF
  348.  
  349.  
  350. III = III + 1
  351. GOTO 80
  352. ************************************************************************
  353. * FIN DE BOUCLE SUR LES MODES *
  354. ************************************************************************
  355. 20 CONTINUE
  356.  
  357.  
  358. ***** MENAGE ***********************************************************
  359.  
  360. SEGSUP, ALFR, ALFI, BETA
  361. SEGSUP, MATZ
  362. IF(IORDRE.GT.0) SEGSUP,MLENTI
  363.  
  364. *
  365. ********************************************
  366. * Creation de la table BASE_MODALE *
  367. ********************************************
  368. *
  369. IF (AFFICH) WRITE (*,*) 'Creation de la table BASE_MODALE...'
  370. *
  371. CALL CRTABL(IPBC)
  372. CALL ECCTAB(IPBC,'MOT',0,0.0D0,'SOUSTYPE',.TRUE.,0,
  373. & 'MOT',0,0.0D0,'BASE_MODALE',.TRUE.,0)
  374. CALL ECCTAB(IPBC,'MOT',0,0.0D0,'MODES',.TRUE.,0,
  375. & 'TABLE',0,0.0D0,' ',.TRUE.,IPTAB2)
  376. cbp : on ajoute un indicateur de succes (en plus du message dans vibrac)
  377. CONV=ERR.eq.0
  378. CALL ECCTAB(IPBC,'MOT',0,0.0D0,'CONVERGENCE',.TRUE.,0,
  379. & 'LOGIQUE ',0,0.0D0,' ',CONV,0)
  380. *
  381. RETURN
  382. END
  383.  
  384.  
  385.  

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