Télécharger la plaquette Cast3M


fr en

Télécharger devini.eso

Retour à la liste

Numérotation des lignes : 

devini

  1. C DEVINI    SOURCE    BP208322  20/09/18    21:15:34     10718          
  2.        SUBROUTINE DEVINI(ITINIT,KTKAM,KTQ,KTFEX,KTPAS,KTNUM,KTLIAA,
  3.      &                  KTLIAB,KTPHI,KCPR,KOCLFA,KOCLB1,REPRIS,
  4.      &                  RIGIDE,lmodyn)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. *                                                                    *
  9. *     Operateur DYNE : algorithme de Fu - de Vogelaere               *
  10. *     ________________________________________________               *
  11. *                                                                    *
  12. *     Initialisation de l'algorithme ou reprise de calcul.           *
  13. *                                                                    *
  14. *     Parametres:                                                    *
  15. *                                                                    *
  16. * e   ITINIT  Table contenant les conditions initiales ou les        *
  17. *             champs necessaires a la reprise du calcul              *
  18. * es  KTKAM   Segment contenant les vecteurs XK, XASM et XM          *
  19. * es  KTQ     Segment des variables de mouvement generalisees        *
  20. *             (et des travaux)
  21. * es  KTFEX   Segment contenant les chargements libres               *
  22. * es  KTPAS   Segment des variables au cours d'un pas de temps       *
  23. * es  KTNUM   Segment contenant les parametres numeriques            *
  24. * e   KTLIAA  Segment descriptif des liaisons en base A              *
  25. * e   KTLIAB  Segment descriptif des liaisons en base B              *
  26. * e   KTPHI   Segment contenant les deformees modales                *
  27. * e   KCPR    Segment des points                                     *
  28. * -   KOCLFA  Segment contenant les tableaux locaux de la subroutine *
  29. *             DEVLFA                                                 *
  30. * -   KOCLB1  Segment contenant les tableaux locaux de la subroutine *
  31. *             DEVLB1                                                 *
  32. * e   REPRIS  Vrai si reprise de calcul                              *
  33. * e   RIGIDE  Vrai si corps rigide
  34. *                                                                    *
  35. *     Auteur, date de creation:                                      *
  36. *                                                                    *
  37. *     Denis ROBERT-MOUGIN, le 25 mai 1989.                           *
  38. *                                                                    *
  39. *--------------------------------------------------------------------*
  40.  
  41. -INC PPARAM
  42. -INC CCOPTIO
  43. -INC CCNOYAU
  44. -INC SMCOORD
  45. -INC SMTABLE
  46. -INC CCASSIS
  47. *
  48.       SEGMENT,MTQ
  49.          REAL*8 Q1(NA1,4),Q2(NA1,4),Q3(NA1,4)
  50.          REAL*8 WEXT(NA1,2),WINT(NA1,2)
  51.       ENDSEGMENT
  52.       SEGMENT,MTFEX
  53.          REAL*8  FEXA(NPFEXA,NPC1,2)
  54.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  55.          REAL*8  FTEXB(NPLB,NPC1,2,IDIM)
  56. *         INTEGER IFEXA(NPFEXA),IFEXB(NPFEXB)
  57.       ENDSEGMENT
  58.       SEGMENT,MTPAS
  59.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  60.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  61.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  62.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  63.       ENDSEGMENT
  64.       SEGMENT,MTNUM
  65.          REAL*8 XDT(NPC1),XTEMPS(NPC1)
  66.       ENDSEGMENT
  67.       SEGMENT,MTKAM
  68.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  69.          REAL*8 XOPER(NB1,NB1,NOPER)
  70.       ENDSEGMENT
  71.       SEGMENT,MTLIAA
  72.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  73.          REAL*8  XPALA(NLIAA,NXPALA)
  74.       ENDSEGMENT
  75.       SEGMENT,MTLIAB
  76.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  77.          REAL*8  XPALB(NLIAB,NXPALB)
  78.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  79.       ENDSEGMENT
  80.       SEGMENT,MTPHI
  81.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  82.          INTEGER IAROTA(NSB)
  83.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  84.       ENDSEGMENT
  85.       SEGMENT,ICPR(nbpts)
  86.  
  87. * Segment "local" pour DEVLFA ...
  88.       SEGMENT,LOCLFA
  89.          REAL*8 FTEST(NA1,4)
  90.       ENDSEGMENT
  91. * Segment "local" pour DEVLB1 ...
  92.       SEGMENT,LOCLB1
  93.          REAL*8 FTEST2(NPLB,6)
  94.       ENDSEGMENT
  95.       segment mwinit
  96.         integer jpdep,jpvit,jrepr
  97.       endsegment
  98. c segment local
  99.       SEGMENT   MAMOR
  100.         REAL*8  FAMOR(NA1,4)
  101.         REAL*8  FTEMP(NA1),QTEMP(NA1)
  102.       ENDSEGMENT
  103. c segment local : calcule un operateur et son inverse
  104.       SEGMENT   MOP
  105.         REAL*8  XOP(NB1,NB1)
  106.         INTEGER INVOP(NB1)
  107.         REAL*8  XOPM1(NB1,NB1)
  108.       ENDSEGMENT
  109.  
  110. *
  111.       LOGICAL L0,L1,REPRIS,RIGIDE,lmodyn
  112.       CHARACTER*8 TYPRET,TYPIND,CHARRE
  113.       CHARACTER*(19) CHAI1
  114.       ILC1 = 19
  115.       DATA CHAI1 /'VARIABLES_LIAISON_A'/
  116. *
  117.       MTFEX  = KTFEX
  118.       MTPAS  = KTPAS
  119.       MTNUM  = KTNUM
  120.       MTKAM  = KTKAM
  121.       MTQ    = KTQ
  122.       MTLIAA = KTLIAA
  123.       MTLIAB = KTLIAB
  124.       MTPHI  = KTPHI
  125.       LOCLFA = KOCLFA
  126.       LOCLB1 = KOCLB1
  127.       IDEPL  = 0
  128.       IVITE  = 0
  129.       ITABL  = 0
  130.       ICH1   = 0
  131.       ICH2   = 0
  132.       ICH3   = 0
  133.       ICH4   = 0
  134.       ICH5   = 0
  135.       ICH6   = 0
  136.       NA1    = Q1(/1)
  137.       NB1K   = XK(/2)
  138.       NB1C   = XASM(/2)
  139.       NB1M   = XM(/2)
  140.       NB1    = XOPER(/1)
  141.       NLIAA  = IPALA(/1)
  142.       NLIAB  = IPALB(/1)
  143.       NPLB   = JPLIB(/1)
  144.       NSB    = XPHILB(/1)
  145.       NPLSB  = XPHILB(/2)
  146.       NA2    = XPHILB(/3)
  147.       IDIMB  = XPHILB(/4)
  148.       NPFEXA = FEXA(/1)
  149.       NPC1   = FEXPSM(/2)
  150.       NIP    = XABSCI(/2)
  151.       IERRD = 0
  152. *
  153. *     S'agit-il d'une initialisation ou d'une reprise de calcul ? ------
  154. *
  155.       IF ( REPRIS ) THEN
  156.          IF (IIMPI.EQ.333) THEN
  157.             WRITE(IOIMP,*)'DEVINI : ---> reprise de calcul'
  158.          ENDIF
  159.          if (lmodyn) then
  160.            mwinit = itinit
  161.            segact mwinit
  162.            ITABL = jrepr
  163.          else
  164.          CALL ACCTAB(ITINIT,'MOT',I0,X0,'REPRISE',L0,IP0,
  165.      &                    'TABLE',I1,X1,' ',L1,ITABL)
  166.          endif
  167. *                                               .             .
  168. *        Reprise du calcul, on remplit: Q       Q        Q    Q
  169. *                                        i,-1/2   i,-1/2  i,0  i,0
  170. *
  171.          CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT',L0,IP0,
  172.      &                 'CHPOINT',I1,X1,' ',L1,ICH1)
  173.          IF (IERR.NE.0) RETURN
  174. *
  175.          CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE',L0,IP0,
  176.      &                 'CHPOINT',I1,X1,' ',L1,ICH2)
  177.          IF (IERR.NE.0) RETURN
  178. *
  179.          CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT_1/2',L0,IP0,
  180.      &                 'CHPOINT',I1,X1,' ',L1,ICH3)
  181.          IF (IERR.NE.0) RETURN
  182. *
  183.          CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE_1/2',L0,IP0,
  184.      &                 'CHPOINT',I1,X1,' ',L1,ICH4)
  185.          IF (IERR.NE.0) RETURN
  186. *
  187.          CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A',L0,IP0,
  188.      &                 'CHPOINT',I1,X1,' ',L1,ICH5)
  189.          IF (IERR.NE.0) RETURN
  190. *
  191.          CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A_1/2',L0,IP0,
  192.      &                 'CHPOINT',I1,X1,' ',L1,ICH6)
  193.          IF (IERR.NE.0) RETURN
  194. *
  195.          CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION',L0,IP0,
  196.      &                 'CHPOINT',I1,X1,' ',L1,ICH7)
  197.          IF (IERR.NE.0) RETURN
  198. *
  199.          CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION_1/2',L0,IP0,
  200.      &                 'CHPOINT',I1,X1,' ',L1,ICH8)
  201.          IF (IERR.NE.0) RETURN
  202. *
  203.          CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_EXTERIEUR',L0,IP0,
  204.      &                 'CHPOINT',I1,X1,' ',L1,ICH9)
  205.          IF (IERR.NE.0) RETURN
  206. *
  207.          CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_INTERIEUR',L0,IP0,
  208.      &                 'CHPOINT',I1,X1,' ',L1,ICH10)
  209.          IF (IERR.NE.0) RETURN
  210. *
  211. c        Q1(3)=q_0   
  212.          IF (ICH1.NE.0) THEN
  213.             CALL DYNE18(ICH1,KTQ,1,3,KCPR)
  214.          ELSE
  215.             CALL ERREUR(487)
  216.             RETURN
  217.          ENDIF
  218. *
  219. c        Q2(3)=\dot{q}_0
  220.          IF (ICH2.NE.0) THEN
  221.             CALL DYNE18(ICH2,KTQ,2,3,KCPR)
  222.          ELSE
  223.             CALL ERREUR(487)
  224.             RETURN
  225.          ENDIF
  226. *
  227.          IF (ICH3.NE.0) THEN
  228.             CALL DYNE18(ICH3,KTQ,1,4,KCPR)
  229.          ELSE
  230.             CALL ERREUR(487)
  231.             RETURN
  232.          ENDIF
  233. *
  234.          IF (ICH4.NE.0) THEN
  235.             CALL DYNE18(ICH4,KTQ,2,4,KCPR)
  236.          ELSE
  237.             CALL ERREUR(487)
  238.             RETURN
  239.          ENDIF
  240.  
  241. *        FTOTA(3)=F_0
  242.          IF (ICH5.NE.0) THEN
  243.             CALL DYNE23(ICH5,KTPAS,3)
  244.          ELSE
  245.             CALL ERREUR(487)
  246.             RETURN
  247.          ENDIF
  248.  
  249. *        FTOTA(4)=F_-1/2
  250.          IF (ICH6.NE.0) THEN
  251.             CALL DYNE23(ICH6,KTPAS,4)
  252.          ELSE
  253.             CALL ERREUR(487)
  254.             RETURN
  255.          ENDIF
  256. *
  257.          IF (ICH7.NE.0) THEN
  258.             CALL DYNE18(ICH7,KTQ,3,3,KCPR)
  259.          ELSE
  260.             CALL ERREUR(487)
  261.             RETURN
  262.          ENDIF
  263. *
  264.          IF (ICH8.NE.0) THEN
  265.             CALL DYNE18(ICH8,KTQ,3,4,KCPR)
  266.          ELSE
  267.             CALL ERREUR(487)
  268.             RETURN
  269.          ENDIF
  270. *
  271.          IF (ICH9.NE.0) THEN
  272.             CALL DYNE18(ICH9,KTQ,4,1,KCPR)
  273.          ENDIF
  274.          IF (ICH10.NE.0) THEN
  275.             CALL DYNE18(ICH10,KTQ,5,1,KCPR)
  276.          ENDIF
  277. *
  278. *
  279.          IF (NLIAA.NE.0) THEN
  280. *
  281. *           l'indice VARIABLES_LIAISON_A n'existe que pour
  282. *           les liaisons POLYNOMIALES 
  283. *           et liaison COUPLAGE_DEPLACEMENT avec CONVOLUTION
  284. *
  285.             IPOLY = 0
  286.             MTABLE = ITABL
  287.             SEGACT,MTABLE
  288.             NIND1 = MLOTAB
  289.             if(nbesc.ne.0)segact ipiloc
  290.             DO 100 I=1,NIND1
  291.                TYPIND = MTABTI(I)
  292.                IF (TYPIND.EQ.'MOT     ') THEN
  293.                   IP = MTABII(I)
  294.                   ID = IPCHAR(IP)
  295.                   IFI = IPCHAR(IP+1)
  296.                   IL1 = IFI - ID
  297.                   IF (IL1.EQ.ILC1) THEN
  298.                      IF (CHAI1.EQ.ICHARA(ID:IFI-1)) THEN
  299.                         IPOLY = 1
  300.                      ENDIF
  301.                   ENDIF
  302.                ENDIF
  303.  100        CONTINUE
  304.             if(nbesc.ne.0)SEGDES,IPILOC
  305.             SEGDES,MTABLE
  306.             IF (IIMPI.EQ.333) THEN
  307.                WRITE(IOIMP,*) 'DEVINI : IPOLY = ',IPOLY
  308.             ENDIF
  309.             IF (IPOLY.NE.0) THEN
  310.                CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_A',
  311.      &                     L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  312.                CALL DYNA14(ITREFR,KTLIAA)
  313.                IF (IERR.NE.0) RETURN
  314.             ENDIF
  315.          ENDIF
  316.          IF (NLIAB.NE.0) THEN
  317. *old            CALL DEVRCO(Q1,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB,
  318.             CALL DEVRCO(Q1,Q2,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB,
  319.      &                  IBASB,IPLSB,INMSB,IORSB,3,IAROTA)
  320. *           --> XPTB(:,3,:)=x_0   XPTB(:,4,:)=\dot{q}_0
  321. *               ok car XPTB(:,4,:) utilisé uniquement pour calcul FnlB
  322.             CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_B',
  323.      &                  L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  324.             IF (IERR.NE.0) RETURN
  325.             CALL DYNE14(ITREFR,KTLIAB)
  326.             IF (IERR.NE.0) RETURN
  327.          ENDIF
  328. *
  329.  
  330.       ELSE
  331. *
  332. *        Chpoints initiaux de deplacement et de vitesse: ---------------
  333. *
  334.          IF (ITINIT.NE.0) THEN
  335.           if (lmodyn) then
  336.            mwinit = itinit
  337.            segact mwinit
  338.            idepl = jpdep
  339.            ivite = jpvit
  340.          else
  341.             TYPRET=' '
  342.             CALL ACCTAB(ITINIT,'MOT',I0,X0,'DEPLACEMENT',L0,IP0,
  343.      &                        TYPRET,I1,X1,CHARRE,L1,IDEPL)
  344. *
  345.             TYPRET=' '
  346.             CALL ACCTAB(ITINIT,'MOT',I0,X0,'VITESSE',L0,IP0,
  347.      &                        TYPRET,I1,X1,CHARRE,L1,IVITE)
  348.           endif
  349.          ENDIF
  350. *
  351. *        Mise des valeurs initiales au pas m (indice 3)
  352. c        Q1(3)=q_0   
  353. c        Q2(3)=\dot{q}_0
  354.          IF (IDEPL.NE.0) THEN
  355.             CALL DYNE18(IDEPL,KTQ,1,3,KCPR)
  356.             MTQ = KTQ
  357.          ENDIF
  358.          IVINIT = 0
  359.          IF (IVITE.NE.0) THEN
  360.             CALL DYNE18(IVITE,KTQ,2,3,KCPR)
  361.             IVINIT = 1
  362.          ENDIF
  363. *
  364. *        Ajout des forces de liaison a t=0
  365. *
  366.          PDT=XDT(1)
  367.          T=XTEMPS(1)
  368.          IF (NLIAA.NE.0) THEN
  369.             CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
  370.      &                  NLIAA,PDT,T,1,3,FINERT,IVINIT,FTEST)
  371.          ENDIF
  372.          IF (NLIAB.NE.0) THEN
  373.              IF (RIGIDE) THEN
  374.                 DO 13 IP=1,NPLB
  375.                   DO 15 ID=1,IDIM
  376.                        FEXB(IP,1,ID) = FTEXB(IP,1,1,ID)
  377.  15               CONTINUE
  378.  13             CONTINUE
  379.             ENDIF
  380. c      SEGMENT,MTPHI
  381. c         INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  382. c         INTEGER IAROTA(NSB)
  383. c         REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  384. c      ENDSEGMENT
  385.             CALL DEVLF2(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB,
  386.      &           NLIAB,XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,
  387.      &           PDT,T,1,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,3,
  388.      &           FEXPSM,NPC1,IERRD,FTEST2,XABSCI,XORDON,
  389.      &           NIP,IAROTA,RIGIDE,FEXB,XCHPFB)
  390.              IF (IERRD.NE.0) THEN
  391.                IF (IERRD.EQ.1) CALL ERREUR(528)
  392.                RETURN
  393.             ENDIF
  394.  
  395.          ENDIF
  396. c        fin distinction reprise/CI  -----------------------------------
  397.  
  398.  
  399.          IF (NLIAA.NE.0 .OR. NLIAB.NE.0) THEN
  400.             DO 5 I=1,NA1
  401.                FINERT(I,4) = FINERT(I,3)
  402.  5             CONTINUE
  403. *           end do
  404.          ENDIF
  405. *
  406. *        Calcul des deplacements et des vitesses au pas -1/2 (indice 4)
  407. *
  408. c          PDT = XDT(1)
  409.          PDTS2 = PDT / 2.D0
  410.          AUX2 = PDT * PDT / 8.D0
  411.       IF (IIMPI.EQ.333) THEN
  412.          DO 91 I=1,NA1
  413.             WRITE(IOIMP,*)'DEVINI :'
  414.             WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  415.             WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  416.             WRITE(IOIMP,*)'DEVINI : Fexa (',I,',1) =',Fexa(I,1)
  417.             WRITE(IOIMP,*)'DEVINI : Fexa (',I,',2) =',Fexa(I,2)
  418.  91         CONTINUE
  419. *        end do
  420.       ENDIF
  421.  
  422. c        Cas selon nature (pleine ou diag) des matrices K,C et M :
  423.  
  424. c------- Cas K,C ou M pleine -----------------------------------
  425.          IF(NB1K.GT.1.OR.NB1C.GT.1.OR.NB1M.GT.1) THEN
  426.           IF (IIMPI.EQ.333) WRITE(IOIMP,*)'pleine K,C,M',NB1K,NB1C,NB1M
  427.  
  428. c         F_0 =  F^liaison_0 + F^ext_0 - K * Q_0
  429.           DO 61 I = 1,NA1
  430.              FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1)
  431.  61       CONTINUE
  432.           CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,3)
  433.  
  434. c         F^AMOR_0 = C * \dot Q_0
  435.           SEGINI,MAMOR
  436.           CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,3)
  437.  
  438. c         \ddot Q_0 = M-1 * (F_0 - F^AMOR_0)
  439.           CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,3)
  440.  
  441. c         Q_-1/2    = Q_0 - dt/2 * \dot Q_0 + dt^2/8 *  \ddot Q_0
  442.           DO 62 I = 1,NA1
  443.              Q1(I,4) = Q1(I,3) - PDTS2*Q2(I,3) + AUX2*Q3(I,3)
  444.  62       CONTINUE
  445.  
  446. c         F_-1/2 =  F^liaison_-1/2 + F^ext_-1/2 - K * Q_-1/2
  447.           DO 63 I = 1,NA1
  448.              FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2)
  449.  63       CONTINUE
  450.           CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,4)
  451.  
  452. c         \dot Q_-1/2 = [4I+dt*M-1*C]-1
  453. c                     * ( 4 \dot Q_0 - dt*M-1*(F_-1/2+F_0-F^AMOR_0) )
  454.           DO 64 I = 1,NA1
  455.              FTEMP(I)   = FTOTA(I,4) + FTOTA(I,3) - FAMOR(I,3)
  456.  64       CONTINUE
  457.           CALL DEVLM1(QTEMP,XM,XOPER,FTEMP,NA1,NB1,3,NB1M)
  458.           DO 65 I = 1,NA1
  459.              FTEMP(I)   = 4.D0*Q2(I,3) - PDT*QTEMP(I)
  460.  65       CONTINUE
  461.           IF(NB1C.GT.1.OR.NB1M.GT.1) THEN
  462. *           initialisation de l'operateur (local) et de son inverse
  463. *           XOP = 4I - dt*M^-1*C
  464.             SEGINI,MOP
  465. c           ...C pleine, M pleine
  466.             IF(NB1C.GT.1.AND.NB1M.GT.1) THEN
  467.                DO 662 IA=1,NA1
  468.                   XOP(IA,IA) = 4.D0
  469.                DO 663 IB=1,NA1
  470.                   AUX = 0.D0
  471.                   DO 664 IC=1,NA1
  472.                      AUX = AUX + XOPER(IA,IC,3)*XASM(IC,IB)
  473.  664              CONTINUE
  474.                   XOP(IA,IB) = XOP(IA,IB) - PDT*AUX
  475.  663           CONTINUE
  476.  662           CONTINUE
  477. c           ...C pleine, M diago
  478.             ELSEIF(NB1C.GT.1) THEN
  479.                DO 665 IA=1,NA1
  480.                   XOP(IA,IA) = 4.D0
  481.                DO 666 IB=1,NA1
  482.                   XOP(IA,IB) = XOP(IA,IB) - PDT*XASM(IA,IB)/XM(IA,1)
  483.  666           CONTINUE
  484.  665           CONTINUE
  485. c           ...C diago, M pleine
  486.             ELSEIF(NB1M.GT.1) THEN
  487.                DO 667 IA=1,NA1
  488.                   XOP(IA,IA) = 4.D0
  489.                DO 668 IB=1,NA1
  490.                   XOP(IA,IB) = XOP(IA,IB) - PDT*XOPM1(IA,IB)*XASM(IB,1)
  491.  668           CONTINUE
  492.  667           CONTINUE
  493.             ENDIF
  494.             CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET)
  495.             IF (IIMPI.EQ.333) THEN
  496.               WRITE(IOIMP,*) 'FTEMP(:)=',(FTEMP(jou),jou=1,NA1)
  497.               do iou=1,NB1
  498.                 WRITE(IOIMP,*) 'XOPM1(',iou,',:)=',
  499.      &                         (XOPM1(iou,jou),jou=1,NB1)
  500.               enddo
  501.             ENDIF
  502. c           \dot Q_-1/2 = [4I-dt*M-1*C]-1 *  FTEMP
  503.             DO 67 IA = 1,NA1
  504.                Q2(IA,4) = 0.D0
  505.             DO 68 IB = 1,NA1
  506.                Q2(IA,4) = Q2(IA,4) + XOPM1(IA,IB)*FTEMP(IB)
  507.  68         CONTINUE
  508.  67         CONTINUE
  509.             SEGSUP,MOP
  510.           ELSE
  511.             DO 69 IA = 1,NA1
  512.                Q2(IA,4) = FTEMP(IA) / (4.D0 - (PDT*XASM(IA,1)/XM(IA,1)))
  513.  69         CONTINUE
  514.           ENDIF
  515.  
  516. c         \ddot Q_-1/2 = M-1 * ( F_-1/2 - F^AMOR_-1/2)
  517.           CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,4)
  518.           CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,4)
  519.  
  520.           SEGSUP,MAMOR
  521.  
  522. c------- Cas K,C et M diagonal -----------------------------------
  523.          ELSE
  524.  
  525.           DO 10 I = 1,NA1
  526.             UNSM3 = 1.D0 / ( XM(I,1) - FINERT(I,3) )
  527.             UNSM4 = 1.D0 / ( XM(I,1) - FINERT(I,4) )
  528. *
  529. *           G    =  F    - K  Q
  530. *            i,0     i,0    i  i,0
  531. *
  532.             FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1) - (XK(I,1)*Q1(I,3))
  533.             AUX1 = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  534.             AUX3 = Q1(I,3) - ( PDTS2 * Q2(I,3) )
  535. *
  536. *                                .       2               .
  537. *           Q       = Q    - h/2 Q    + h /8 ( G    - A  Q    )
  538. *            i,-1/2    i,0        i,0           i,0    i  i,0
  539. *
  540.             Q1(I,4) = AUX3 + ( AUX2 * AUX1 )
  541. *
  542. *           G       =  F       - K  Q
  543. *            i,-1/2     i,-1/2    i  i,-1/2
  544. *
  545.             FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2) - (XK(I,1)*Q1(I,4))
  546.             AUX4 = 4.D0 + ( XASM(I,1) * PDT * UNSM3 )
  547.             AUX5 = ( FTOTA(I,4) * UNSM4 ) + ( FTOTA(I,3) * UNSM3 )
  548.             AUX6 = 1.D0 / ( 4.D0 - ( XASM(I,1) * PDT * UNSM3 ) )
  549. *
  550. *        .            1                    .
  551. *        Q       = ________ ( ( 4 + A  h ) Q    - h ( G       + G    ) )
  552. *         i,-1/2   4 - A  h          i      i,0        i,-1/2    i,0
  553. *                       i
  554. *
  555.             Q2(I,4) = AUX6 * ( (AUX4*Q2(I,3))  - (PDT*AUX5) )
  556. *
  557. *           "               .
  558. *           Q  =  ( G  - A  Q  ) / M
  559. *            i       i    i  i      i
  560. *
  561.             Q3(I,3) = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  562.             Q3(I,4) = (FTOTA(I,4)*UNSM4) - (XASM(I,1)*UNSM4*Q2(I,4))
  563.  10       CONTINUE
  564.  
  565.          ENDIF
  566. *
  567.       ENDIF
  568.       IF (IIMPI.EQ.333) THEN
  569.          DO 30 I=1,NA1
  570.             WRITE(IOIMP,*)'DEVINI :'
  571.             WRITE(IOIMP,*)'DEVINI :  Q1(',I,',3) =',Q1(I,3)
  572.             WRITE(IOIMP,*)'DEVINI :  Q2(',I,',3) =',Q2(I,3)
  573.             WRITE(IOIMP,*)'DEVINI :  Q3(',I,',3) =',Q3(I,3)
  574.             WRITE(IOIMP,*)'DEVINI :'
  575.             WRITE(IOIMP,*)'DEVINI :  Q1(',I,',4) =',Q1(I,4)
  576.             WRITE(IOIMP,*)'DEVINI :  Q2(',I,',4) =',Q2(I,4)
  577.             WRITE(IOIMP,*)'DEVINI :  Q3(',I,',4) =',Q3(I,4)
  578.  30         CONTINUE
  579. *        end do
  580.       ENDIF
  581.       IF (IIMPI.EQ.333) THEN
  582.          DO 31 I=1,NA1
  583.             WRITE(IOIMP,*)'DEVINI :'
  584.             WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  585.             WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  586.  31         CONTINUE
  587. *        end do
  588.       ENDIF
  589.  
  590. *
  591.       ICPR = KCPR
  592.       SEGSUP,ICPR
  593. *
  594.       END
  595.  
  596.  
  597.  
  598.  
  599.  
  600.  
  601.  

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