Télécharger la plaquette Cast3M


fr en

Télécharger hbmalo.eso

Retour à la liste

Numérotation des lignes : 

hbmalo

  1. C HBMALO    SOURCE    CB215821  24/04/12    21:16:13     11897          
  2. c
  3. C DEVALO    SOURCE    BP208322  18/01/30    21:15:12     9719           
  4.       SUBROUTINE HBMALO(ITBAS,ITKM,ITA,ITLIA,ITCHAR,ITINIT,NINS,ITREDU,
  5.      &                  IPARNUM,KPREF,KTQ,KTKAM,KTPHI,KTLIAA,KTEMP,
  6.      &                  KTLIAB,KTFEX,KTPAS,KTRES,KTNUM,IPMAIL,REPRIS,
  7.      &                  KPARNUM,KSORT,ICHAIN,KOCLFA,KOCLB1,NHBM,NFFT)
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8(A-H,O-Z) 
  10. *--------------------------------------------------------------------*
  11. *                                                                    *
  12. *     Operateur DYNC : continuation par longueur d'arc               *
  13. *     ________________________________________________               *
  14. *                                                                    *
  15. *     Dimensionnement des tableaux de travail ( allocation de la     *
  16. *     memoire ).                                                     *
  17. *                                                                    *
  18. *     Parametres:                                                    *
  19. *                                                                    *
  20. * e   ITBAS   Table representant une base modale                     *
  21. * e   ITKM    Table contenant les matrices XK et XM                  *
  22. * e   ITA     Table contenant la matrice XASM                        *
  23. * e   ITLIA   Table rassemblant la description des liaisons          *
  24. * e   ITCHAR  Table contenant les chargements                        *
  25. * e   ITINIT  Table donnant les conditions initiales                 *
  26. * e   NINS    On veut une sortie tous les NINS pas de calcul         *
  27. * e   ITREDU  Table contenant les noms d'inconnues de la base B      *
  28. *             auxquelles on se restreint                             *  
  29. * e   IPARNUM Table des les parametres numeriques (continuation)     *  
  30. * e   KPREF   Segment des points de reference                        *
  31. * e   NHBM    Nombre d'harmoniques de l'approximation                *
  32. * s   MTQ     Segment contenant les coefficients de Fourier          *
  33. * s   MTKAM   Segment contenant les vecteurs XK, XASM et XM          *
  34. * s   MTPHI   Segment contenant les deformees modales                *
  35. * s   MTLIAA  Segment descriptif des liaisons en base A              *
  36. * s   MTLIAB  Segment descriptif des liaisons en base B              *
  37. * s   MTFEX   Segment contenant les chargements libres               *
  38. * s   MTPAS   Segment des variables au cours d'un pas de temps       *
  39. * s   PARNUM  Segment des parametres numeriques (continuation)       *
  40. * s   MTRES   Segment de sauvegarde des resultats                    *
  41. * s   MTNUM   Segment contenant les parametres temporels             *
  42. * s   REPRIS  Vrai si reprise de calcul, faux sinon                  *
  43. * s   ICHAIN  Segment MLENTI (ACTIF) contenant les adresses des      *
  44. *             chaines dans la pile des mots de CCNOYAU               *
  45. * s   KOCLFA  Segment contenant les tableaux locaux de la subroutine *
  46. *             DEVLFA                                                 *
  47. * s   KOCLB1  Segment contenant les tableaux locaux de la subroutine *
  48. *             DEVLB1                                                 *
  49. *                                                                    *
  50. *     Auteur, date de creation:                                      *
  51. *                                                                    *
  52. *     Roberto ALCORTA, le 18 juin 2019.                              *
  53. *--------------------------------------------------------------------*
  54. -INC PPARAM
  55. -INC CCOPTIO
  56. -INC SMCOORD
  57. -INC SMMODEL
  58. -INC SMELEME
  59. -INC SMCHAML
  60. -INC SMLENTI
  61. -INC CCREEL
  62. -INC SMCHPOI
  63. *-INC TMDYNC.INC
  64. ************************** debut TMDYNC.INC ****************************
  65.  
  66. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  67. *    TODO   : a extraire dans un include des que stabilise
  68. *
  69. *     Segment des variables generalisees:
  70. *     -----------------------------------
  71.       SEGMENT MTQ
  72.          REAL*8 Q1(NT1)
  73.          REAL*8 OMEG,XPARA
  74.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  75.          REAL*8 dX(NT1), dw, dv
  76.       ENDSEGMENT
  77. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  78. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  79. *            avec q_i vecteur de dimension n ou n=nombre de modes
  80. *     OMEG : frequence fondamentale de l'approximation
  81. *     XPARA: parametre de continuation (par defaut la frequence) 
  82. *            \in [PARINI,PARFIN]
  83. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  84. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  85. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  86. *            lineaires et constantes
  87. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  88. *
  89. *
  90. *     Segment contenant les matrices XK, XASM et XM:
  91. *     ---------------------------------------------
  92.       SEGMENT MTKAM
  93.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  94.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  95. *         REAL*8 GAMFIN(NPC2,nl1)
  96.       ENDSEGMENT
  97. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  98. *     GAM et IGAM : matrices pour la FFT et son inverse
  99. *     GAMFIN      : 
  100. *
  101. *     Segment des deformees modales:
  102. *     ------------------------------
  103. *     (idem DYNE)
  104.       SEGMENT MTPHI
  105.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  106.          INTEGER IAROTA(NSB)
  107.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  108.       ENDSEGMENT
  109. *
  110. *     Segment descriptif des liaisons en base A:
  111. *     ------------------------------------------
  112. *     (idem DYNE)
  113.       SEGMENT MTLIAA
  114.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  115.          REAL*8  XPALA(NLIAA,NXPALA)
  116.       ENDSEGMENT
  117. *
  118. *     Segment descriptif des liaisons en base B:
  119. *     ------------------------------------------
  120. *     (idem DYNE)
  121.       SEGMENT MTLIAB
  122.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  123.          REAL*8  XPALB(NLIAB,NXPALB)
  124.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  125.       ENDSEGMENT
  126. *
  127. *     Segment representant les chargements exterieurs:
  128. *     -----------------------------------------------
  129.       SEGMENT MTFEX
  130.          REAL*8  FEXA(NT1)
  131.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  132.          INTEGER BAL
  133.       ENDSEGMENT
  134. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  135. *            Fourier et exprimes en base A
  136. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  137. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  138. *             compatibilite avec calcul des Fnl.
  139. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  140. *            (cad proportionnel a OMEG**2)
  141. *
  142. *     Segment "local" pour DEVLFA:
  143. *     ----------------------------
  144.       SEGMENT LOCLFA
  145.          REAL*8 FTEST(NA1,4)
  146.       ENDSEGMENT
  147. *      
  148. *     Segment "local" pour DEVLB1:
  149. *     ----------------------------
  150.       SEGMENT LOCLB1
  151.          REAL*8 FTEST2(NPLB,6)
  152.       ENDSEGMENT
  153. *      
  154. *     Segment contenant les variables au cours d un pas de temps:
  155. *     ----------------------------------------------------------
  156.       SEGMENT MTPAS
  157.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  158.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  159.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  160.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  161.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  162.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  163.       ENDSEGMENT
  164. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  165. *     XPTB       : deplacement du point d'une liaison en base B
  166. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  167. *     FEXB       : forces exterieures en base B (a priori uniquement 
  168. *                  pour les moments appliques aux rotations rigides ?)
  169. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  170. *                  maillage de contact dans certaines liaisons)
  171. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  172. *                       base A/B (= contributions a dFnl/dX)
  173. *
  174. *
  175. *     Segment des points de reference des modes (base A):
  176. *     --------------------------------------------------
  177.       SEGMENT MPREF
  178.          INTEGER IPOREF(NPREF)
  179.       ENDSEGMENT
  180. *
  181. *     Segment des points en base B:
  182. *     -----------------------------
  183.       SEGMENT NCPR(NBPTS)
  184. *     NCRP(#global) = #local dans XPTB (1er indice)      
  185. *
  186. *     Segment des parametres numeriques pour la continuation: 
  187. *     ------------------------------------------------------
  188.       SEGMENT PARNUM
  189.          CHARACTER*4 TYPS
  190.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  191.          REAL*8   PARINI,PARFIN
  192.          INTEGER  ITERMAX,NBPAS
  193.          LOGICAL  JANAL
  194.       ENDSEGMENT
  195. *
  196. *     Segment des resultats:
  197. *     ---------------------
  198.       SEGMENT  PSORT
  199.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  200.         REAL*8 VSAVE(NPAS)
  201.         LOGICAL ZSAVE(NPAS)
  202.         CHARACTER*2 TYPBIF(NBIFU)
  203.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  204.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  205.         INTEGER CBIF
  206.       ENDSEGMENT
  207. *     QSAVE(i,j) = Q harmonique i au pas j
  208. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  209. *     ZSAVE(j)   = stabilite au j-eme pas
  210. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  211. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  212. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  213. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  214. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  215. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  216.  
  217. *     Segment des tableaux de travail:
  218. *     -------------------------------
  219.       SEGMENT  MTEMP
  220.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  221.          REAL*8 T02(NT1+2), TP2(NT1+2)
  222.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  223.          REAL*8 res
  224.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  225.          REAL*8 QOLD(NT1),OMEGOLD
  226.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  227.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  228.       ENDSEGMENT
  229. *     Jacobiennes augmentees 
  230. *     Ja : [ RX  Rw ; dX dw]
  231. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  232.  
  233. *      SEGMENT  NNNN
  234. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  235. *      ENDSEGMENT
  236.  
  237. *************************** fin TMDYNC.INC *****************************
  238.  
  239. *
  240. *     Segment de travail pour reprise force POLYNOMIALE base A
  241. *     a metrre dans TMDYNC ?
  242.       SEGMENT,MTRA
  243.         INTEGER IPLA(NTRA)
  244.       ENDSEGMENT
  245. *      
  246.       LOGICAL L0,L1,REPRIS
  247.       CHARACTER*6  MO2
  248.       CHARACTER*8  TYPRET,CHARRE,CHARR0
  249.       CHARACTER*10 MO1
  250.       CHARACTER*40 MONMOT
  251.  
  252.       ITREP  = 0
  253.       MTQ    = 0
  254.       MTKAM  = 0
  255.       MTPHI  = 0
  256.       MTLIAA = 0
  257.       MTEMP  = 0
  258.       MTLIAB = 0
  259.       MTFEX  = 0
  260.       MTPAS  = 0
  261.       MTNUM  = 0
  262.       MTRA   = 0
  263.       PARNUM = 0
  264.       PSORT  = 0
  265.       XTINI  = 0.D0
  266.       ITLA = 0
  267.       ITLB = 0
  268. c       NNNN = 0
  269.       REPRIS = .FALSE.
  270.  
  271.  
  272. ************************************************************************
  273. *     Recherche du nombre de modes: autant que de points de reference
  274. ************************************************************************
  275.       IF (IIMPI.EQ.333) write(IOIMP,*) 'HBMALO: recup des modes' 
  276. *
  277.       MPREF = KPREF
  278.       NPREF = IPOREF(/1)
  279.       NA1   = NPREF
  280.  
  281.  
  282. c     on intialise NB1 a 1; le segment sera eventuellement ajuste
  283. c     lors du remplissage par D2VTRA
  284.       NB1  = 1
  285.       NB1K = 1
  286.       NB1C = 1
  287.       NB1M = 1
  288.       nl1 = 2*NHBM+1
  289.       NT1 = nl1*na1
  290. ** FORMULE INTELLIGENTE A TROUVER
  291.       NPC1 = NFFT
  292.       NPC2 = 40*NFFT
  293.  
  294. **===============================
  295.       SEGINI,MTQ,MTKAM
  296.       SEGINI,LOCLFA
  297.       SEGINI,MTEMP
  298. c       SEGINI,NNNN
  299.       KTEMP = MTEMP
  300.       KOCLFA = LOCLFA
  301.       KTQ   = MTQ
  302.       KTKAM = MTKAM  
  303.  
  304. *     Transformees de Fourier Inverse/Directe
  305. *     Pour la continuation: GAM,IGAM 
  306.       CALL AFTM(NPC1,NHBM,GAM,IGAM,DL)
  307. **     Pour la stabilite: GAMFIN si besoin de + de points
  308. *      CALL AFTM(NPC2,NHBM,GAMFIN,IGAM2,DL2)
  309. c       SEGDES,NNNN
  310.  
  311.  
  312. **********************************************************************
  313. *     Segment des parametres numeriques:
  314. *     on renseigne ici seulement TYPS 
  315. *     d'apres 'TYPE' = FORCe ou AUTOnome ou NonNormalMode?
  316. **********************************************************************
  317.       SEGINI, PARNUM
  318.       KPARNUM = PARNUM
  319.       IF (IPARNUM.gt.0) THEN   
  320.         CALL ACCTAB(IPARNUM,'MOT',I0,X0,'TYPE',L0,IP0,
  321.      &                      'MOT',I1,X1,CHARRE,L1,IP1)    
  322.         IF(IERR.NE.0) RETURN
  323.         TYPS(1:4) = CHARRE(1:4) 
  324.       ELSE 
  325. c       valeur par defaut (pas super)
  326.         TYPS(1:4) = 'FORC'
  327.       ENDIF
  328. c     verif
  329.       IF(TYPS.NE.'FORC' .AND. TYPS.NE.'AUTO' .AND.TYPS.NE.'NNM') THEN 
  330. c         On a lu %m1:8, alors qu'on attend un des mots cles suivant :
  331. c         %m9:16 %m17:24 %m25:32 %m33:40. Consulter la notice.
  332.           MOTERR(1:4)=TYPS(1:4)
  333.           MOTERR(5:8)='    '
  334.           MOTERR(9:40)='FORC    AUTO    NNM     '
  335.           CALL ERREUR(930)
  336.           RETURN      
  337.       ENDIF
  338.  
  339.  
  340. ************************************************************************
  341. *     Table INITIAL --> Initialisation du vecteur des inconnues 
  342. ************************************************************************
  343.  
  344.       IF (ITINIT.GT.0) THEN   
  345.  
  346. ******* Cas d'une reponse FORCee ou d'un systeme AUTOnome *******
  347.         IF (TYPS.EQ.'FORC' .OR. TYPS.EQ.'AUTO') THEN   
  348.  
  349. c         TINIT . jharm . imode  = valeur (flottant)  -> idem sortie 
  350. c         TINIT . jharm  = chpoint                    -> idem chargement
  351. c         boucle sur les harmoniques
  352.           DO KHBM=-1*NHBM,NHBM
  353.             TYPRET = ' '
  354.             CALL ACCTAB(ITINIT,'ENTIER',KHBM,X0,CHARR0,L0,IP0,
  355.      &                  TYPRET,I1,X1,CHARRE,L1,IP1)
  356.             IF(IERR.NE.0) RETURN
  357. *          -cas d'une table de CHPOINT     
  358.             IF(TYPRET.EQ.'CHPOINT ') THEN            
  359.               MCHPOI=IP1
  360.               SEGACT,MCHPOI 
  361.               NSOUPO=IPCHP(/1)
  362. c             boucle sur les zones (definies a partir des noms de composantes)
  363.               DO I=1,NSOUPO
  364.                  MSOUPO = IPCHP(I)
  365.                  SEGACT,MSOUPO
  366. c                on recupere le maillage
  367.                  MELEME = IGEOC
  368.                  SEGACT,MELEME
  369. c                NC = nombre de composantes
  370.                  NC = NOCOMP(/2)
  371.                  IF(NC.NE.1) THEN
  372. c                  Il faut specifier un champ par point avec une seule composante
  373.                    CALL ERREUR(180)
  374.                    RETURN
  375.                  ENDIF
  376.                  MPOVAL = IPOVAL
  377.                  SEGACT,MPOVAL
  378. c                VPOCHA(i,j) = valeur au noeud i de la composante j
  379.                  N = VPOCHA(/1)
  380.                  DO J=1,N
  381. c                   on cherche le #noeud dans IPOREF
  382.                     KNOE = NUM(1,J)
  383.                     CALL PLACE2(IPOREF,NPREF,IPOS,KNOE)
  384. c                   IPOS est le numero (local) de mode
  385.                     IF (IPOS.NE.0) THEN
  386. c                       DO K=1,NC
  387.                         K=1
  388.                           IF (KHBM.LE.0) THEN
  389.                             IND1=2*ABS(KHBM)
  390.                           ELSE
  391.                             IND1=2*ABS(KHBM)-1
  392.                           ENDIF
  393.                           IND=NA1*IND1+IPOS
  394. c                         normalement, on ne peut/doit definir qu'1 fois Q1 initial
  395.                           Q1(IND) = VPOCHA(J,K)
  396. c                       ENDDO
  397.                     ENDIF
  398.                  ENDDO
  399.                  SEGDES,MPOVAL,MELEME,MSOUPO
  400.               ENDDO
  401.               SEGDES,MCHPOI              
  402. c *          -cas d'une table de TABLE     
  403. c             ELSEIF(TYPRET.EQ.'TABLE   ') THEN
  404. c             TODO
  405.             ELSEIF(TYPRET.NE.' ') THEN
  406. c              On ne veut pas d'objet de type %m1:8
  407.                MOTERR(1:8)=TYPRET
  408.                CALL ERREUR(39)
  409.                RETURN               
  410.             ENDIF
  411.           ENDDO 
  412.  
  413. *         frequence initiale        
  414.           CALL ACCTAB(ITINIT,'MOT',I0,X0,'FREQUENCE',L0,IP0,
  415.      &                                 'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  416.           IF(IERR.NE.0) RETURN
  417.           OMEG   = X1
  418.  
  419. ******* Cas d'une reponse FORCee ou d'un systeme AUTOnome *******
  420.         ELSEIF (TYPS.EQ.'NNM') THEN
  421.  
  422. *         mode lineaire considere comme point initial
  423.           CALL ACCTAB(ITINIT,'MOT',I0,X0,'MODE',L0,IP0,
  424.      &                       'ENTIER',I1,X1,CHARRE,L1,IP1)         
  425.           IF(IERR.NE.0) RETURN
  426.           JMODE=I1
  427. c         OMEG sera rempli lors du parcours de la base ....
  428. *         energie initiale du mode (doit etre assez faible ~quasi-nulle)
  429.           CALL ACCTAB(ITINIT,'MOT',I0,X0,'ENERGIE',L0,IP0,
  430.      &                       'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  431.           XPARA  = ABS(X1)
  432.           IF(XPARA.le.XPETIT) THEN
  433. c           La valeur de %M1:8 doit etre positive
  434.             MOTERR(1:8)='ENERGIE '
  435.             CALL ERREUR(549)
  436.             RETURN
  437.           ENDIF
  438. *         Par defaut, on initialise selon la composante cos Wt (les autres=0)
  439. *         (on aurait pu initialiser sin Wt)
  440.           Q1(NA1+JMODE)=SQRT(XPARA)
  441.  
  442.         ELSE 
  443.  
  444.         ENDIF
  445.  
  446.  
  447.       ENDIF
  448.  
  449.  
  450. ************************************************************************
  451. *     Parametres temporels: pas de temps constant
  452. ************************************************************************
  453. *
  454.       NPC=1    
  455.  
  456.  
  457. ************************************************************************
  458. *     Gestion des segments descriptifs des liaisons
  459. ************************************************************************
  460.       IF (IIMPI.EQ.333) write(IOIMP,*) 'HBMALO: recup des liaisons' 
  461. *
  462.       NLIAA  = 0
  463.       NIPALA = 0
  464.       NXPALA = 0
  465.       NPLAA  = 0
  466.       NPLA   = 0
  467.       NLIAB  = 0
  468.       NIPALB = 0
  469.       NXPALB = 0
  470.       NIP    = 0
  471.       NPLBB  = 0
  472.       NPLB   = 0
  473.       NPLSB  = 0
  474.       IDIMB  = 0
  475.       NA2    = 0
  476.       NSB    = 0
  477.       KCPR   = 0
  478.       NTVAR  = 6 + 4 * IDIM
  479. *
  480. *     MTRA indiquera la presence de liaisons POLYNOMIALEs
  481. *     (on suppose un maximum de 100 liaisons en base A)
  482. *+*   passe a 10000 le 28/1/93
  483.       NTRA   = 10000
  484.       SEGINI,MTRA
  485. *
  486.       IF (ITLIA.NE.0) THEN
  487. *
  488.          TYPRET = ' '
  489.          CALL ACCTAB(ITLIA,'MOT',I0,X0,'LIAISON_A',L0,IP0,
  490.      &                    TYPRET,I1,X1,CHARRE,L1,ITLA)
  491.          IF (IERR.NE.0) RETURN
  492. *
  493. *        Liaisons sur la base A : determination des parametres
  494. *
  495.          IF (ITLA.NE.0) THEN
  496.             IF (TYPRET.EQ.'TABLE   ') THEN
  497.                CALL DYNE21(ITLA,PDT,MTRA,KLIAA,KXPALA,KPLAA,KIPALA)
  498. *               WRITE(*,*) 'KXPALA=',KXPALA
  499.                IF (IERR.NE.0) RETURN
  500.                NLIAA  = KLIAA
  501.                NIPALA = KIPALA
  502.                NXPALA = KXPALA
  503.                NPLAA  = KPLAA
  504.                NPLA   = NA1
  505.             ELSE
  506.                CALL ERREUR(492)
  507.                RETURN
  508.             ENDIF
  509.          ENDIF
  510. *
  511. *        Liaisons sur la base B : determination des parametres
  512. *
  513.          TYPRET = ' '
  514.          CALL ACCTAB(ITLIA,'MOT',I0,X0,'LIAISON_B',L0,IP0,
  515.      &                    TYPRET,I1,X1,CHARRE,L1,ITLB)
  516.          IF (IERR.NE.0) RETURN
  517.          IF (ITLB.NE.0) THEN
  518.             IF (TYPRET.EQ.'TABLE   ') THEN
  519.                CALL DYNE22(ITLB,KLIAB,KXPALB,KPLBB,KPLB,KDIMB,KCPR,
  520.      &                     KIPALB,KNIP)
  521. *           WRITE(*,*) 'KXPALB=',KXPALB
  522.                IF (IERR.NE.0) RETURN
  523.                NLIAB  = KLIAB
  524.                NIPALB = KIPALB
  525.                NXPALB = KXPALB
  526.                NPLBB  = KPLBB
  527.                NPLB   = KPLB
  528.                IDIMB  = KDIMB
  529.                NIP = KNIP
  530.             ELSE
  531.                CALL ERREUR(493)
  532.                RETURN
  533.             ENDIF
  534.          ENDIF
  535.       ENDIF       
  536.       SEGINI,LOCLB1
  537.       KOCLB1 = LOCLB1
  538. *
  539. *     Les segments seront remplis dans le s-p DEVLIA:
  540. *
  541.       SEGINI,MTLIAA
  542.       SEGINI,MTLIAB
  543.       KTLIAA = MTLIAA
  544.       KTLIAB = MTLIAB
  545.       IF (NLIAB.NE.0) THEN
  546.          NCPR = KCPR
  547.          IN = 0
  548.          DO I = 1,NBPTS
  549.             IF (NCPR(I).NE.0) THEN
  550.                IN = IN + 1
  551.                JPLIB(IN) = I
  552.             ENDIF
  553.          ENDDO
  554.          SEGSUP,NCPR
  555.       ENDIF
  556.  
  557.  
  558. ************************************************************************
  559. *     Segment des deformees modales:
  560. ************************************************************************
  561.       IF (IIMPI.EQ.333) write(IOIMP,*) 'HBMALO: recup deformees modales' 
  562.  
  563. ***** Cas table BASE_MODALE et RAIDEUR_ET_MASSE *****
  564.  
  565.       IF (ITKM.GT.0) THEN
  566.          TYPRET = ' '
  567.          CALL ACCTAB(ITKM,'MOT',I0,X0,'BASE_MODALE',L0,IP0,
  568.      &               TYPRET,I1,X1,CHARRE,L1,ITBAS2)
  569.       ELSE
  570.          ITBAS2=ITBAS
  571.       ENDIF
  572.  
  573.       IF (ITBAS2.NE.0) THEN
  574.         CALL ACCTAB(ITBAS2,'MOT',I0,X0,'SOUSTYPE',L0,IP0,
  575.      &                     'MOT',I1,X1,MONMOT,L1,IP1)
  576.         IF (IERR.NE.0) RETURN
  577.         IF (IIMPI.EQ.333) write(IOIMP,*) ITBAS2,'de SOUSTYPE',MONMOT
  578. *
  579. *      -Cas ou la base est unique
  580.         IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  581.           NSB = 1
  582.           NA2 = NA1
  583. *         changement de dimension en cas de corps rigide
  584.           CALL ACCTAB(ITBAS2,'MOT',I0,X0,'MODES',L0,IP0,
  585.      &                    'TABLE',I1,X1,' ',L1,IBAS)
  586.           IP = 0
  587.  22       CONTINUE
  588.             IP = IP + 1
  589.             TYPRET = ' '
  590.             CALL ACCTAB(IBAS,'ENTIER',IP,X0,' ',L0,IP0,
  591.      &                        TYPRET,I1,X1,CHARRE,L1,ITP1)
  592.             IF (IERR.NE.0) RETURN
  593.             IF (TYPRET.NE.'TABLE') GOTO 23
  594.             IF (ITP1.LE.0)         GOTO 23
  595. c           s'agit-il d'un corps rigide ?
  596.             TYPRET = ' '
  597.             CALL ACCTAB(ITP1,'MOT',I0,X0,'CORPS_RIGIDE',L0,IP0,
  598.      &                        TYPRET,I1,X1,MONMOT,L1,IP1)
  599.             IF (IERR.NE.0) RETURN
  600.             IF (TYPRET.EQ.'MOT') THEN
  601.               IF (MONMOT(1:4).EQ.'VRAI') THEN
  602.                 if (IDIM.eq.2 .and. IDIMB.lt.3) IDIMB = 3
  603.                 if (IDIM.eq.3 .and. IDIMB.lt.6) IDIMB = 6
  604.                 GOTO 23
  605.               ENDIF
  606.             ENDIF
  607. c           cas NNM : on veut recuperer la frequence du mode lineaire 
  608. c                     initial -> OMEG
  609.             IF(TYPS.EQ.'NNM ' .AND. IP.EQ.JMODE) THEN
  610.               CALL ACCTAB(ITP1,'MOT',I0,X0,'FREQUENCE',L0,IP0,
  611.      &                         'FLOTTANT',I1,X1,MONMOT,L1,IP1)
  612.               IF (IERR.NE.0) RETURN
  613.               OMEG = 2.D0*XPI * X1
  614.             ENDIF
  615.           GOTO 22
  616.  23       CONTINUE
  617. *      -Cas ou la base est un ensemble de bases
  618.         ELSE
  619.           IB = 0
  620.           NA2 = 0
  621. *         changement de dimension en cas de corps rigide
  622.           IR = 0
  623.  30       CONTINUE
  624.             IB = IB + 1
  625. c             write(ioimp,*) IB,'ieme base de l ensemble'
  626.             TYPRET = ' '
  627.             CALL ACCTAB(ITBAS2,'ENTIER',IB,X0,' ',L0,IP0,
  628.      &                         TYPRET,I1,X1,CHARRE,L1,ITBB)
  629.             IF (IERR.NE.0) RETURN
  630. c         --cas lecture table de la IB ieme base modale ok
  631.             IF (ITBB.NE.0) THEN
  632.             IF (TYPRET.EQ.'TABLE   ') THEN
  633.                CALL ACCTAB(ITBB,'MOT',I0,X0,'MODES',L0,IP0,
  634.      &                        'TABLE',I1,X1,' ',L1,IBAS)
  635.                IF (IERR.NE.0) RETURN
  636.                NNA2 = 0
  637.                IP = 0
  638.  32            CONTINUE
  639.                  IP = IP + 1
  640. c             write(ioimp,*) '  +',IP,'ieme mode de la ',IB,'ieme base'
  641.                  TYPRET = ' '
  642.                  CALL ACCTAB(IBAS,'ENTIER',IP,X0,' ',L0,IP0,
  643.      &                            TYPRET,I1,X1,CHARRE,L1,ITPP)
  644.                  IF (IERR.NE.0) RETURN
  645. c              --cas lecture table du IP ieme mode ok
  646.                  IF (ITPP.NE.0) THEN
  647.                  IF (TYPRET.EQ.'TABLE   ') THEN
  648. *                 changement de dimension en cas de corps rigide
  649.                   IF (IR.GT.1) GOTO 24
  650.                   TYPRET = ' '
  651.                   CALL ACCTAB(ITPP,'MOT',I0,X0,'CORPS_RIGIDE',L0,IP0,
  652.      &                    TYPRET,I1,X1,MONMOT,L1,IP1)
  653.                   IF (IERR.NE.0) RETURN
  654.                   IF (TYPRET.EQ.'MOT') THEN
  655.                     IF (MONMOT(1:4).EQ.'VRAI') THEN
  656.                       if (IDIM.eq.2 .and. IDIMB.lt.3) IDIMB = 3
  657.                       if (IDIM.eq.3 .and. IDIMB.lt.6) IDIMB = 6
  658.                     ENDIF
  659.                   ENDIF
  660.  24               CONTINUE
  661.                   NNA2 = NNA2 + 1
  662.                   GOTO 32
  663. c              --fin du cas lecture table du IP ieme mode ok
  664.                  ELSE
  665.                     CALL ERREUR(491)
  666.                     RETURN
  667.                  ENDIF
  668.                  ENDIF
  669. c              --fin du cas lecture table du IP ieme mode non ok
  670.                NA2 = MAX(NNA2,NA2)
  671.                GOTO 30
  672. c         --fin du cas lecture table de la IB ieme base modale ok
  673.             ELSE
  674.                CALL ERREUR(491)
  675.                RETURN
  676.             ENDIF
  677.             ENDIF
  678. c         --fin du cas lecture table de la IB ieme base modale non ok
  679.           NSB = IB - 1
  680.         ENDIF
  681. *      -fin distinction base modale simple / ensemble de bases
  682.         NPLSB = NPLB      
  683.       ENDIF
  684.       NPLSB=1
  685.       SEGINI,MTPHI
  686.       KTPHI = MTPHI
  687. *
  688. ************************************************************************
  689. *     Variables au cours d'un pas de temps:
  690. ************************************************************************
  691.       IF (IIMPI.EQ.333) write(IOIMP,*) 'HBMALO: SEGINI,MTPAS' 
  692. *      
  693.       MTPAS  = 0
  694. *       write(*,*) NA1,NPLB,IDIMB,NLIAA,NLIAB,NTVAR,IDIM,NPLB
  695.       SEGINI,MTPAS
  696.       KTPAS = MTPAS
  697.  
  698.  
  699. ************************************************************************
  700. *     Initialisation du segment representant les chargements ( bases A
  701. *     et B ):
  702. ************************************************************************
  703. *           
  704.       IF (ITCHAR.GT.0) THEN     
  705. c       NT1 deja dimensionne
  706.         SEGINI,MTFEX
  707.         KTFEX = MTFEX
  708.         TYPRET = ' '
  709.         CALL ACCTAB(ITCHAR,'MOT',I0,X0,'BALOURD',L0,IP0,
  710.      &                                      TYPRET,ICH,X1,CHARRE,L1,IP1)
  711.         IF (TYPRET.EQ.'ENTIER') THEN
  712.               BAL = ICH
  713.         ELSE
  714.               BAL = 0   
  715.         ENDIF      
  716. c       boucle sur les harmoniques
  717.         DO KHBM=-1*NHBM,NHBM
  718.             TYPRET = ' '
  719.             CALL ACCTAB(ITCHAR,'ENTIER',KHBM,X0,CHARR0,L0,IP0,
  720.      &                  TYPRET,I1,X1,CHARRE,L1,IP1)
  721. c            WRITE(IOIMP,*) 'HBMALO, Chargement. KHBM:',KHBM 
  722. c            WRITE(IOIMP,*) '-------------------------',TYPRET
  723.             IF (TYPRET.EQ.'CHARGEME') THEN
  724. c On ne veut pas d'objet de type %m1:8
  725.                 MOTERR(1:8)='CHARGEME'
  726.                 CALL ERREUR(39)
  727.                 RETURN
  728.             ELSEIF(TYPRET.EQ.'CHPOINT ') THEN
  729.               MCHPOI=IP1
  730.               SEGACT,MCHPOI 
  731.               NSOUPO=IPCHP(/1)
  732. c             boucle sur les zones (definies a partir des noms de composantes)
  733.               DO I=1,NSOUPO
  734.                  MSOUPO = IPCHP(I)
  735.                  SEGACT,MSOUPO
  736. c                on recupere le maillage
  737.                  MELEME = IGEOC
  738.                  SEGACT,MELEME
  739. c                NC = nombre de composantes
  740.                  NC = NOCOMP(/2)
  741.                  MPOVAL = IPOVAL
  742.                  SEGACT,MPOVAL
  743. c                VPOCHA(i,j) = valeur au noeud i de la composante j
  744.                  N = VPOCHA(/1)
  745.                  DO J=1,N
  746.                     DO K=1,NC
  747. c                      on cherche le #noeud dans IPOREF
  748.                        KNOE = NUM(1,J)
  749. c                      WRITE(IOIMP,*) 'KNOE=',KNOE
  750.                        CALL PLACE2(IPOREF,NPREF,IPOS,KNOE)
  751. *                     write(IOIMP,*) J,'eme noeud #',KNOE,' = mode',IPOS
  752. c                    IPOS est le numero de mode
  753. c                    WRITE(IOIMP,*) 'IPOS=',IPOS
  754. c                    WRITE(IOIMP,*) 'VPOCHA=',VPOCHA(J,K),'J=',J,',K=',K
  755.                        IF (IPOS.NE.0) THEN
  756.                           XFORCA = VPOCHA(J,K)
  757.                           IF (KHBM.le.0) then
  758.                             IND1=2*ABS(KHBM)
  759.                           else
  760.                             IND1=2*ABS(KHBM)-1
  761.                           endif
  762.                           IND=NA1*IND1+IPOS
  763. c                         WRITE(IOIMP,*) 'FEXA(',IND,')=',XFORCA
  764. c                         comme c'est constant on somme
  765.                           FEXA(IND) = FEXA(IND) + XFORCA
  766.                        ENDIF
  767.                     ENDDO
  768.                  ENDDO
  769.                  SEGDES,MPOVAL,MELEME,MSOUPO
  770.               ENDDO
  771.               SEGDES,MCHPOI
  772.             ENDIF
  773. c               On ne veut pas d'objet de type %m1:8
  774. c                MOTERR(1:8)=TYPRET
  775. c                CALL ERREUR(39)
  776. c                RETURN              
  777. c            ENDIF                
  778.         ENDDO
  779. *       mise a 0 de FEXPSM
  780.         do i1=1,NPLB
  781.         do i2=1,NPC1
  782.         do i3=1,2
  783.         do i4=1,IDIMB
  784.           FEXPSM(i1,i2,i3,i4)=0.D0
  785.         enddo
  786.         enddo
  787.         enddo
  788.         enddo
  789.       ELSE
  790.         SEGINI,MTFEX
  791.         KTFEX = MTFEX
  792.         DO I=1,Q1(/1)
  793.             FEXA(I) = 0.
  794.         ENDDO
  795.       ENDIF
  796. c     fin du cas ITCHAR existe ou pas
  797.  
  798.  
  799. **********************************************************************
  800. *     Segment des parametres numeriques:
  801. **********************************************************************
  802. c       SEGINI, PARNUM
  803. c       KPARNUM = PARNUM
  804. c deplace + haut
  805.       IF (IPARNUM.GT.0) THEN
  806.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'DS0',L0,IP0,
  807.      &                                  'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  808.           DS = X1
  809.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'DSMAX',L0,IP0,
  810.      &                               'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  811.           DSMAX=X1
  812.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'DSMIN',L0,IP0,
  813.      &                              'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  814.           DSMIN=X1
  815.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'ANGLE_MIN',L0,IP0,
  816.      &                              'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  817.           ANGMIN=X1
  818.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'ANGLE_MAX',L0,IP0,
  819.      &                               'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  820.           ANGMAX=X1
  821.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'ITERMOY',L0,IP0,
  822.      &                               'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  823.           ITERMOY=X1
  824.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'ITERMAX',L0,IP0,
  825.      &                               'ENTIER',I1,X1,CHARRE,L1,IP1)
  826.           ITERMAX=I1
  827.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'NBPAS',L0,IP0,
  828.      &                               'ENTIER',I1,X1,CHARRE,L1,IP1)
  829.           NBPAS=I1
  830.           NPAS =I1
  831.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'ISENS',L0,IP0,
  832.      &                               'FLOTTANT',I1,X1,CHARRE,L1,IP1)
  833.           ISENS=X1
  834.           TYPRET = ' '
  835.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'CALCUL_JAC',L0,IP0,
  836.      &                                 TYPRET,I1,X1,CHARRE,L1,IP1)
  837.           IF (TYPRET.EQ.'LOGIQUE') THEN 
  838.              JANAL = L1 
  839.           ELSE
  840.              JANAL = .FALSE.
  841.           ENDIF
  842.           TYPRET = ' '
  843.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'TOLERANCE',L0,IP0,
  844.      &                                 TYPRET,I1,X1,CHARRE,L1,IP1)
  845.           IF (TYPRET.EQ.'FLOTTANT') THEN 
  846.              TOLMIN = X1 
  847.           ELSE
  848.              TOLMIN = 1.E-8
  849.           ENDIF
  850.       ENDIF
  851. c 
  852. c *     Type = FORCe ou AUTOnome ou NonNormalMode?
  853. c       CALL ACCTAB(IPARNUM,'MOT',I0,X0,'TYPE',L0,IP0,
  854. c      &                                        'MOT',I1,X1,CHARRE,L1,IP1)    
  855. c       TYPS(1:4) = CHARRE(1:4)   
  856. c ==> deplace + haut
  857.  
  858.       IF (TYPS.EQ.'FORC') THEN   
  859. *     Le parametre de continuation sera la frequence de forcage
  860. *     (option par defaut pour l'instant).    
  861.           PARINI = OMEG
  862.           XPARA  = PARINI
  863.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'VAL_FIN',L0,IP0,
  864.      &                        'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  865.           PARFIN = X1           
  866.       ELSEIF (TYPS.EQ.'AUTO') THEN
  867. *         Continuation par rapport à un autre parametre.   
  868. *         On utilise comme frequence celle des conditions initiales.
  869. *         On a besoin des bornes initiales et finales
  870.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'VAL_INI',L0,IP0,
  871.      &                        'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  872.           PARINI = X1           
  873.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'VAL_FIN',L0,IP0,
  874.      &                        'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  875.           PARFIN = X1           
  876.       ELSEIF (TYPS.EQ.'NNM') THEN
  877. *         Calcul de mode non lineaire.
  878. *         La branche est suivie jusqu'a un certain niveau d'energie.
  879.           CALL ACCTAB(IPARNUM,'MOT',I0,X0,'VAL_FIN',L0,IP0,
  880.      &                        'FLOTTANT',I1,X1,CHARRE,L1,IP1)         
  881.           PARINI = XPARA
  882.           PARFIN = X1
  883.       ENDIF 
  884. *     test pas cher 
  885.       IF (abs(PARINI-PARFIN).le.(XZPREC**0.5)) THEN
  886. c         Les valeurs des parametres doivent etre distinctes
  887.           CALL ERREUR(402)
  888.           RETURN
  889.       ENDIF
  890.  
  891. c *     valeur initiale (selon sens de parcours)
  892. c       IF (ISENS.LT.0.) THEN
  893. c         XPARA = PARFIN
  894. c       ELSE
  895. c         XPARA = PARINI
  896. c       ENDIF
  897. cbp : ci-dessus desormais inutile car on utilise la table "INITIALE"
  898. cbp : PAR_INI et PAR_FIN ne sont plus necessairement ordonnes 
  899.  
  900.  
  901. **********************************************************************
  902. *     Tableau des resultats :
  903. **********************************************************************
  904.        NBIFU=MIN(NPAS/2+1,100)
  905.        SEGINI, PSORT
  906.        KSORT= PSORT
  907.  
  908.  
  909. ************************************************************************
  910. *    Impressions :
  911. ************************************************************************
  912.       IF (IIMPI.GE.333) THEN
  913.  
  914.          WRITE(IOIMP,*) 'Frequence initiale: W0 =',OMEG
  915.          WRITE(IOIMP,*) 'Le systeme est de type: ',TYPS
  916. *        DO i=1,Q1(/1)
  917. *          write(*,*) 'Q1 initial=',(Q1(i))
  918. *        ENDDO
  919.  
  920.          WRITE(IOIMP,*)'   segment MTLIAB'
  921.          WRITE(IOIMP,*)'     NLIAB  =',IPALB(/1)
  922.          WRITE(IOIMP,*)'     NIPALB =',IPALB(/2)
  923.          WRITE(IOIMP,*)'     NXPALB =',XPALB(/2)
  924.          WRITE(IOIMP,*)'     NPLBB  =',IPLIB(/2)
  925.          WRITE(IOIMP,*)'     NPLB   =',JPLIB(/1)
  926.          WRITE(IOIMP,*)'     NIP    =',XABSCI(/2)
  927. *
  928.          WRITE(IOIMP,*)'   segment MTLIAA'
  929.          WRITE(IOIMP,*)'     NLIAA  =',IPALA(/1)
  930.          WRITE(IOIMP,*)'     NIPALA =',IPALA(/2)
  931.          WRITE(IOIMP,*)'     NXPALA =',XPALA(/2)
  932.          WRITE(IOIMP,*)'     NPLAA  =',IPLIA(/2)
  933.          WRITE(IOIMP,*)'     NPLA   =',JPLIA(/1)
  934. *
  935.          WRITE(IOIMP,*)'   segment MTFEX : chargement frequentiel'
  936.          IF (BAL.EQ.1)WRITE(IOIMP,*)'de type balourd'
  937.          DO i=1,FEXA(/1)
  938.             write(IOIMP,*) 'FEXA(',i,',:)=',(FEXA(i))
  939.          ENDDO
  940.       ENDIF
  941. *
  942.       END
  943.  
  944.  
  945.  
  946.  
  947.  
  948.  
  949.  
  950.  
  951.  
  952.  

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