Télécharger la plaquette Cast3M


fr en

Télécharger hbmtra.eso

Retour à la liste

Numérotation des lignes : 

hbmtra

  1. C HBMTRA    SOURCE    CB215821  23/01/25    21:15:23     11573          
  2.  
  3. C DEVTRA    SOURCE    BP208322  15/10/16    21:15:01     8685
  4.       SUBROUTINE HBMTRA(ITBAS,ITKM,ITA,KTKAM,IPMAIL,NHBM,KTRES,KTNUM,
  5.      &                  KPREF,KTPHI,KTLIAB,RIGIDE)
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8(A-H,O-Z)
  8. *--------------------------------------------------------------------*
  9. *                                                                    *
  10. *     Operateur DYNC                                                 *
  11. *     ________________________________________________               *
  12. *                                                                    *
  13. *     Transpose l'information des objets de Castem2000 dans des      *
  14. *     tableaux de travail.                                           *
  15. *                                                                    *
  16. *     Parametres:                                                    *
  17. *                                                                    *
  18. * e   ITBAS   Table representant une base modale                     *
  19. * e   ITKM    Table contenant les matrices XK et XM                  *
  20. * e   ITA     Table contenant la matrice XASM                        *
  21. * es  KTKAM   Segment contenant les matrices XK, XASM et XM          *
  22. * e   IPMAIL  Maillage de reference pour les CHPOINTs resultats      *
  23. * es  KTRES   Segment de sauvegarde des resultats                    *
  24. * e   KPREF   Segment des points de reference                        *
  25. * es  KTPHI   Segment des deformees modales                          *
  26. * e   KTLIAB  Segment des liaisons sur base B                        *
  27. * e   RIGIDE  Vrai si corps rigide, faux sinon                       *
  28. *                                                                    *
  29. *--------------------------------------------------------------------*
  30. -INC SMRIGID
  31. -INC SMELEME
  32. -INC PPARAM
  33. -INC CCOPTIO
  34. -INC CCREEL
  35. *
  36. *-INC TMDYNC.INC
  37. ************************** debut TMDYNC.INC ****************************
  38.  
  39. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  40. *    TODO   : a extraire dans un include des que stabilise
  41. *
  42. *     Segment des variables generalisees:
  43. *     -----------------------------------
  44.       SEGMENT MTQ
  45.          REAL*8 Q1(NT1)
  46.          REAL*8 OMEG,XPARA
  47.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  48.          REAL*8 dX(NT1), dw, dv
  49.       ENDSEGMENT
  50. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  51. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  52. *            avec q_i vecteur de dimension n ou n=nombre de modes
  53. *     OMEG : frequence fondamentale de l'approximation
  54. *     XPARA: parametre de continuation (par defaut la frequence) 
  55. *            \in [PARINI,PARFIN]
  56. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  57. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  58. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  59. *            lineaires et constantes
  60. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  61. *
  62. *
  63. *     Segment contenant les matrices XK, XASM et XM:
  64. *     ---------------------------------------------
  65.       SEGMENT MTKAM
  66.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  67.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  68. *         REAL*8 GAMFIN(NPC2,nl1)
  69.       ENDSEGMENT
  70. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  71. *     GAM et IGAM : matrices pour la FFT et son inverse
  72. *     GAMFIN      : 
  73. *
  74. *     Segment des deformees modales:
  75. *     ------------------------------
  76. *     (idem DYNE)
  77.       SEGMENT MTPHI
  78.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  79.          INTEGER IAROTA(NSB)
  80.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  81.       ENDSEGMENT
  82. *
  83. *     Segment descriptif des liaisons en base A:
  84. *     ------------------------------------------
  85. *     (idem DYNE)
  86.       SEGMENT MTLIAA
  87.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  88.          REAL*8  XPALA(NLIAA,NXPALA)
  89.       ENDSEGMENT
  90. *
  91. *     Segment descriptif des liaisons en base B:
  92. *     ------------------------------------------
  93. *     (idem DYNE)
  94.       SEGMENT MTLIAB
  95.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  96.          REAL*8  XPALB(NLIAB,NXPALB)
  97.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  98.       ENDSEGMENT
  99. *
  100. *     Segment representant les chargements exterieurs:
  101. *     -----------------------------------------------
  102.       SEGMENT MTFEX
  103.          REAL*8  FEXA(NT1)
  104.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  105.          INTEGER BAL
  106.       ENDSEGMENT
  107. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  108. *            Fourier et exprimes en base A
  109. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  110. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  111. *             compatibilite avec calcul des Fnl.
  112. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  113. *            (cad proportionnel a OMEG**2)
  114. *
  115. *     Segment "local" pour DEVLFA:
  116. *     ----------------------------
  117.       SEGMENT LOCLFA
  118.          REAL*8 FTEST(NA1,4)
  119.       ENDSEGMENT
  120. *      
  121. *     Segment "local" pour DEVLB1:
  122. *     ----------------------------
  123.       SEGMENT LOCLB1
  124.          REAL*8 FTEST2(NPLB,6)
  125.       ENDSEGMENT
  126. *      
  127. *     Segment contenant les variables au cours d un pas de temps:
  128. *     ----------------------------------------------------------
  129.       SEGMENT MTPAS
  130.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  131.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  132.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  133.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  134.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  135.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  136.       ENDSEGMENT
  137. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  138. *     XPTB       : deplacement du point d'une liaison en base B
  139. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  140. *     FEXB       : forces exterieures en base B (a priori uniquement 
  141. *                  pour les moments appliques aux rotations rigides ?)
  142. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  143. *                  maillage de contact dans certaines liaisons)
  144. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  145. *                       base A/B (= contributions a dFnl/dX)
  146. *
  147. *
  148. *     Segment des points de reference des modes (base A):
  149. *     --------------------------------------------------
  150.       SEGMENT MPREF
  151.          INTEGER IPOREF(NPREF)
  152.       ENDSEGMENT
  153. *
  154. *     Segment des points en base B:
  155. *     -----------------------------
  156.       SEGMENT NCPR(NBPTS)
  157. *     NCRP(#global) = #local dans XPTB (1er indice)      
  158. *
  159. *     Segment des parametres numeriques pour la continuation: 
  160. *     ------------------------------------------------------
  161.       SEGMENT PARNUM
  162.          CHARACTER*4 TYPS
  163.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  164.          REAL*8   PARINI,PARFIN
  165.          INTEGER  ITERMAX,NBPAS
  166.          LOGICAL  JANAL
  167.       ENDSEGMENT
  168. *
  169. *     Segment des resultats:
  170. *     ---------------------
  171.       SEGMENT  PSORT
  172.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  173.         REAL*8 VSAVE(NPAS)
  174.         LOGICAL ZSAVE(NPAS)
  175.         CHARACTER*2 TYPBIF(NBIFU)
  176.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  177.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  178.         INTEGER CBIF
  179.       ENDSEGMENT
  180. *     QSAVE(i,j) = Q harmonique i au pas j
  181. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  182. *     ZSAVE(j)   = stabilite au j-eme pas
  183. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  184. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  185. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  186. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  187. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  188. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  189.  
  190. *     Segment des tableaux de travail:
  191. *     -------------------------------
  192.       SEGMENT  MTEMP
  193.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  194.          REAL*8 T02(NT1+2), TP2(NT1+2)
  195.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  196.          REAL*8 res
  197.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  198.          REAL*8 QOLD(NT1),OMEGOLD
  199.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  200.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  201.       ENDSEGMENT
  202. *     Jacobiennes augmentees 
  203. *     Ja : [ RX  Rw ; dX dw]
  204. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  205.  
  206. *      SEGMENT  NNNN
  207. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  208. *      ENDSEGMENT
  209.  
  210. *************************** fin TMDYNC.INC *****************************
  211.  
  212.  
  213.       LOGICAL L0,L1,RIGIDE
  214.       CHARACTER*4 CMOT,MOINC
  215.       CHARACTER*8 TYPRET,CHARRE
  216.       CHARACTER*40 MONMOT
  217. *
  218.       MTKAM  = KTKAM
  219.       MTPHI  = KTPHI
  220.       MTLIAB = KTLIAB
  221.       MPREF  = KPREF
  222.       MTNUM  = KTNUM
  223.  
  224. *     dimensions de MTPHI
  225.       NPLB  = IBASB(/1)
  226.       NSB   = INMSB(/1)
  227.       NA2   = XPHILB(/3)
  228.       IDIMB = XPHILB(/4)
  229.       NLIAB = IPALB(/1)
  230.  
  231. *     dimensions de MTKAM
  232.       NA1  = XASM(/1)
  233.       NB1K = XK(/2)
  234.       NB1C = XASM(/2)
  235.       NB1M = XM(/2)
  236.  
  237. *     dimensions de MTQ     
  238. *      NT1 = NA1*(2*NHBM+1)
  239.  
  240.       NPREF=IPOREF(/1)
  241. *
  242.       IA1 = 0
  243.       DEUXPI = 2.D0 * XPI
  244.       RIGIDE =.FALSE.
  245. *
  246. *     Traitement des matrices de variables generalisees:
  247. *
  248.       IF (ITBAS.NE.0 .AND.ITKM.EQ.0) THEN
  249.          IF (IIMPI.EQ.333)
  250.      &   WRITE(IOIMP,*) 'HBMTRA: cas table BASE_DE_MODES, quel type?'
  251.          CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'SOUSTYPE',L0,IP0,
  252.      &                     'MOT',I1,X1,MONMOT,L1,IP1)
  253.          IF (IERR.NE.0) RETURN
  254. *
  255. *        Cas ou la base est unique
  256. *
  257.          IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  258.             IF (IIMPI.EQ.333)
  259.      &      WRITE(IOIMP,*) 'HBMTRA: lecture table BASE_MODALE'
  260. *
  261. *           On recupere la base de modes
  262. *
  263.             CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  264.      &                        'TABLE',I1,X1,' ',L1,IBAS)
  265.             IF (IERR.NE.0) RETURN
  266.             CALL HBM26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,RIGIDE,ITKM)
  267. *            CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,RIGIDE,
  268. *     &                  0,.false.,ITKM)
  269.             IF (RIGIDE) THEN
  270.                RIGIDE =.FALSE.
  271.                DO 80 ILIA =1,NLIAB
  272.                    ITYP = IPALB(ILIA,1)
  273.                    IF (ITYP.EQ.35.OR.ITYP.EQ.36) THEN
  274.                        RIGIDE =.TRUE.
  275.                    ENDIF
  276. 80             CONTINUE
  277.             ENDIF
  278.             IF (IERR.NE.0) RETURN
  279. *
  280. *        Cas ou on a un ensemble de bases
  281. *
  282.          ELSE IF (MONMOT(1:17).EQ.'ENSEMBLE_DE_BASES') THEN
  283.             IF (IIMPI.EQ.333)
  284.      &      WRITE(IOIMP,*) 'HBMTRA: lecture table ENSEMBLE_DE_BASES'
  285. *
  286. *           On boucle sur le nombre de bases
  287. *
  288.             IT = 0
  289.             NPLSB = 0
  290.  10         CONTINUE
  291.             TYPRET = ' '
  292.             IT = IT + 1
  293.             CALL ACCTAB(ITBAS,'ENTIER',IT,X0,' ',L0,IP0,
  294.      &                          TYPRET,I1,X1,CHARRE,L1,ITTBAS)
  295.             IF (IERR.NE.0) RETURN
  296.             IF (ITTBAS.NE.0) THEN
  297.                IF (TYPRET.EQ.'TABLE   ') THEN
  298.                   CALL ACCTAB(ITTBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  299.      &                             'TABLE',I1,X1,' ',L1,IBAS)
  300.                   IF (IERR.NE.0) RETURN
  301.                   CALL HBM26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  302.      &                              RIGIDE,ITKM)           
  303. *                  CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  304. *     &                              RIGIDE,0,.false.,ITKM)
  305.                   IF (IERR.NE.0) RETURN
  306.                   NPLSB = MAX(NPLSB,ICOMP)
  307.                   GOTO 10
  308.                ELSE
  309.                   CALL ERREUR(491)
  310.                   RETURN
  311.                ENDIF
  312.             ENDIF
  313.          ENDIF
  314. *
  315.       ELSE IF (ITKM.NE.0) THEN
  316. *        cas table RAIDEUR_ET_MASSE non prevu pour l'instant
  317.          CALL ERREUR(491)
  318.          RETURN
  319.       ENDIF
  320. *
  321. *     Traitement de la matrice d'amortissement
  322. *
  323.       IF (ITA.NE.0) THEN
  324.          IF (IIMPI.EQ.333)
  325.      &   WRITE(IOIMP,*) 'HBMTRA: cas table AMORTISSEMENT'
  326.            TYPRET = ' '
  327.            CALL ACCTAB(ITA,'MOT',I0,X0,'AMORTISSEMENT',L0,IP0,
  328.      &                    TYPRET,I1,X1,CHARRE,L1,IAMOR)
  329.            IF (IERR.NE.0) RETURN
  330.          IF (IAMOR.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  331.             IF (IIMPI.EQ.333)
  332.      &      WRITE(IOIMP,*) 'HBMTRA: lecture table AMORTISSEMENT ok'
  333.             MRIGID = IAMOR
  334.             SEGACT,MRIGID
  335.             NAMOR = IRIGEL(/2)
  336.             DO 60 I=1,NAMOR
  337.                COEF = COERIG(I)
  338. c          write(ioimp,*) 'HBMTRA: sous rigidite ',I,'/',NAMOR,COEF
  339.                MELEME = IRIGEL(1,I)
  340.                DESCR  = IRIGEL(3,I)
  341.                XMATRI = IRIGEL(4,I)
  342.                SEGACT,DESCR,MELEME,XMATRI
  343.                NRIG = RE(/3)
  344.                LVAL = RE(/1)
  345.                DO 70 IRIG=1,NRIG
  346. c          write(ioimp,*) 'HBMTRA: + element',IRIG,'/',NRIG
  347. c                 boucle sur les lignes (ddls duals)
  348.                   DO 75 IN=1,LVAL
  349.                     INODE=NOELED(IN)
  350.                     IF(INODE.ne.NOELEP(IN)) THEN
  351.                        WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  352.      &                 'primales et duales de la matrice AMORTISSEMENT'
  353.                        CALL ERREUR(47)
  354.                        RETURN
  355.                     ENDIF
  356.                     NNODE=NUM(INODE,IRIG)
  357. c          write(ioimp,*) 'HBMTRA:   + noeud dual',IN,'/',LVAL,' #',NNODE
  358. c                   position de cette inconnue dans IPOREF de MPREF
  359.                     DO 76 IA=1,NPREF
  360.                       IF (IPOREF(IA).EQ.NNODE) GOTO 79
  361.  76                 CONTINUE
  362.                     write(ioimp,*) 'HBMTRA: Incoherence entre les ',
  363.      &              'points de reference et la matrice AMORTISSEMENT'
  364.                     CALL ERREUR(504)
  365.  79                 CONTINUE
  366. c          write(ioimp,*) 'HBMTRA:   + noeud dual trouvé en position',IA
  367. *                     Partie diagonale seulement ...
  368.                       XASM(IA,1) = XASM(IA,1) + (RE(IN,IN,IRIG) * COEF)
  369.  75               CONTINUE
  370.  70            CONTINUE
  371.                SEGDES,XMATRI,MELEME,DESCR
  372.  60         CONTINUE
  373.             SEGDES,MRIGID
  374.          ELSE
  375.             CALL ERREUR(485)
  376.             RETURN
  377.          ENDIF
  378.       ENDIF
  379. *
  380.       IF (IIMPI.EQ.333) THEN
  381.          WRITE(IOIMP,*)'     segment MTPHI'
  382.          WRITE(IOIMP,*)'HBMTRA : valeur de NPLB  :',IBASB(/1)
  383.          WRITE(IOIMP,*)'HBMTRA : valeur de NSB   :',XPHILB(/1)
  384.          WRITE(IOIMP,*)'HBMTRA : valeur de NPLSB :',XPHILB(/2)
  385.          WRITE(IOIMP,*)'HBMTRA : valeur de NA2   :',XPHILB(/3)
  386.          WRITE(IOIMP,*)'HBMTRA : valeur de IDIMB :',XPHILB(/4)
  387.          WRITE(IOIMP,*)'     segment MTKAM'
  388.          WRITE(IOIMP,*)'NA1,NB1K,NB1C,NB1M=',NA1,NB1K,NB1C,NB1M
  389.          if(NB1K.gt.1) then
  390.             do iou=1,NA1
  391.                WRITE(IOIMP,*) 'XK=',(XK(iou,jou),jou=1,NB1K)
  392.             enddo
  393.          else
  394.             do iou=1,NA1
  395.                WRITE(IOIMP,*) 'XK(',iou,',1)=',XK(iou,1)
  396.             enddo
  397.          endif
  398.          if(NB1C.gt.1) then
  399.             do iou=1,NA1
  400.               WRITE(IOIMP,*) 'XASM=',(XASM(iou,jou),jou=1,NB1C)
  401.             enddo
  402.          else
  403.             do iou=1,NA1
  404.               WRITE(IOIMP,*) 'XASM(',iou,',1)=',XASM(iou,1)
  405.             enddo
  406.          endif
  407.          if(NB1M.gt.1) then
  408.             do iou=1,NA1
  409.                WRITE(IOIMP,*) 'XM=',(XM(iou,jou),jou=1,NB1M)
  410.             enddo
  411.          else
  412.             do iou=1,NA1
  413.                WRITE(IOIMP,*) 'XM(',iou,',1)=',XM(iou,1)
  414.             enddo
  415.          endif
  416.        ENDIF
  417. *
  418.       END
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  

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