Télécharger la plaquette Cast3M


fr en

Télécharger funcv.eso

Retour à la liste

Numérotation des lignes : 

funcv

  1. C FUNCV     SOURCE    FANDEUR   22/05/02    21:15:19     11359          
  2.  
  3.         SUBROUTINE funcv(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,KTPHI,KTLIAA,
  4.      &             KTLIAB,KTPAS,KTFEX,KOCLFA,KOCLB1,JANA2,fvec,fref,NNM)
  5. *
  6. *=======================================================================  
  7. *       Calcul du residu R(X) = ZX-Fnl-Fext -> fvec 
  8. *       et de la Jacobienne dFnl/dX         -> JAC(*,*) du segment KTQ
  9. *=======================================================================
  10.  
  11. *----- Declarations ----------------------------------------------------
  12.  
  13.       IMPLICIT INTEGER(I-N)
  14.       IMPLICIT REAL*8(A-H,O-Z)
  15.  
  16.       INTEGER NDDL,NFFT,NT,NHBM
  17.       SEGMENT mwork
  18.       REAL*8 ZX(NT),VCTCS(7),LAMBD(NDDL),AiDi(2),Di(2),bi(2)
  19.       ENDSEGMENT
  20.       REAL*8 FNLP(NT),FNL(NT),fvec(NT),PP(NT)
  21.       REAL*8 EPS, summ, Q1SAVE 
  22. c     NNM: cas des modes non-lineaires  
  23.       LOGICAL JANA2,NNM        
  24.       REAL*8 fref
  25. *       
  26. *-INC TMDYNC.INC
  27. ************************** debut TMDYNC.INC ****************************
  28.  
  29. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  30. *    TODO   : a extraire dans un include des que stabilise
  31. *
  32. *     Segment des variables generalisees:
  33. *     -----------------------------------
  34.       SEGMENT MTQ
  35.          REAL*8 Q1(NT1)
  36.          REAL*8 OMEG,XPARA
  37.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  38.          REAL*8 dX(NT1), dw, dv
  39.       ENDSEGMENT
  40. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  41. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  42. *            avec q_i vecteur de dimension n ou n=nombre de modes
  43. *     OMEG : frequence fondamentale de l'approximation
  44. *     XPARA: parametre de continuation (par defaut la frequence) 
  45. *            \in [PARINI,PARFIN]
  46. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  47. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  48. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  49. *            lineaires et constantes
  50. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  51. *
  52. *
  53. *     Segment contenant les matrices XK, XASM et XM:
  54. *     ---------------------------------------------
  55.       SEGMENT MTKAM
  56.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  57.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  58. *         REAL*8 GAMFIN(NPC2,nl1)
  59.       ENDSEGMENT
  60. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  61. *     GAM et IGAM : matrices pour la FFT et son inverse
  62. *     GAMFIN      : 
  63. *
  64. *     Segment des deformees modales:
  65. *     ------------------------------
  66. *     (idem DYNE)
  67.       SEGMENT MTPHI
  68.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  69.          INTEGER IAROTA(NSB)
  70.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  71.       ENDSEGMENT
  72. *
  73. *     Segment descriptif des liaisons en base A:
  74. *     ------------------------------------------
  75. *     (idem DYNE)
  76.       SEGMENT MTLIAA
  77.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  78.          REAL*8  XPALA(NLIAA,NXPALA)
  79.       ENDSEGMENT
  80. *
  81. *     Segment descriptif des liaisons en base B:
  82. *     ------------------------------------------
  83. *     (idem DYNE)
  84.       SEGMENT MTLIAB
  85.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  86.          REAL*8  XPALB(NLIAB,NXPALB)
  87.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  88.       ENDSEGMENT
  89. *
  90. *     Segment representant les chargements exterieurs:
  91. *     -----------------------------------------------
  92.       SEGMENT MTFEX
  93.          REAL*8  FEXA(NT1)
  94.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  95.          INTEGER BAL
  96.       ENDSEGMENT
  97. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  98. *            Fourier et exprimes en base A
  99. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  100. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  101. *             compatibilite avec calcul des Fnl.
  102. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  103. *            (cad proportionnel a OMEG**2)
  104. *
  105. *     Segment "local" pour DEVLFA:
  106. *     ----------------------------
  107.       SEGMENT LOCLFA
  108.          REAL*8 FTEST(NA1,4)
  109.       ENDSEGMENT
  110. *      
  111. *     Segment "local" pour DEVLB1:
  112. *     ----------------------------
  113.       SEGMENT LOCLB1
  114.          REAL*8 FTEST2(NPLB,6)
  115.       ENDSEGMENT
  116. *      
  117. *     Segment contenant les variables au cours d un pas de temps:
  118. *     ----------------------------------------------------------
  119.       SEGMENT MTPAS
  120.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  121.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  122.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  123.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  124.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  125.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  126.       ENDSEGMENT
  127. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  128. *     XPTB       : deplacement du point d'une liaison en base B
  129. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  130. *     FEXB       : forces exterieures en base B (a priori uniquement 
  131. *                  pour les moments appliques aux rotations rigides ?)
  132. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  133. *                  maillage de contact dans certaines liaisons)
  134. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  135. *                       base A/B (= contributions a dFnl/dX)
  136. *
  137. *
  138. *     Segment des points de reference des modes (base A):
  139. *     --------------------------------------------------
  140.       SEGMENT MPREF
  141.          INTEGER IPOREF(NPREF)
  142.       ENDSEGMENT
  143. *
  144. *     Segment des points en base B:
  145. *     -----------------------------
  146.       SEGMENT NCPR(XCOOR(/1)/(IDIM+1))
  147. *     NCRP(#global) = #local dans XPTB (1er indice)      
  148. *
  149. *     Segment des parametres numeriques pour la continuation: 
  150. *     ------------------------------------------------------
  151.       SEGMENT PARNUM
  152.          CHARACTER*4 TYPS
  153.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  154.          REAL*8   PARINI,PARFIN
  155.          INTEGER  ITERMAX,NBPAS
  156.          LOGICAL  JANAL
  157.       ENDSEGMENT
  158. *
  159. *     Segment des resultats:
  160. *     ---------------------
  161.       SEGMENT  PSORT
  162.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  163.         REAL*8 VSAVE(NPAS)
  164.         LOGICAL ZSAVE(NPAS)
  165.         CHARACTER*2 TYPBIF(NBIFU)
  166.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  167.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  168.         INTEGER CBIF
  169.       ENDSEGMENT
  170. *     QSAVE(i,j) = Q harmonique i au pas j
  171. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  172. *     ZSAVE(j)   = stabilite au j-eme pas
  173. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  174. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  175. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  176. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  177. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  178. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  179.  
  180. *     Segment des tableaux de travail:
  181. *     -------------------------------
  182.       SEGMENT  MTEMP
  183.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  184.          REAL*8 T02(NT1+2), TP2(NT1+2)
  185.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  186.          REAL*8 res
  187.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  188.          REAL*8 QOLD(NT1),OMEGOLD
  189.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  190.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  191.       ENDSEGMENT
  192. *     Jacobiennes augmentees 
  193. *     Ja : [ RX  Rw ; dX dw]
  194. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  195.  
  196. *      SEGMENT  NNNN
  197. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  198. *      ENDSEGMENT
  199.  
  200. *************************** fin TMDYNC.INC *****************************
  201.  
  202. -INC CCREEL
  203. -INC SMLREEL
  204.  
  205. C Fonctions BLAS/LAPACK
  206.        REAL*8    DNRM2
  207.        EXTERNAL  DNRM2
  208.  
  209. *----- Recup + Init ----------------------------------------------------
  210.         segini,mwork
  211. *       Variables generalisees: coefficients de Fourier et frequence
  212.         MTQ=KTQ
  213. *       Matrices XK,KASM,XM,GAM,IGAM,DL      
  214.         MTKAM=KTKAM        
  215. *       Deformees modales
  216.         MTPHI = KTPHI   
  217. *       Liaisons sur base A
  218.         MTLIAA = KTLIAA
  219. *       Liaisons sur base B
  220.         MTLIAB = KTLIAB
  221. *       Forces externes
  222.         MTFEX=KTFEX    
  223. *       Truc local base A
  224.         LOCLFA = KOCLFA
  225. *       Truc local base B
  226.         LOCLB1 = KOCLB1         
  227. *       Inconnues sur un pas de temps (AFT)         
  228.         MTPAS=KTPAS  
  229.  
  230. *------ Calcul du produit Z*X ------------------------------------------
  231.         CALL HBMZX(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,ZX,.not.NNM)
  232.  
  233. ****** Cas Modes non-lineaires ****************************************
  234.        IF (NNM) THEN
  235. c         forces ext = 0
  236.           DO I= 1,NT
  237.             PP(I) = 0.D0
  238.           ENDDO        
  239.  
  240. ****** Cas Reponse Forcee et Autonome *********************************
  241.        ELSE        
  242.  
  243. c         CALL HBMZX(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,ZX)
  244. *       Si presence de forces de couplage, on enrichit Z*X 
  245.         DO J = 1,IPALA(/1)
  246.           IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
  247.            DO IJ = 1,7
  248.              VCTCS(IJ) = XPALA(J,IJ)
  249.            ENDDO
  250.            JG = NDDL
  251.            MLREE1 = IPALA(J,7)
  252. c            SEGACT, MLREE1
  253.            DO IJ = 1,JG
  254.              LAMBD(IJ) = MLREE1.PROG(IJ)
  255.            ENDDO
  256.            MLREE2 = IPALA(J,4)
  257.            MLREE3 = IPALA(J,8)
  258. c            SEGACT, MLREE2,MLREE3
  259. c *          Nombre de termes fixe, à generaliser si besoin.            
  260. c               JG = 2
  261. *          tentative de generalisation
  262.            JG = MLREE2.PROG(/1)
  263.            DO IJ = 1,JG
  264.              AiDi(IJ) = MLREE2.PROG(IJ)
  265.              Di(IJ) = MLREE3.PROG(IJ)
  266.            ENDDO
  267. c            SEGDES,MLREE1,MLREE2,MLREE3     
  268.              CALL HBMZFX(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,ZX)
  269.           ENDIF
  270.         ENDDO    
  271. *        DO I = 1, NT
  272. *          WRITE(*,*) 'ZZ(',I,',:)=',(ZZ(I,IOU),IOU=1,NT)
  273. *        ENDDO
  274. *       P: donnee dans MTFEX --> FEXA 
  275. *       Si excitation de type balourd, on met P à jour: P = w²FEXA
  276.         DO I=1,NT
  277.           IF (BAL.EQ.1) THEN
  278.             PP(I) = FEXA(I)*(OMEG**2)
  279.           ELSE
  280.             PP(I) = FEXA(I)
  281.           END IF
  282.         ENDDO
  283.  
  284.        ENDIF
  285. ****** Fin distinction Modes non-lineaires / autres cas  ***************
  286.  
  287.  
  288. *------ Calcul de FNL(X) et dFNLdX (AFT) -------------------------------
  289.  
  290. *  --- Calcul analytique ---
  291.        IF (JANA2) THEN
  292. *       FNL(X) et JAC du segment MTQ 
  293.         CALL HBMFN2(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  294.      &                    MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNL)
  295.  
  296. *  --- Calcul par difference finies ---
  297.        ELSE
  298. *       FNL(X)  
  299.         CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  300.      &                   MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNL)
  301. *       dFNLFdX (Differences finies)
  302. c       Magnitude de la perturbation
  303.         summ = dnrm2(NT,Q1,1)
  304.         IF (summ.NE.0.) THEN
  305. ctest       EPS = summ*((XZPREC)**0.3)
  306.             EPS = 1.D-6 * summ
  307.         ELSE 
  308.             EPS = (XZPREC)**0.3
  309.         ENDIF
  310. c         if(OMEG.gt.0.153.and.OMEG.lt.0.170) write(*,*)'>>> EPS=',EPS
  311. c       EPS = max(XZPREC,summ*1.E-6)
  312. c       Boucle sur les inconnues
  313.         DO J = 1,NT
  314. c          Perturbation des coefficients de Fourier
  315.            Q1SAVE=Q1(J)
  316. c          forward differentation -> JAC+
  317.            Q1(J) =Q1SAVE+EPS
  318.            CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  319.      &                      MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNLP)
  320. cdebug
  321. c             if(OMEG.gt.0.153.and.OMEG.lt.0.170) 
  322. c      &      write(*,123) J,(FNLP(iou),iou=1,NT)
  323.            Q1(J) =Q1SAVE
  324.            DO K=1,NT
  325.               JAC(K,J) =  FNLP(K)
  326.            ENDDO
  327. c          backward differentation -> JAC+
  328.            Q1(J) =Q1SAVE-EPS
  329.            CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  330.      &                      MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNLP)
  331.            Q1(J) =Q1SAVE
  332. c          central differentation -> JAC = (JAC+ + JAC-) / 2
  333.            DO K=1,NT
  334.               JAC(K,J) =  (JAC(K,J) - FNLP(K))/(2.*EPS)
  335.            ENDDO
  336.         ENDDO  
  337.        ENDIF
  338.  
  339.  
  340. *------ Residu d'equilibre : R(X) = ZX-Fnl-Fext ------------------------
  341.  
  342.         DO N=1,NT 
  343.            fvec(N) = ZX(N)-FNL(N)-PP(N)
  344.         ENDDO 
  345. *       reference (maxi entre la norme 2 des 3 termes du residu)
  346.         fref = dnrm2(NT,ZX,1)
  347.         fref = MAX(fref,dnrm2(NT,FNL,1))
  348.         fref = MAX(fref,dnrm2(NT,PP,1))
  349. c         fref = dnrm2(NT,ZX,1) + dnrm2(NT,FNL,1) + dnrm2(NT,PP,1)
  350.  
  351. *       Menage        
  352.         segsup,mwork
  353. cdebug
  354. c  123    FORMAT('Fnl(Q+eps*e_',I2,')=',9(E10.3E1,1X))
  355.  
  356.         RETURN
  357.         END
  358.  
  359.  
  360.  

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