C KMDMT SOURCE CB215821 20/11/25 13:31:31 10792 SUBROUTINE KMDMT(MTABP,MCHB,MCHR,IMPR,CBETA,KDPDQ,KPIMP,PIMP,NIMP) C *********************************************************************** C -1 T C CE SOUS PROGRAMME CALCULE LA MATRICE C D C CONNAISSANT C C C POINTEURS : C C IZL CONTIENT LE PROFIL ET LA RENUMEROTATION (CF PROFCH) C MELEME OBJET MAILLAGE SUR LEQUEL EST FAIT LE CALCUL C MATRAK MATRICES ELEMENTAIRES DE LA DIVERGENCE (ALIAS"C") C C IZDV DIAGONALE D C IZIPAD CORRESPONDANCE NUMER. GLOBALE --> NUMER. LOCALE C (DOMAINE SUR LEQUEL PORTE AP ET NON LA PRESSION) C C EN SORTIE : C C LA DEMI-MATRICE INFERIEUR EST RANGEE PAR LIGNE C (UN SEGMENT IZA PAR LIGNE) C LA DIAGONALE DE LA MATRICE EST RANGEE A PART DANS IZD C KZA CONTIENT TOUT LES POINTEURS IZA C IZD=KZA(1) C C *********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMELEME POINTEUR IPTJ.MELEME,IPTK.MELEME,MELEM1.MELEME POINTEUR IGEOM.MELEME,IGEOMI.MELEME,IGEOMC.MELEME POINTEUR MELEMK.MELEME,MELEMC.MELEME POINTEUR MELSTB.MELEME C? POINTEUR ELTFA.MELEME,FACEL.MELEME,MELEMC.MELEME C? POINTEUR MELEMF.MELEME,MELEMK.MELEME,MELELI.MELEME -INC SMLENTI POINTEUR IZIPAD.MLENTI,MLENTP.MLENTI,IZIPAK.MLENTI -INC SMCOORD -INC SMCHPOI POINTEUR IZDV.MCHPOI POINTEUR IZDP.MCHPOI,IZDDP.MPOVAL POINTEUR MCHB.MCHPOI,MCHR.MCHPOI C? POINTEUR IZB.MPOVAL,IZR.MPOVAL,IZPOF.MPOVAL,IZPOC.MPOVAL POINTEUR IZB.MPOVAL,IZR.MPOVAL,IZPOC.MPOVAL C? POINTEUR IZP.MPOVAL,IZF.MPOVAL POINTEUR IZP.MPOVAL C? POINTEUR IZFAC.MPOVAL,IZVOL.MPOVAL C-INC SMMATRAKANC -INC SMMATRK1 C************************************************************************* C C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees C * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange * (points CENTRE ) pour chaque operateur de contrainte * KGEOC SPG pour la totalite des points CENTRE. * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse) * KLEMC Connectivites de l'ensemble des contraintes * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones SEGMENT MATRAK INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP) INTEGER LIZAFM(NBSOUS) INTEGER IKAM0 (NBSOUS) INTEGER IMEM (NBELC) INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC ENDSEGMENT SEGMENT IZAFM REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX) ENDSEGMENT POINTEUR IPMJ.IZAFM,IPMK.IZAFM C******************************************************************* SEGMENT/IZA/(A(LA)*D) SEGMENT/IZD/(D(NELP)*D) -INC SMTABLE POINTEUR MTABP.MTABLE C CHARACTER*8 TYPE,TYPC * SAVE BETA C*** C KDPDQ=0 C write(6,*)' Entree dans kmdmt KDPDQ=',kdpdq TYPE=' ' CALL ACMO(MTABP,'MATC',TYPE,MATRAK) C WRITE(6,*)' KMDMT : MATC=',MATRAK IF(TYPE.NE.'MATRAK')THEN WRITE(6,*)' IL N''Y A PAS D''ENTREE MATC DANS LA TABLE ' GO TO 90 ENDIF SEGACT MATRAK * initialisations pour plaire a l'optimiseur IBM MELEMK=MATRAK IZPOC=MATRAK IZIPAK=MATRAK MELSTB=MATRAK IZA=MATRAK IZL=MATRAK MLENTI=MATRAK IZD=MATRAK IDBLK=MATRAK IPMJ=MATRAK NBELC=IMEM(/1) NIMP=NBELC/2 IF(KPIMP.EQ.2)NIMP=NBELC C write(6,*)' CAS KPIMP=',kpimp,' PIMP=',pimp CALL PROFCH(MATRAK,IMPR,KPIMP,NIMP) TYPE=' ' CALL ACMO(MTABP,'DIAGV',TYPE,IZDV) IF(TYPE.NE.'CHPOINT ')THEN WRITE(6,*)' IL N''Y A PAS D''ENTREE DIAGV DANS LA TABLE ' GO TO 90 ENDIF SEGACT MATRAK*MOD IZL=KIZCL C write(6,*)' KIZCL=',KIZCL C WRITE(6,*)' KDIAG=',KDIAG,' IZDV=',IZDV,' KCAC=',KCAC IF(KDIAG.EQ.IZDV.AND.KCAC.EQ.1)THEN GO TO 500 ENDIF IF(KDPDQ.NE.0)THEN C DPDQ TYPE=' ' CALL ACMO(MTABP,'DOMAINE',TYPE,MTABD) IF(TYPE.NE.'TABLE')GO TO 90 TYPE=' ' CALL ACMO(MTABD,'MELSTB',TYPE,MELSTB) IF(TYPE.NE.'MAILLAGE') GO TO 90 CALL LEKTAB(MTABD,'CENTRE',MELEMC) TYPE=' ' CALL ACMO(MTABD,'MCHPOC',TYPE,MCHPOC) IF(TYPE.NE.'CHPOINT') GO TO 90 CALL LICHT(MCHPOC,IZPOC,TYPC,IGEOMC) NBK=IZPOC.VPOCHA(/2) SEGACT MELEMC,IGEOMC C DPDQ ENDIF KDIAG=IZDV CALL LICHT(IZDV,MPOVAL,TYPC,IGEOM) C WRITE(6,*)' DIAGONALE : ' C LDG=VPOCHA(/1) C WRITE(6,*)' SUR X ' C WRITE(6,1002)(VPOCHA(II,1),II=1,LDG) C WRITE(6,*)' SUR Y ' C WRITE(6,1002)(VPOCHA(II,2),II=1,LDG) C C ON COMMENCE PAR ACTIVER C SEGACT IZL C. NJAN=NUAN(/1) C. write(6,*)' NUAN=' C. write(6,1001)(NUAN(II),II=1,NJAN) C. write(6,*)' NUNA=' C. write(6,1001)(NUNA(II),II=1,NJAN) C write(6,*)' IMEL=' C write(6,1001)(IMEL(II),II=1,NJAN) C write(6,*)' IMJ=' C write(6,1001)(IMJ(II),II=1,NJAN) MELEM1=KGEOS CALL KRIPAD(MELEM1,IZIPAD) MELEME=KLEMC SEGACT MELEME IKAS=0 IF(IKAS.EQ.1)THEN NIZ=LIZAFM(/1) DO 1940 L=1,NIZ IZAFM=LIZAFM(L) WRITE(6,*)' MATRICE AM BLOC N° ',L SEGACT IZAFM NP=AM(/2) IES=AM(/3) NNELP=AM(/1) WRITE(6,*)'NP,IES,NNELP=',NP,IES,NNELP DO 1941 K=1,NNELP WRITE(6,*)'K=',K WRITE(6,1008)((AM(K,I,J),I=1,NP),J=1,IES) 1941 CONTINUE 1940 CONTINUE C CALL ARRET(0) 1008 FORMAT(8(1X,1PE11.4)) ENDIF IES=VPOCHA(/2) NPT=VPOCHA(/1) NELP=B(/1) C IDMAT=KZA1 SEGACT IDMAT*MOD NBLK=IDESCR(/1) WRITE(6,*) '* NOMBRE DE BLOCS DE LA MATRICE = ',NBLK SEGINI IZD IDIAG=IZD C DO 71 NS=1,LIZAFM(/1) IZAFM=LIZAFM(NS) SEGACT IZAFM C C NP1=AM(/1) C NP2=AM(/3) C WRITE(6,*)' NP=',NP1,' NELB=',NP2 C DO 1834 I=1,NP1 C WRITE(6,*)' I1=',I C WRITE(6,1002)(AM(I,1,KK),KK=1,NP2) C WRITE(6,*)' I2=',I C WRITE(6,1002)(AM(I,2,KK),KK=1,NP2) C1834 CONTINUE 71 CONTINUE NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 DO 10 JIS=1,NBSOUS IF(NBSOUS.EQ.1)IPTJ=MELEME IF(NBSOUS.NE.1)IPTJ=LISOUS(JIS) SEGACT IPTJ 10 CONTINUE IF(KDPDQ.NE.0)THEN C C OPERATEUR DPDQ C CALL KRIPAD(IGEOMC,MLENTI) SEGACT MELSTB,IGEOMC MELEMK=KGEOC CALL KRIPAD(MELEMK,IZIPAK) SEGACT MELEMK ENDIF C Calcul de la diagonale AD=0.D0 ADB=1.D0 NBELP=NUNA(/1) DO 7710 KJ=1,NBELP KJJ=NUNA(KJ) IPMJ=LIZAFM(IMEM(NUNA(KJ))) JKAM=IKAM0 (IMEM(NUNA(KJ))) KJJ=KJJ-JKAM+1 J=IMJ(KJ) IPTJ=IMEL(KJ) NPJ=IPTJ.NUM(/1) K1=KJ-KZA(KJ)+1 KK=KJ KKK=NUNA(KK) IPMK=LIZAFM(IMEM(NUNA(KK))) KKAM=IKAM0 (IMEM(NUNA(KK))) KKK=KKK-KKAM+1 K=IMJ(KK) IPTK=IMEL(KK) NPK=IPTK.NUM(/1) IF(KK.NE.KJ) KA=KA+1 AUX=0.D0 DO 7731 I=1,NPJ DO 7732 N=1,NPK IF(IPTJ.NUM(I,J).NE.IPTK.NUM(N,K))GO TO 7732 JD=IZIPAD.LECT(IPTK.NUM(N,K)) AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) IF(IES.EQ.3) &AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 7732 CONTINUE 7731 CONTINUE C AD=AD+AUX D(KJ)=D(KJ)+AUX 7710 CONTINUE AD=AD/FLOAT(NBELP) IF(KDPDQ.NE.0)THEN C C OPERATEUR DPDQ C NBELP=NUNA(/1) ADB=0.D0 DO 7711 KJ=1,NBELP KC=NUNA(KJ) LPC=LECT(MELEMK.NUM(1,KC)) IF(LPC.NE.0)THEN ADB=ADB+ABS(IZPOC.VPOCHA(LPC,1)) ENDIF 7711 CONTINUE ADB=ADB/FLOAT(NBELP) BETA0=AD/ADB*5.E-4 BETA=CBETA*BETA0 C WRITE(6,7271)AD,ADB,BETA0,CBETA C7271 FORMAT(1X,' AD=',1PE11.4,' ADB=',1PE11.4,' BETAopt=',1PE11.4, WRITE(6,7271)BETA0,CBETA 7271 FORMAT(1X,' BETAopt=',1PE11.4,' CBETA=',1PE11.4) KJ=1 KC=NUNA(KJ) LPC=LECT(MELEMK.NUM(1,KC)) IF(LPC.NE.0)THEN P1=BETA*IZPOC.VPOCHA(LPC,1) D(KJ)=D(KJ)+ABS(P1) ENDIF ENDIF C FIN DPDQ C C CALCUL DES BLOCS C DO 100 IBLK=1,NBLK C NUMÉROS DES LIGNES DE DÉBUT ET DE FIN DE BLOC KJD,KJF KJD=NLDBLK(IBLK) KJF=NLDBLK(IBLK+1)-1 C ALLOCATION DE LA MÉMOIRE POUR LE BLOC IDBLK=IDESCR(IBLK) SEGACT IDBLK*MOD LA=ILON WRITE(6,*) '* BLOC NUMÉRO : ',IBLK,' TAILLE = ',ILON,' MOTS' SEGINI IZA IMAT=IZA KA=0 DO 1 KJ=KJD,KJF KJJ=NUNA(KJ) IPMJ=LIZAFM(IMEM(NUNA(KJ))) JKAM=IKAM0 (IMEM(NUNA(KJ))) KJJ=KJJ-JKAM+1 J=IMJ(KJ) IPTJ=IMEL(KJ) NPJ=IPTJ.NUM(/1) K1=KJ-KZA(KJ)+1 DO 2 KK=K1,KJ KKK=NUNA(KK) IPMK=LIZAFM(IMEM(NUNA(KK))) KKAM=IKAM0 (IMEM(NUNA(KK))) KKK=KKK-KKAM+1 K=IMJ(KK) IPTK=IMEL(KK) NPK=IPTK.NUM(/1) IF(KK.NE.KJ) KA=KA+1 AUX=0.D0 DO 31 I=1,NPJ DO 32 N=1,NPK IF(IPTJ.NUM(I,J).NE.IPTK.NUM(N,K))GO TO 32 JD=IZIPAD.LECT(IPTK.NUM(N,K)) AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1) & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2) IF(IES.EQ.3) &AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3) 32 CONTINUE 31 CONTINUE C IF(KK.EQ.KJ)GO TO 21 A(KA)=A(KA)+AUX C 21 CONTINUE 2 CONTINUE C?. ad=ad+aux C?. D(KJ)=D(KJ)+AUX IF(KDPDQ.NE.0)THEN C OPERATEUR DPDQ C KJJ=NUNA(KJ) LPJ=LECT(MELEMK.NUM(1,KJJ)) IF(LPJ.NE.0)THEN P1=BETA*IZPOC.VPOCHA(LPJ,1) D(KJ)=D(KJ)+ABS(P1) DO 3819 JC=2,NBK NPC=MELSTB.NUM(JC,LPJ) P=BETA*IZPOC.VPOCHA(LPJ,JC) KPI=IZIPAK.LECT(NPC) KPJN=NUAN(KPI) IF(KPJN.GT.KJ)GO TO 3819 IDECI=IDEBLK(KJ-KJD+1) LLI=IDEBLK(KJ-KJD+2)-IDEBLK(KJ-KJD+1) KJAA=LLI-KJ+KPJN+IDECI A(KJAA)=A(KJAA)+P 3819 CONTINUE ENDIF ENDIF C FIN DPDQ 1 CONTINUE SEGDES IZA SEGDES IDBLK 100 CONTINUE SEGDES IDMAT IF(KDPDQ.NE.0)THEN SEGSUP IZIPAK,MLENTI ENDIF IF(KPIMP.EQ.1)THEN D(NIMP)=1.D30 ELSEIF(KPIMP.EQ.2)THEN SEGACT IZA NELP=D(/1) LGLA=A(/1) LGLI=KZA(NIMP) CALL INITD(A(LGLA-LGLI+1),LGLI,0.D0) D(NELP)=1.D-10 CALL LEKTAB(MTABD,'XXVOLUM',MCHPOI) CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) VOLT=0.D0 DO 482 I=1,NELP VOLT=VOLT+VPOCHA(I,1) 482 CONTINUE DO 483 I=1,NELP XVOL=VPOCHA(NUNA(I),1) A(LGLA-LGLI+I)=-XVOL/VOLT 483 CONTINUE ENDIF SEGSUP IZIPAD DO 20 JIS=1,NBSOUS IF(NBSOUS.EQ.1)IPTJ=MELEME IF(NBSOUS.NE.1)IPTJ=LISOUS(JIS) C write(6,*)' On desactive JIS,IPTJ=',JIS,IPTJ SEGDES IPTJ IPMJ=LIZAFM(JIS) JKAM=IKAM0 (JIS) SEGDES IPMJ 20 CONTINUE SEGDES MPOVAL,IZDV KCAC=1 TYPE=' ' CALL ACMO(MTABP,'DIAGP',TYPE,IZDP) IF(TYPE.EQ.'CHPOINT ')THEN CALL LICHT(IZDP,IZDDP,TYPE,IGEOMI) IGEOM=KGEOC SEGACT IZD SEGACT IGEOM,IGEOMI NPT=D(/1) NPTI=IZDDP.VPOCHA(/1) DO 51 I=1,NPTI XVALI=IZDDP.VPOCHA(I,1) IPOI=IGEOMI.NUM(1,I) DO 52 J=1,NPT IF(IGEOM.NUM(1,J).EQ.IPOI)THEN JJ=NUAN(J) D(JJ)=D(JJ)+XVALI GO TO 53 ENDIF 52 CONTINUE 53 CONTINUE 51 CONTINUE SEGDES IZDP,IZDDP SEGDES IGEOM,IGEOMI ENDIF DO 72 NS=1,LIZAFM(/1) IZAFM=LIZAFM(NS) SEGDES IZAFM 72 CONTINUE SEGDES IZD SEGDES MATRAK SEGDES MELEME C WRITE(6,*)' APPEL A TRIAKS' CALL TRIAKS(IZL) C WRITE(6,*)' RETOUR DE TRIAKS' 500 CONTINUE CALL LICHT(MCHB,IZB,TYPC,IGEOM) CALL LICHT(MCHR,IZR,TYPC,IGEOM) CALL KAL2P(MTABP,MCHB) C write(6,*)' APPEL A RESOCK' IF(KPIMP.NE.0)THEN SEGACT IZL CALL LICHT(MCHB,IZB,TYPC,IGEOM) IF(KPIMP.EQ.1)IZB.VPOCHA(NUNA(NIMP),1)=PIMP*1.D30 IF(KPIMP.EQ.2)IZB.VPOCHA(NUNA(NIMP),1)=-PIMP ENDIF CALL RESOCK(IZB,IZL,IZR) C write(6,*)' Retour de resock ' SEGDES MCHB,IZB SEGDES MCHR,IZR RETURN 90 CONTINUE WRITE(6,*)' ARRET ANORMAL DE KMDMT' RETURN 1002 FORMAT(10(1X,1PE11.4)) 1001 FORMAT(20(1X,I5)) END