Télécharger la plaquette Cast3M


fr en

Télécharger qzbasc.eso

Retour à la liste

Numérotation des lignes : 

qzbasc

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

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