C DEVINI SOURCE BP208322 20/09/18 21:15:34 10718 SUBROUTINE DEVINI(ITINIT,KTKAM,KTQ,KTFEX,KTPAS,KTNUM,KTLIAA, & KTLIAB,KTPHI,KCPR,KOCLFA,KOCLB1,REPRIS, & RIGIDE,lmodyn) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *--------------------------------------------------------------------* * * * Operateur DYNE : algorithme de Fu - de Vogelaere * * ________________________________________________ * * * * Initialisation de l'algorithme ou reprise de calcul. * * * * Parametres: * * * * e ITINIT Table contenant les conditions initiales ou les * * champs necessaires a la reprise du calcul * * es KTKAM Segment contenant les vecteurs XK, XASM et XM * * es KTQ Segment des variables de mouvement generalisees * * (et des travaux) * es KTFEX Segment contenant les chargements libres * * es KTPAS Segment des variables au cours d'un pas de temps * * es KTNUM Segment contenant les parametres numeriques * * e KTLIAA Segment descriptif des liaisons en base A * * e KTLIAB Segment descriptif des liaisons en base B * * e KTPHI Segment contenant les deformees modales * * e KCPR Segment des points * * - KOCLFA Segment contenant les tableaux locaux de la subroutine * * DEVLFA * * - KOCLB1 Segment contenant les tableaux locaux de la subroutine * * DEVLB1 * * e REPRIS Vrai si reprise de calcul * * e RIGIDE Vrai si corps rigide * * * Auteur, date de creation: * * * * Denis ROBERT-MOUGIN, le 25 mai 1989. * * * *--------------------------------------------------------------------* -INC PPARAM -INC CCOPTIO -INC CCNOYAU -INC SMCOORD -INC SMTABLE -INC CCASSIS * SEGMENT,MTQ REAL*8 Q1(NA1,4),Q2(NA1,4),Q3(NA1,4) REAL*8 WEXT(NA1,2),WINT(NA1,2) ENDSEGMENT SEGMENT,MTFEX REAL*8 FEXA(NPFEXA,NPC1,2) REAL*8 FEXPSM(NPLB,NPC1,2,IDIMB) REAL*8 FTEXB(NPLB,NPC1,2,IDIM) * INTEGER IFEXA(NPFEXA),IFEXB(NPFEXB) ENDSEGMENT SEGMENT,MTPAS REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1) REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4) REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR) REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB) ENDSEGMENT SEGMENT,MTNUM REAL*8 XDT(NPC1),XTEMPS(NPC1) ENDSEGMENT SEGMENT,MTKAM REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M) REAL*8 XOPER(NB1,NB1,NOPER) ENDSEGMENT SEGMENT,MTLIAA INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA) REAL*8 XPALA(NLIAA,NXPALA) ENDSEGMENT SEGMENT,MTLIAB INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB) REAL*8 XPALB(NLIAB,NXPALB) REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP) ENDSEGMENT SEGMENT,MTPHI INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB) INTEGER IAROTA(NSB) REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB) ENDSEGMENT SEGMENT,ICPR(nbpts) * Segment "local" pour DEVLFA ... SEGMENT,LOCLFA REAL*8 FTEST(NA1,4) ENDSEGMENT * Segment "local" pour DEVLB1 ... SEGMENT,LOCLB1 REAL*8 FTEST2(NPLB,6) ENDSEGMENT segment mwinit integer jpdep,jpvit,jrepr endsegment c segment local SEGMENT MAMOR REAL*8 FAMOR(NA1,4) REAL*8 FTEMP(NA1),QTEMP(NA1) ENDSEGMENT c segment local : calcule un operateur et son inverse SEGMENT MOP REAL*8 XOP(NB1,NB1) INTEGER INVOP(NB1) REAL*8 XOPM1(NB1,NB1) ENDSEGMENT * LOGICAL L0,L1,REPRIS,RIGIDE,lmodyn CHARACTER*8 TYPRET,TYPIND,CHARRE CHARACTER*(19) CHAI1 ILC1 = 19 DATA CHAI1 /'VARIABLES_LIAISON_A'/ * MTFEX = KTFEX MTPAS = KTPAS MTNUM = KTNUM MTKAM = KTKAM MTQ = KTQ MTLIAA = KTLIAA MTLIAB = KTLIAB MTPHI = KTPHI LOCLFA = KOCLFA LOCLB1 = KOCLB1 IDEPL = 0 IVITE = 0 ITABL = 0 ICH1 = 0 ICH2 = 0 ICH3 = 0 ICH4 = 0 ICH5 = 0 ICH6 = 0 NA1 = Q1(/1) NB1K = XK(/2) NB1C = XASM(/2) NB1M = XM(/2) NB1 = XOPER(/1) NLIAA = IPALA(/1) NLIAB = IPALB(/1) NPLB = JPLIB(/1) NSB = XPHILB(/1) NPLSB = XPHILB(/2) NA2 = XPHILB(/3) IDIMB = XPHILB(/4) NPFEXA = FEXA(/1) NPC1 = FEXPSM(/2) NIP = XABSCI(/2) IERRD = 0 * * S'agit-il d'une initialisation ou d'une reprise de calcul ? ------ * IF ( REPRIS ) THEN IF (IIMPI.EQ.333) THEN WRITE(IOIMP,*)'DEVINI : ---> reprise de calcul' ENDIF if (lmodyn) then mwinit = itinit segact mwinit ITABL = jrepr else CALL ACCTAB(ITINIT,'MOT',I0,X0,'REPRISE',L0,IP0, & 'TABLE',I1,X1,' ',L1,ITABL) endif * . . * Reprise du calcul, on remplit: Q Q Q Q * i,-1/2 i,-1/2 i,0 i,0 * CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH1) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH2) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT_1/2',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH3) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE_1/2',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH4) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH5) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A_1/2',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH6) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH7) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION_1/2',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH8) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_EXTERIEUR',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH9) IF (IERR.NE.0) RETURN * CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_INTERIEUR',L0,IP0, & 'CHPOINT',I1,X1,' ',L1,ICH10) IF (IERR.NE.0) RETURN * c Q1(3)=q_0 IF (ICH1.NE.0) THEN CALL DYNE18(ICH1,KTQ,1,3,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * c Q2(3)=\dot{q}_0 IF (ICH2.NE.0) THEN CALL DYNE18(ICH2,KTQ,2,3,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * IF (ICH3.NE.0) THEN CALL DYNE18(ICH3,KTQ,1,4,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * IF (ICH4.NE.0) THEN CALL DYNE18(ICH4,KTQ,2,4,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * FTOTA(3)=F_0 IF (ICH5.NE.0) THEN CALL DYNE23(ICH5,KTPAS,3) ELSE CALL ERREUR(487) RETURN ENDIF * FTOTA(4)=F_-1/2 IF (ICH6.NE.0) THEN CALL DYNE23(ICH6,KTPAS,4) ELSE CALL ERREUR(487) RETURN ENDIF * IF (ICH7.NE.0) THEN CALL DYNE18(ICH7,KTQ,3,3,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * IF (ICH8.NE.0) THEN CALL DYNE18(ICH8,KTQ,3,4,KCPR) ELSE CALL ERREUR(487) RETURN ENDIF * IF (ICH9.NE.0) THEN CALL DYNE18(ICH9,KTQ,4,1,KCPR) ENDIF IF (ICH10.NE.0) THEN CALL DYNE18(ICH10,KTQ,5,1,KCPR) ENDIF * * IF (NLIAA.NE.0) THEN * * l'indice VARIABLES_LIAISON_A n'existe que pour * les liaisons POLYNOMIALES * et liaison COUPLAGE_DEPLACEMENT avec CONVOLUTION * IPOLY = 0 MTABLE = ITABL SEGACT,MTABLE NIND1 = MLOTAB if(nbesc.ne.0)segact ipiloc DO 100 I=1,NIND1 TYPIND = MTABTI(I) IF (TYPIND.EQ.'MOT ') THEN IP = MTABII(I) ID = IPCHAR(IP) IFI = IPCHAR(IP+1) IL1 = IFI - ID IF (IL1.EQ.ILC1) THEN IF (CHAI1.EQ.ICHARA(ID:IFI-1)) THEN IPOLY = 1 ENDIF ENDIF ENDIF 100 CONTINUE if(nbesc.ne.0)SEGDES,IPILOC SEGDES,MTABLE IF (IIMPI.EQ.333) THEN WRITE(IOIMP,*) 'DEVINI : IPOLY = ',IPOLY ENDIF IF (IPOLY.NE.0) THEN CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_A', & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR) CALL DYNA14(ITREFR,KTLIAA) IF (IERR.NE.0) RETURN ENDIF ENDIF IF (NLIAB.NE.0) THEN *old CALL DEVRCO(Q1,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB, CALL DEVRCO(Q1,Q2,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB, & IBASB,IPLSB,INMSB,IORSB,3,IAROTA) * --> XPTB(:,3,:)=x_0 XPTB(:,4,:)=\dot{q}_0 * ok car XPTB(:,4,:) utilisé uniquement pour calcul FnlB CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_B', & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR) IF (IERR.NE.0) RETURN CALL DYNE14(ITREFR,KTLIAB) IF (IERR.NE.0) RETURN ENDIF * ELSE * * Chpoints initiaux de deplacement et de vitesse: --------------- * IF (ITINIT.NE.0) THEN if (lmodyn) then mwinit = itinit segact mwinit idepl = jpdep ivite = jpvit else TYPRET=' ' CALL ACCTAB(ITINIT,'MOT',I0,X0,'DEPLACEMENT',L0,IP0, & TYPRET,I1,X1,CHARRE,L1,IDEPL) * TYPRET=' ' CALL ACCTAB(ITINIT,'MOT',I0,X0,'VITESSE',L0,IP0, & TYPRET,I1,X1,CHARRE,L1,IVITE) endif ENDIF * * Mise des valeurs initiales au pas m (indice 3) c Q1(3)=q_0 c Q2(3)=\dot{q}_0 IF (IDEPL.NE.0) THEN CALL DYNE18(IDEPL,KTQ,1,3,KCPR) MTQ = KTQ ENDIF IVINIT = 0 IF (IVITE.NE.0) THEN CALL DYNE18(IVITE,KTQ,2,3,KCPR) IVINIT = 1 ENDIF * * Ajout des forces de liaison a t=0 * PDT=XDT(1) T=XTEMPS(1) IF (NLIAA.NE.0) THEN CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA, & NLIAA,PDT,T,1,3,FINERT,IVINIT,FTEST) ENDIF IF (NLIAB.NE.0) THEN IF (RIGIDE) THEN DO 13 IP=1,NPLB DO 15 ID=1,IDIM FEXB(IP,1,ID) = FTEXB(IP,1,1,ID) 15 CONTINUE 13 CONTINUE ENDIF c SEGMENT,MTPHI c INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB) c INTEGER IAROTA(NSB) c REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB) c ENDSEGMENT CALL DEVLF2(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB, & NLIAB,XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB, & PDT,T,1,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,3, & FEXPSM,NPC1,IERRD,FTEST2,XABSCI,XORDON, & NIP,IAROTA,RIGIDE,FEXB,XCHPFB) IF (IERRD.NE.0) THEN IF (IERRD.EQ.1) CALL ERREUR(528) RETURN ENDIF ENDIF c fin distinction reprise/CI ----------------------------------- IF (NLIAA.NE.0 .OR. NLIAB.NE.0) THEN DO 5 I=1,NA1 FINERT(I,4) = FINERT(I,3) 5 CONTINUE * end do ENDIF * * Calcul des deplacements et des vitesses au pas -1/2 (indice 4) * c PDT = XDT(1) PDTS2 = PDT / 2.D0 AUX2 = PDT * PDT / 8.D0 IF (IIMPI.EQ.333) THEN DO 91 I=1,NA1 WRITE(IOIMP,*)'DEVINI :' WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3) WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4) WRITE(IOIMP,*)'DEVINI : Fexa (',I,',1) =',Fexa(I,1) WRITE(IOIMP,*)'DEVINI : Fexa (',I,',2) =',Fexa(I,2) 91 CONTINUE * end do ENDIF c Cas selon nature (pleine ou diag) des matrices K,C et M : c------- Cas K,C ou M pleine ----------------------------------- IF(NB1K.GT.1.OR.NB1C.GT.1.OR.NB1M.GT.1) THEN IF (IIMPI.EQ.333) WRITE(IOIMP,*)'pleine K,C,M',NB1K,NB1C,NB1M c F_0 = F^liaison_0 + F^ext_0 - K * Q_0 DO 61 I = 1,NA1 FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1) 61 CONTINUE CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,3) c F^AMOR_0 = C * \dot Q_0 SEGINI,MAMOR CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,3) c \ddot Q_0 = M-1 * (F_0 - F^AMOR_0) CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,3) c Q_-1/2 = Q_0 - dt/2 * \dot Q_0 + dt^2/8 * \ddot Q_0 DO 62 I = 1,NA1 Q1(I,4) = Q1(I,3) - PDTS2*Q2(I,3) + AUX2*Q3(I,3) 62 CONTINUE c F_-1/2 = F^liaison_-1/2 + F^ext_-1/2 - K * Q_-1/2 DO 63 I = 1,NA1 FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2) 63 CONTINUE CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,4) c \dot Q_-1/2 = [4I+dt*M-1*C]-1 c * ( 4 \dot Q_0 - dt*M-1*(F_-1/2+F_0-F^AMOR_0) ) DO 64 I = 1,NA1 FTEMP(I) = FTOTA(I,4) + FTOTA(I,3) - FAMOR(I,3) 64 CONTINUE CALL DEVLM1(QTEMP,XM,XOPER,FTEMP,NA1,NB1,3,NB1M) DO 65 I = 1,NA1 FTEMP(I) = 4.D0*Q2(I,3) - PDT*QTEMP(I) 65 CONTINUE IF(NB1C.GT.1.OR.NB1M.GT.1) THEN * initialisation de l'operateur (local) et de son inverse * XOP = 4I - dt*M^-1*C SEGINI,MOP c ...C pleine, M pleine IF(NB1C.GT.1.AND.NB1M.GT.1) THEN DO 662 IA=1,NA1 XOP(IA,IA) = 4.D0 DO 663 IB=1,NA1 AUX = 0.D0 DO 664 IC=1,NA1 AUX = AUX + XOPER(IA,IC,3)*XASM(IC,IB) 664 CONTINUE XOP(IA,IB) = XOP(IA,IB) - PDT*AUX 663 CONTINUE 662 CONTINUE c ...C pleine, M diago ELSEIF(NB1C.GT.1) THEN DO 665 IA=1,NA1 XOP(IA,IA) = 4.D0 DO 666 IB=1,NA1 XOP(IA,IB) = XOP(IA,IB) - PDT*XASM(IA,IB)/XM(IA,1) 666 CONTINUE 665 CONTINUE c ...C diago, M pleine ELSEIF(NB1M.GT.1) THEN DO 667 IA=1,NA1 XOP(IA,IA) = 4.D0 DO 668 IB=1,NA1 XOP(IA,IB) = XOP(IA,IB) - PDT*XOPM1(IA,IB)*XASM(IB,1) 668 CONTINUE 667 CONTINUE ENDIF CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET) IF (IIMPI.EQ.333) THEN WRITE(IOIMP,*) 'FTEMP(:)=',(FTEMP(jou),jou=1,NA1) do iou=1,NB1 WRITE(IOIMP,*) 'XOPM1(',iou,',:)=', & (XOPM1(iou,jou),jou=1,NB1) enddo ENDIF c \dot Q_-1/2 = [4I-dt*M-1*C]-1 * FTEMP DO 67 IA = 1,NA1 Q2(IA,4) = 0.D0 DO 68 IB = 1,NA1 Q2(IA,4) = Q2(IA,4) + XOPM1(IA,IB)*FTEMP(IB) 68 CONTINUE 67 CONTINUE SEGSUP,MOP ELSE DO 69 IA = 1,NA1 Q2(IA,4) = FTEMP(IA) / (4.D0 - (PDT*XASM(IA,1)/XM(IA,1))) 69 CONTINUE ENDIF c \ddot Q_-1/2 = M-1 * ( F_-1/2 - F^AMOR_-1/2) CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,4) CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,4) SEGSUP,MAMOR c------- Cas K,C et M diagonal ----------------------------------- ELSE DO 10 I = 1,NA1 UNSM3 = 1.D0 / ( XM(I,1) - FINERT(I,3) ) UNSM4 = 1.D0 / ( XM(I,1) - FINERT(I,4) ) * * G = F - K Q * i,0 i,0 i i,0 * FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1) - (XK(I,1)*Q1(I,3)) AUX1 = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3)) AUX3 = Q1(I,3) - ( PDTS2 * Q2(I,3) ) * * . 2 . * Q = Q - h/2 Q + h /8 ( G - A Q ) * i,-1/2 i,0 i,0 i,0 i i,0 * Q1(I,4) = AUX3 + ( AUX2 * AUX1 ) * * G = F - K Q * i,-1/2 i,-1/2 i i,-1/2 * FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2) - (XK(I,1)*Q1(I,4)) AUX4 = 4.D0 + ( XASM(I,1) * PDT * UNSM3 ) AUX5 = ( FTOTA(I,4) * UNSM4 ) + ( FTOTA(I,3) * UNSM3 ) AUX6 = 1.D0 / ( 4.D0 - ( XASM(I,1) * PDT * UNSM3 ) ) * * . 1 . * Q = ________ ( ( 4 + A h ) Q - h ( G + G ) ) * i,-1/2 4 - A h i i,0 i,-1/2 i,0 * i * Q2(I,4) = AUX6 * ( (AUX4*Q2(I,3)) - (PDT*AUX5) ) * * " . * Q = ( G - A Q ) / M * i i i i i * Q3(I,3) = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3)) Q3(I,4) = (FTOTA(I,4)*UNSM4) - (XASM(I,1)*UNSM4*Q2(I,4)) 10 CONTINUE ENDIF * ENDIF IF (IIMPI.EQ.333) THEN DO 30 I=1,NA1 WRITE(IOIMP,*)'DEVINI :' WRITE(IOIMP,*)'DEVINI : Q1(',I,',3) =',Q1(I,3) WRITE(IOIMP,*)'DEVINI : Q2(',I,',3) =',Q2(I,3) WRITE(IOIMP,*)'DEVINI : Q3(',I,',3) =',Q3(I,3) WRITE(IOIMP,*)'DEVINI :' WRITE(IOIMP,*)'DEVINI : Q1(',I,',4) =',Q1(I,4) WRITE(IOIMP,*)'DEVINI : Q2(',I,',4) =',Q2(I,4) WRITE(IOIMP,*)'DEVINI : Q3(',I,',4) =',Q3(I,4) 30 CONTINUE * end do ENDIF IF (IIMPI.EQ.333) THEN DO 31 I=1,NA1 WRITE(IOIMP,*)'DEVINI :' WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3) WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4) 31 CONTINUE * end do ENDIF * ICPR = KCPR SEGSUP,ICPR * END