Télécharger la plaquette Cast3M


fr en

Télécharger hbmfn2.eso

Retour à la liste

Numérotation des lignes : 

hbmfn2

  1. C HBMFN2    SOURCE    OF166741  26/05/11    21:15:10     12538          
  2.  
  3. *=======================================================================
  4. *      Calcul de FNL(X) et dFNLdX par AFT (Alternating Frequency Time)
  5. *=======================================================================
  6.  
  7.       SUBROUTINE HBMFN2(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,KTPHI,KTLIAA,
  8.      &                  KTLIAB,KTPAS,KTFEX,KOCLFA,KOCLB1,MFNL)
  9.  
  10. *----- Declarations ----------------------------------------------------
  11.  
  12.         IMPLICIT INTEGER(I-N)
  13.         IMPLICIT REAL*8(A-H,O-Z)
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17.  
  18. -INC SMLREEL
  19.       POINTEUR MFNL.MLREEL
  20.  
  21. -INC TMDYNC
  22.       POINTEUR JACV.MATWRK
  23.  
  24.       SEGMENT mwork
  25.         REAL*8 XT(NDDL,NFFT),FT(NDDL,NFFT),VT(NDDL,NFFT)
  26.         REAL*8 XAUX(NDDL,4),VAUX(NDDL,4),VCX(NFFT),VCV(NFFT)
  27.         REAL*8 MATJX(NDDL,NDDL,NFFT), MATJV(NDDL,NDDL,NFFT)
  28.         REAL*8 auxDL(2*NHBM+1,NFFT), JFR(2*NHBM+1,2*NHBM+1)
  29.         REAL*8 auxVDL(2*NHBM+1,NFFT), JVFR(2*NHBM+1,2*NHBM+1)
  30.         REAL*8 FNL2(NT),V1(NT)
  31.       ENDSEGMENT
  32.  
  33. c     segment to call fftpack5.1
  34.       SEGMENT wfft51
  35.         real*8 work(lenwrk)
  36.         real*8 wwsave(lensav)
  37.         real*8 XX51(NCOU)
  38.         real*8 XV51(NCOU)
  39.       ENDSEGMENT
  40.  
  41.  
  42.       INTEGER NP,NDDL,NT,NFFT
  43.       REAL*8 PDT
  44.       LOGICAL RIGIDE
  45.  
  46.       REAL*8 ZERO,ONE,TWO
  47.       PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0)
  48.  
  49. *   Variables generalisees: coefficients de Fourier et frequence
  50.       MTQ=KTQ
  51. *   Matrices XK,KASM,XM,GAM,IGAM,DL
  52.       MTKAM=KTKAM
  53. *   Deformees modales
  54.       MTPHI = KTPHI
  55.       NSB   = XPHILB(/1)
  56.       NPLSB = XPHILB(/2)
  57.       NA2   = XPHILB(/3)
  58.       IDIMB = XPHILB(/4)
  59. *   Liaisons sur base A
  60.       MTLIAA = KTLIAA
  61.       NLIAA = IPALA(/1)
  62. *   Liaisons sur base B
  63.       MTLIAB = KTLIAB
  64.       NLIAB  = IPALB(/1)
  65.       NIP    = XABSCI(/2)
  66.       NPLB   = JPLIB(/1)
  67. *   Forces externes
  68.       MTFEX=KTFEX
  69. *   Truc local base A
  70.       LOCLFA = KOCLFA
  71. *   Truc local base B
  72.       LOCLB1 = KOCLB1
  73. *   Inconnues sur un pas de temps (AFT)
  74.       MTPAS=KTPAS
  75.       RIGIDE =.FALSE.
  76.  
  77.       SEGINI,MWORK
  78.       nmc = NT
  79.       SEGINI,JACV
  80.  
  81. c   Initialisation FFT
  82.       IERRD=0
  83.       NCOU=NFFT
  84.       lenwrk = 2*NFFT
  85.       lensav = NFFT + INT(LOG(REAL(NFFT))/LOG(TWO)) + 4
  86.       SEGINI,wfft51
  87.  
  88. *------ Time-domain displacements and velocities: XT, VT ---------------
  89.  
  90.       CALL HBMDVEC(NT,NHBM,NDDL,Q1,OMEG,V1)
  91.       CALL IFT(Q1,NT,NDDL,NHBM,NFFT,GAM,XT)
  92.       CALL IFT(V1,NT,NDDL,NHBM,NFFT,GAM,VT)
  93.  
  94. *------ Time-domain nonlinear forces, tangent matrix: FT, KTOTBA -------
  95.  
  96. c   >>> Loop over time steps >>>
  97.         PDT = ONE/NFFT
  98.         DO I = 1,NFFT
  99.  
  100. c          initialisation des vecteurs du pas de temps
  101.            DO L=1,NDDL
  102.              XAUX(L,1) = XT(L,I)
  103.              VAUX(L,1) = VT(L,I)
  104.              FTOTA(L,1)=ZERO
  105.            ENDDO
  106.            DO JJ=1,NDDL
  107.              DO IJ=1,NDDL
  108.                KTOTXA(IJ,JJ) = ZERO
  109.                KTOTVA(IJ,JJ) = ZERO
  110.              ENDDO
  111.            ENDDO
  112.  
  113. c        Fnl(x,v) en base A (ddls modaux)
  114.            IF (NLIAA.NE.0) THEN
  115.          CALL D2VLFA(XAUX,VAUX,FTOTA,NDDL,IPALA,IPLIA,XPALA,XVALA,
  116.      &               NLIAA,PDT,ZERO,I,1,FINERT,1,LOCLFA.FTESTA,
  117.      &               KTOTXA,KTOTVA,.TRUE.)
  118.            ENDIF
  119.  
  120. c        Fnl(x,-) en base B (ddls physique)
  121.            IF (NLIAB.NE.0) THEN
  122.          CALL D2VLFB(XAUX,VAUX,FTOTA,NDDL,IPALB,IPLIB,XPALB,XVALB,NLIAB,
  123.      & XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,PDT,ZERO,I,IBASB,IPLSB,
  124.      &   INMSB,IORSB,NSB,NPLSB,NA2,1,FEXPSM,NPC1,IERRD,LOCLB1.FTESTB,
  125.      &   XABSCI,XORDON,NIP,FEXB,RIGIDE,IAROTA,XCHPFB,
  126.      &   KTOTXA,KTOTVA,KTOTXB,KTOTVB,.TRUE.)
  127.            ENDIF
  128.  
  129. c          stockage de la force + JAC au pas de temps I
  130.            DO J=1,NDDL
  131.                FT(J,I) = FTOTA(J,1)
  132.                DO K = 1,NDDL
  133.                    MATJX(J,K,I) = KTOTXA(J,K)
  134.                    MATJV(J,K,I) = KTOTVA(J,K)
  135.                ENDDO
  136.            ENDDO
  137.  
  138.         ENDDO
  139. c   <<< end of Loop over time steps <<<
  140.  
  141.  
  142. c------ Fourier coefficients  ------------------
  143. c       Initialize the WWSAVE array.
  144.         CALL rfft1i (NFFT,wwsave,lensav,ier)
  145.         IF (IERR.ne.0) RETURN
  146.  
  147. c------ Fourier coefficients of nonlinear forces: FNL ------------------
  148. *        CALL DFT(FT,NDDL,NHBM,NFFT,IGAM,MFNL.PROG)
  149.         DO IM = 1,NDDL
  150.           DO KI = 1,NFFT
  151.             XX51(KI) = FT(IM,KI)
  152.           ENDDO
  153.           CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  154.           IA = 0
  155.           DO KI = IM,NT-(NDDL-IM),NDDL
  156.            IA = IA+1
  157.            MFNL.PROG(KI) = XX51(IA)
  158.           ENDDO
  159.         ENDDO
  160.  
  161. c------ Fourier coefficients of nonlinear tangent matrix: dFNLdX -------
  162.  
  163. c   >>> Loop over modes I,J >>>
  164.         DO I = 1,NDDL
  165.           DO J = 1,NDDL
  166.  
  167. c           Extract row vectors for the (i,j)-th derivative term
  168.             DO IK = 1,NFFT
  169.               VCX(IK) = MATJX(I,J,IK)
  170.               VCV(IK) = MATJV(I,J,IK)
  171.             ENDDO
  172. c           First series of FFTs on the columns of diagonal dfdx matrix
  173.             DO II = 1,NFFT
  174.               DO KI = 1,NFFT
  175.                 XX51(KI) = ZERO
  176.                 XV51(KI) = ZERO
  177.               ENDDO
  178.               XX51(II) = VCX(II)
  179.               XV51(II) = VCV(II)
  180.              CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  181.              CALL rfft1f(NFFT,1,XV51,NFFT,wwsave,lensav,work,lenwrk,ier)
  182.               DO JK = 1,2*NHBM+1
  183.                 auxDL(JK,II)  = NFFT*XX51(JK)
  184.                 auxVDL(JK,II) = NFFT*XV51(JK)
  185.               ENDDO
  186.             ENDDO
  187. c           Second series of FFTs on the rows of auxDL
  188.             DO II = 1,2*NHBM+1
  189.               DO KI = 1,NFFT
  190.                   XX51(KI) = auxDL(II,KI)
  191.                   XV51(KI) = auxVDL(II,KI)
  192.               ENDDO
  193.              CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  194.              CALL rfft1f(NFFT,1,XV51,NFFT,wwsave,lensav,work,lenwrk,ier)
  195.                  JFR(II,1)  = XX51(1)
  196.                 JVFR(II,1)  = XV51(1)
  197.               DO JK = 2,2*NHBM+1
  198.                 JFR(II,JK)  = XX51(JK)*0.5D0
  199.                 JVFR(II,JK) = XV51(JK)*0.5D0
  200.               ENDDO
  201.             ENDDO
  202.  
  203. c           Assemble the total matrices JAC, JACV
  204.             IL = 0
  205.             DO NI = I,NDDL*(2*NHBM)+I,NDDL
  206.               IL = IL + 1
  207.               IC = 0
  208.               DO NJ = J,NDDL*(2*NHBM)+J,NDDL
  209.                 IC = IC + 1
  210.                 JAC(NI,NJ) = JFR(IL,IC)
  211.                 JACV.MATRC(NI,NJ)= JVFR(IL,IC)
  212.               ENDDO
  213.             ENDDO
  214.  
  215.           ENDDO
  216.         ENDDO
  217. c   <<< fin de la boucle sur les modes I,J <<<
  218.  
  219. *       Derivative of the velocity term
  220.         CALL HBMDMAT(NT,NDDL,OMEG,JACV)
  221.  
  222. *       Total nonlinear tangent stiffness: dFNLdX = JAC + w*JACV*Nabla
  223.         DO I = 1,NT
  224.           DO J = 1,NT
  225.               JAC(I,J) = JAC(I,J) + JACV.MATRC(I,J)
  226.           ENDDO
  227.         ENDDO
  228.  
  229. c------ Menage ------------------
  230.       SEGSUP,MWORK,JACV
  231.  
  232.       RETURN
  233.       END
  234.  
  235.  
  236.  

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