Sources Esope | Cast3M

devpas
C DEVPAS    SOURCE    BP208322  20/09/18    21:15:40     10718          
*
      SUBROUTINE DEVPAS(NA1,NPC1,XK,XASM,XM,PDT,T,NPAS,FTOTA,
     &   FEXA,NPFEXA,NLIAA,NLSA,IPALA,IPLIA,XPALA,XVALA,
     &   NLIAB,NLSB,NPLB,IDIMB,IPALB,IPLIB,JPLIB,XPALB,XVALB,FTOTB,
     &   FTOTBA,XPTB,FEXPSM,
     &   FINERT,IERRD,FTEST,FTEST2,KTQ,
     &   XABSCI,XORDON,NIP,FTEXB,FEXB,RIGIDE,KTPHI,XCHPFB,XOPM1,NB1,
     &   NB1K,NB1C,NB1M)
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*--------------------------------------------------------------------*
*                                                                    *
*     Opérateur DYNE : algorithme de Fu - de Vogelaere               *
*     ________________________________________________               *
*                                                                    *
*     Calcul d'un pas de temps, appel aux s-p spécifiques.           *
*                                                                    *
*     Paramètres:                                                    *
*                                                                    *
* es  Q1(,)   Vecteur des déplacements généralisés                   *
* es  Q2(,)   Vecteur des vitesses généralisées                      *
* es  Q3(,)   Vecteur des accélérations généralisées                 *
* es  NA1     Nombre total d'inconnues en base A                     *
* es  NPC1    Nombre de pas de calcul                                *
* es  XK      Vecteur des raideurs généralisées                      *
* es  XASM    Vecteur des amortissements généralisés                 *
* es  XM      Vecteur des masses généralisées                        *
* e   PDT     pas de temps                                           *
* e   T       temps                                                  *
* e   NPAS    Numéro du pas de temps                                 *
* es  FTOTA   Forces extérieures totalisées, sur la base A           *
* es  FEXA    Evolution des forces extérieures en base A             *
* e   FTEXB   Evolution des forces extérieures en base B             *
* e   FEXB    Forces extérieures sur la base B, servant au calcul    *
*             des moments pour les corps rigides.                    *
* e   RIGIDE  Vrai si corps rigide, faux sinon                       *
* es  IFEXA   Numéro du mode correspondant au point de chargement    *
*             (supprime le 2018-12-14 par bp)
* es  NPFEXA  Nombre de points de chargement                         *
* e   NLIAA   Nombre de liaisons sur la base A                       *
* e   NLSA    Nombre de liaisons A en sortie                         *
* e   IPALA   Tableau renseignant sur le type de liaison (base A)    *
* e   IPLIA   Tableau contenant les numéros "DYNE" des points        *
* e   XPALA   Tableau contenant les paramètres des liaisons          *
* es  XVALA   Tableau contenant les variables internes de liaison A  *
*     XPHILB  Vecteur des deformees modales                          *
* e   NLIAB   Nombre de liaisons sur la base B                       *
* e   NLSB    Nombre de liaisons base B en sortie                    *
* e   NPLB    Nombre total de points de liaisons (base B)            *
* e   IDIMB   Nombre de directions                                   *
* e   IPALB   Tableau renseignant sur le type de liaison             *
* e   IPLIB   Tableau contenant les numeros "DYNE" des points        *
* e   JPLIB   Tableau contenant les numeros "GIBI" des points        *
* e   XPALB   Tableau contenant les parametres des liaisons (base B) *
* es  XVALB   Tableau contenant les variables internes de liaison B  *
*     FTOTB   Forces exterieures totalisees sur la base B            *
* e   XABSCI  Tableau contenant les abscisses de la loi plastique    *
* e   XORDON  Tableau contenant les ordonnees de la loi plastique    *
* e   NIP     Nbr de points dans l'evolution de la loi plastique     *
*     FTOTBA  Forces totales base B projetees base A                 *
*     XPTB    Deplacements des points de liaison                     *
*     IBASB   Appartenance des points de liaison a une sous-base     *
*     IPLSB   Position du point de liaison dans la sous-base         *
*     INMSB   Nombre de modes dans la sous-base                      *
*     IORSB   Position du 1er mode de la sous-base dans ens. modes   *
*     IAROTA  Indique la position des modes de rotation              *
*     NSB     Nombre de sous-bases                                   *
*     NPLSB   Nombre de points de liaison par sous-base              *
*     NA2     Nombre d'inconnues dans la sous-base                   *
*     FEXPSM  Pseudo-Modes base B                              *
*     FINERT  Forces d'inertie base B                                *
*     IERRD   Indicateur d'erreur                                    *
* -   FTEST   Tableau local FTEST de la subroutine DEVLFA            *
* -   FTEST2  Tableau local FTEST de la subroutine DEVLB1            *
* e,s WEXT   travail des forces exterieures
* e,s WINT   travail des forces interieures (rigidite et
*            amortissement et forces de liaison )
*                                                                    *
*     Auteur, date de création:                                      *
*                                                                    *
*     Denis ROBERT-MOUGIN, le 26 mai 1989.                           *
*                                                                    *
*--------------------------------------------------------------------*
*
      SEGMENT,MTPHI
         INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
         INTEGER IAROTA(NSB)
         REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
      ENDSEGMENT
*
      SEGMENT,MTQ
         REAL*8 Q1(NA1,4),Q2(NA1,4),Q3(NA1,4)
         REAL*8 WEXT(NA1,2),WINT(NA1,2)
      ENDSEGMENT
*
c       INTEGER IFEXA(*),IPALA(NLIAA,*),IPLIA(NLIAA,*)
      INTEGER IPALA(NLIAA,*),IPLIA(NLIAA,*)
      INTEGER IPALB(NLIAB,*),IPLIB(NLIAA,*),JPLIB(*)
*      REAL*8 Q1(NA1,*),Q2(NA1,*),Q3(NA1,*)
      REAL*8 XVALA(NLIAA,4,*),XPALA(NLIAA,*),XM(NA1,*),XK(NA1,*)
      REAL*8 XPALB(NLIAB,*),XVALB(NLIAB,4,*),FEXPSM(NPLB,NPC1,2,*)
      REAL*8 XASM(NA1,*),FTOTA(NA1,*),FEXA(NPFEXA,NPC1,*)
      REAL*8 FTOTB(NPLB,*),FTOTBA(*),XPTB(NPLB,2,*),FINERT(NA1,*)
*      REAL*8 WEXT(NA1,2),WINT(NA1,2)
      REAL*8 XABSCI(NLIAB,*),XORDON(NLIAB,*),FEXB(NPLB,2,*)
      REAL*8 FTEST(NA1,4), FTEST2(NPLB,6)
      REAL*8 FTEXB(NPLB,NPC1,2,*),XCHPFB(2,NLIAB,4,NPLB)
      REAL*8 XOPM1(NB1,NB1,*),FAMOR(NA1,4)
*
      LOGICAL RIGIDE
c       LOGICAL LWRITE
c       LWRITE=.FALSE.
c       LWRITE=(NPAS.le.10).or.(mod(NPAS,100).eq.0)
c       if(LWRITE) write(*,*) '---- DEVPAS : PAS',NPAS
c       if(NPAS.eq.1) write(*,*) 'NB1*=',NB1K,NB1C,NB1M,NB1
*
      MTQ = KTQ
      MTPHI  = KTPHI
      NSB = XPHILB(/1)
      NPLSB = XPHILB(/2)
      NA2 =  XPHILB(/3)
      IVINIT = 0
      PDTS2 = 0.5D0 * PDT
 
 
*--------------------------- 1er demi-pas --------------------------*
      III = 2
*
*     Force amortissement généralisées pour le premier demi-pas de temps
*     FAMOR(,4) = C * \dot{q}_-1/2
      CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,4)
*     FAMOR(,3) = C * \dot{q}_0
      CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,3)
*
*     Déplacements généralisés pour le premier demi-pas de temps
*     Q_1/2  =  Q_0 + dt/2 * \dot Q_0
*            + dt^2/24 * M^-1 * (4F_0 - F_-1/2 - 4FAMOR_0 + FAMOR_-1/2)
      IF(NB1.NE.1) THEN
        CALL DEVEQ5(Q1,Q2,NA1,XM,PDT,NPAS,FTOTA,FAMOR,XOPM1,NB1,NB1M)
      ELSE
        CALL DEVEQ1(Q1,Q2,NA1,XASM,XM,PDT,NPAS,FTOTA,FINERT)
      ENDIF
c       if(LWRITE) write(*,*) 'Q1(:,2) =',(Q1(iou,2),iou=1,NA1)
*
*     Totalisation des forces extérieures pour la base A
*     pour la fin du pas :
*     F_1/2 = FEXT_1/2  et F_1 = FEXT_1
*      CALL DEVFXA(FEXA,IFEXA,FTOTA,NPFEXA,NA1,NPC1,NPAS,FTEXB,FEXB,
      CALL DEVFXA(FEXA,FTOTA,NPFEXA,NA1,NPC1,NPAS,FTEXB,FEXB,
     &                  NPLB,IDIMB,RIGIDE)
*
*     Ajout des forces de raideur a l'issue du premier demi-pas
*     F_1/2 = FEXT_1/2 - K Q_1/2
      CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,III)
*
*     estimation de la vitesse (ici plutot que dans DEVLF*, bp,2020-09):
*     \dot{q}_1/2 ~ ({q}_1/2 - {q}_0) / dt/2
*     on utilise Q2(I,III=2) temporairement (il sera "bien" rempli par DEVEQ2/6)
      DO I=1,NA1
        Q2(I,III)=(Q1(I,III)-Q1(I,3))/PDTS2
      ENDDO
c       if(LWRITE) write(*,*) 'Q2(:,2) ~',(Q2(iou,2),iou=1,NA1)
 
*     Ajout des forces de liaison
*     F_1/2 = FEXT_1/2 + FLIAI_1/2
      IF (NLIAA.NE.0) THEN
         CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
     &               NLIAA,PDT,T,NPAS,III,FINERT,IVINIT,FTEST)
      ENDIF
c       if(LWRITE) write(*,*) 'FTOTA_a(:,2) =',(FTOTA(iou,2),iou=1,NA1)
      IF (NLIAB.NE.0) THEN
         CALL DEVLFB(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB,NLIAB,
     &               XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,PDT,T,
     &               NPAS,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,III,
     &               FEXPSM,NPC1,IERRD,FTEST2,
     &               XABSCI,XORDON,NIP,FEXB,RIGIDE,IAROTA,XCHPFB)
         IF (IERRD.NE.0) RETURN
      ENDIF
c       if(LWRITE) write(*,*) 'FTOTA_a+b(:,2) =',(FTOTA(iou,2),iou=1,NA1)
*
*     Vitesses généralisées pour le premier demi-pas de temps
*     \dot{q}_1/2 = ...
      IF(NB1.NE.1) THEN
        CALL DEVEQ6(Q2,NA1,XM,PDT,NPAS,FTOTA,FAMOR,
     &              XOPM1,NB1,NB1M)
      ELSE
        CALL DEVEQ2(Q2,NA1,XASM,XM,PDT,NPAS,FTOTA,FINERT)
      ENDIF
 
*     FAMOR(,2) = C * \dot{q}_1/2
      CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,2)
 
*     accelerations généralisées pour le premier demi-pas de temps
*     \ddot{q}_1/2 = M^-1 * (F_1/2 - FAMOR_0)
c    -cas M pleine
      IF(NB1M.NE.1) THEN
        DO 12 I  = 1,NA1
           Q3(I,2) = 0.D0
        DO 13 IB = 1,NB1
           Q3(I,2) = Q3(I,2)
     &     + XOPM1(I,IB,3)*(FTOTA(IB,III) - FAMOR(IB,III))
 13     CONTINUE
 12     CONTINUE
c    -cas M diagonal
      ELSE
        DO 10 I = 1,NA1
           Q3(I,2) = ( FTOTA(I,III) - FAMOR(I,III) )
     &             / ( XM(I,1) - FINERT(I,III) )
 10     CONTINUE
      ENDIF
c       if(LWRITE) write(*,*) 'Q3(:,2) =',(Q3(iou,2),iou=1,NA1)
 
c     calcul des travaux pour le premier demi-pas de temps
      CALL DEVENE (NA1,III,NPAS,FEXA,Q1,Q2,FTOTA,WEXT,WINT,
     &              FAMOR,NPC1)
 
 
*--------------------------- 2nd demi-pas --------------------------*
      III = 1
 
*     Déplacements généralisés pour le second demi-pas de temps
*     q_1 = ...
      IF(NB1.NE.1) THEN
        CALL DEVEQ7(Q1,Q2,NA1,XM,PDT,NPAS,FTOTA,FAMOR,XOPM1,NB1,NB1M)
      ELSE
        CALL DEVEQ3(Q1,Q2,NA1,XASM,XM,PDT,NPAS,FTOTA,FINERT)
      ENDIF
c       if(LWRITE) write(*,*) 'Q1(:,1) =',(Q1(iou,1),iou=1,NA1)
*
*     Ajout des forces de raideur a l'issue du deuxième demi-pas
*     F_1 = FEXT_1 - K Q_1
      CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,III)
*
*     estimation de la vitesse (ici plutot que dans DEVLF*, bp,2020-09):
*     \dot{q}_1 ~ ({q}_1 - {q}_1/2) / dt/2
*     on utilise Q2(I,III=1) temporairement (il sera "bien" rempli par DEVEQ4/8)
      DO I=1,NA1
        Q2(I,III)=(Q1(I,III)-Q1(I,2))/PDTS2
      ENDDO
c       if(LWRITE) write(*,*) 'Q2(:,1) ~',(Q2(iou,1),iou=1,NA1)
 
*     Ajout des forces de liaison
*     F_1 = FEXT_1 + FLIAI_1
      IF (NLIAA.NE.0) THEN
         CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
     &               NLIAA,PDT,T,NPAS,III,FINERT,IVINIT,FTEST)
      ENDIF
c       if(LWRITE) write(*,*) 'FTOTA_a(:,1) =',(FTOTA(iou,1),iou=1,NA1)
      IF (NLIAB.NE.0) THEN
         CALL DEVLFB(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB,NLIAB,
     &               XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,PDT,T,
     &               NPAS,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,III,
     &               FEXPSM,NPC1,IERRD,FTEST2,
     &               XABSCI,XORDON,NIP,FEXB,RIGIDE,IAROTA,XCHPFB)
         IF (IERRD.NE.0) RETURN
      ENDIF
c       if(LWRITE) write(*,*) 'FTOTA_a+b(:,1) =',(FTOTA(iou,1),iou=1,NA1)
*
*     Vitesses généralisées pour le second demi-pas de temps
*     \dot{q}_1 = ...
      IF(NB1.NE.1) THEN
        CALL DEVEQ8(Q2,NA1,XM,PDT,NPAS,FTOTA,FAMOR,XOPM1,NB1,NB1M)
      ELSE
        CALL DEVEQ4(Q2,NA1,XASM,XM,PDT,NPAS,FTOTA,FINERT)
      ENDIF
*
*     accelerations généralisees pour le second demi-pas de temps
*     \dodt{q}_1 = M 1 * (F_1 - FAMOR_1)
      CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,III)
c    -cas M (ou C) pleine
      IF(NB1M.NE.1) THEN
        DO 22 I  = 1,NA1
           Q3(I,1) = 0.D0
        DO 23 IB = 1,NB1
           Q3(I,1) = Q3(I,1)
     &     + XOPM1(I,IB,3)*(FTOTA(IB,III) - FAMOR(IB,III))
 23     CONTINUE
 22     CONTINUE
c    -cas M diagonal
      ELSE
        DO 20 I = 1,NA1
           Q3(I,1) = ( FTOTA(I,III) - FAMOR(I,III) )
     &             / ( XM(I,1) - FINERT(I,III) )
 20     CONTINUE
      ENDIF
c       if(LWRITE) write(*,*) 'Q3(:,1) =',(Q3(iou,1),iou=1,NA1)
 
c     calcul des travaux pour le premier demi-pas de temps
      CALL DEVENE (NA1,III,NPAS,FEXA,Q1,Q2,FTOTA,WEXT,WINT,
     &              FAMOR,NPC1)
*
      END