Sources Esope | Cast3M

funcv
C FUNCV     SOURCE    OF166741  26/05/11    21:15:06     12538          
 
      SUBROUTINE funcv(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,KTPHI,KTLIAA,
     &             KTLIAB,KTPAS,KTFEX,KOCLFA,KOCLB1,JANA2,fvec,fref,NNM)
 
*=======================================================================
*       Calcul du residu R(X) = ZX-Fnl-Fext -> fvec = mlreel.prog(*)
*       et de la Jacobienne dFnl/dX         -> JAC(*,*) du segment KTQ
*=======================================================================
 
*----- Declarations ----------------------------------------------------
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC CCREEL
 
-INC SMLREEL
      POINTEUR FNL.MLREEL,FNLP.MLREEL,fvec.MLREEL
 
-INC TMDYNC
 
      SEGMENT mwork
        REAL*8 VCTCS(7),AiDi(2),Di(2),ZX(NT),LAMBD(NDDL),PP(NT)
      ENDSEGMENT
 
      INTEGER NDDL,NFFT,NT,NHBM
c     NNM: cas des modes non-lineaires
 
      REAL*8 EPS, summ, Q1SAVE, fref
      LOGICAL JANA2,NNM
 
C Fonctions BLAS/LAPACK
      REAL*8    DNRM2
      EXTERNAL  DNRM2
 
*----- Recup -----------------------------------------------------------
*  Variables generalisees: coefficients de Fourier et frequence
      MTQ = KTQ
*  Matrices XK,KASM,XM,GAM,IGAM,DL
      MTKAM = KTKAM
*  Deformees modales
      MTPHI = KTPHI
*  Liaisons sur base A
      MTLIAA = KTLIAA
*  Liaisons sur base B
      MTLIAB = KTLIAB
*  Forces externes
      MTFEX = KTFEX
*  Truc local base A
      LOCLFA = KOCLFA
*  Truc local base B
      LOCLB1 = KOCLB1
*  Inconnues sur un pas de temps (AFT)
      MTPAS = KTPAS
 
*----- Init ------------------------------------------------------------
      SEGINI,mwork
      JG = NT
      SEGINI,FNL,FNLP
 
*------ Calcul du produit Z*X ------------------------------------------
      CALL HBMZX(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,mwork.ZX,.not.NNM)
 
****** Cas Modes non-lineaires ****************************************
      IF (NNM) THEN
c       forces ext = 0
        DO I= 1,NT
          PP(I) = 0.D0
        ENDDO
 
****** Cas Reponse Forcee et Autonome *********************************
      ELSE
 
c        CALL HBMZX(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,ZX)
*       Si presence de forces de couplage, on enrichit Z*X
        DO J = 1,IPALA(/1)
          IF (IPALA(J,1).EQ.5 .AND. IPALA(J,3).EQ.101) THEN
           DO IJ = 1,7
             VCTCS(IJ) = XPALA(J,IJ)
           ENDDO
           JG = NDDL
           MLREE1 = IPALA(J,7)
c            SEGACT, MLREE1
           DO IJ = 1,JG
             LAMBD(IJ) = MLREE1.PROG(IJ)
           ENDDO
           MLREE2 = IPALA(J,4)
           MLREE3 = IPALA(J,8)
c            SEGACT, MLREE2,MLREE3
c *          Nombre de termes fixe, a generaliser si besoin.
c               JG = 2
*          tentative de generalisation
           JG = MLREE2.PROG(/1)
           DO IJ = 1,JG
             AiDi(IJ) = MLREE2.PROG(IJ)
             Di(IJ) = MLREE3.PROG(IJ)
           ENDDO
c            SEGDES,MLREE1,MLREE2,MLREE3
             CALL HBMZFX(NT,NDDL,OMEG,Q1,AiDi,Di,LAMBD,VCTCS,ZX)
          ENDIF
        ENDDO
c        DO I = 1, NT
c          WRITE(*,*) 'ZZ(',I,',:)=',(ZZ(I,IOU),IOU=1,NT)
c        ENDDO
*       P: donnee dans MTFEX --> FEXA
*       Si excitation de type balourd, on met P a jour: P = w^2.FEXA
        DO I=1,NT
          IF (BAL.EQ.1) THEN
            PP(I) = FEXA(I)*(OMEG**2)
          ELSE
            PP(I) = FEXA(I)
          END IF
        ENDDO
 
      ENDIF
****** Fin distinction Modes non-lineaires / autres cas  ***************
 
*------ Calcul de FNL(X) et dFNLdX (AFT) -------------------------------
 
*  --- Calcul analytique ---
      IF (JANA2) THEN
*       FNL(X) et JAC du segment MTQ
        CALL HBMFN2(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
     &                    MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNL)
 
*  --- Calcul par difference finies ---
      ELSE
*       FNL(X)
        CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
     &                   MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNL)
*       dFNLFdX (Differences finies)
c       Magnitude de la perturbation
        summ = dnrm2(NT,Q1,1)
        IF (summ.NE.0.) THEN
ctest       EPS = summ*((XZPREC)**0.3)
            EPS = 1.D-6 * summ
        ELSE
            EPS = (XZPREC)**0.3
        ENDIF
c       EPS = max(XZPREC,summ*1.E-6)
c*        if(OMEG.gt.0.153.and.OMEG.lt.0.170) write(*,*)'>>> EPS=',EPS
c*        if(OMEG.gt.0.153.and.OMEG.lt.0.170)
c*     &    write(*,*) 'FUNCV Q1 =',(Q1(iou),iou=1,NT)
        r_eps = 0.5D0 / EPS
c       Boucle sur les inconnues
        DO J = 1,NT
c          Perturbation des coefficients de Fourier
           Q1SAVE=Q1(J)
c          forward differentation -> JAC+
           Q1(J) =Q1SAVE+EPS
           CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
     &                      MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNLP)
c*           if(OMEG.gt.0.153.and.OMEG.lt.0.170)
c*     &       write(*,123) J,(FNLP.PROG(iou),iou=1,NT)
 123    FORMAT('Fnl(Q+eps*e_',I2,')=',9(E13.5,1X))
           Q1(J) =Q1SAVE
           DO K=1,NT
              JAC(K,J) =  FNLP.PROG(K)
           ENDDO
c          backward differentation -> JAC+
           Q1(J) =Q1SAVE-EPS
           CALL HBMFNL(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
     &                      MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,FNLP)
           Q1(J) =Q1SAVE
c          central differentation -> JAC = (JAC+ + JAC-) / 2
           DO K=1,NT
              JAC(K,J) =  (JAC(K,J) - FNLP.PROG(K))*r_eps
           ENDDO
        ENDDO
      ENDIF
 
*------ Residu d'equilibre : R(X) = ZX-Fnl-Fext ------------------------
 
      DO N=1,NT
        fvec.PROG(N) = ZX(N)-FNL.PROG(N)-PP(N)
      ENDDO
*       reference (maxi entre la norme 2 des 3 termes du residu)
      fref = dnrm2(NT,ZX,1)
      fref = MAX(fref,dnrm2(NT,FNL.PROG(1),1))
      fref = MAX(fref,dnrm2(NT,PP,1))
c      fref = dnrm2(NT,ZX,1) + dnrm2(NT,FNL.PROG(1),1) + dnrm2(NT,PP,1)
 
*   Menage
      SEGSUP,mwork
 
      RETURN
      END