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

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