C KSPRUS SOURCE CB215821 20/11/25 13:33:12 10792 SUBROUTINE KSPRUS &(MELEME,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0, & ANUKK,IDCEN,IKOMP,NPTU,IPADU,IZTU1,TETAN,NP,NEL,CMD) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C***************************************************************************** C C Ce SP calcule les matrices elementaires de divergence alias C C C IKAS=1 KMCT calcul de Ct (Div U) C IKAS=2 KMAC calcul de C uniquement (Grad p) C IKAS=3 KCCT calcul de C assemblage pour C et Ct C C***************************************************************************** CHARACTER*8 NOM0 -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMCOORD -INC SIZFFB POINTEUR IZF1.IZFFM -INC SMMATRIK -INC SMLENTI POINTEUR IPADU.MLENTI -INC SMELEME -INC SMCHPOI POINTEUR IZTGG1.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL -INC CCREEL PARAMETER (LRV=64) SEGMENT PETROV REAL*8 WT(LRV,NP,NPG,IDIM),WS(LRV,NP,NPG,IDIM),HK(LRV,IDIM,NP,NPG) REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG) REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV) REAL*8 AJK(LRV,IDIM,IDIM,NPG) ENDSEGMENT DIMENSION KIPM(3),XYZ1(24) C C OPERATEUR PRESSION C DEUPI=1.D0 IF(IAXI.NE.0)DEUPI=2.D0*XPI IF(IKAS.EQ.3)DEUPI=-DEUPI NBEL = NEL C write(6,*)' DEUPI=',deupi IF(IDIM.EQ.2)IPM3=IPM1 KIPM(1)=IPM1 KIPM(2)=IPM2 KIPM(3)=IPM3 SEGACT MELEME SEGACT IPM1*MOD,IPM2*MOD,IPM3*MOD SEGACT IPADU,IZTU1,TETAN,IZTGG1 IF(KPRE.NE.2)THEN C C ON TRAITE LE CAS DIV U -> T sommet C IF(INEFMD.EQ.3)THEN IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'PRP0' IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'PRP1' IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'PFP1' ELSEIF(INEFMD.EQ.2)THEN IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'MCP0' IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'MCP1' IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'MCF1' ELSEIF(INEFMD.EQ.4)THEN NOM0=NOMS(ITYPEL)//' ' ENDIF c write(*,*)' Kpruss : NOM0=',nom0,' ikas=',ikas,' nel=',nel c write(*,*)' IK1',IK1,' COEF',IZTGG1.VPOCHA(1,1),'ANUK',ANUKK CALL KALPBG(NOM0,'FONFORM ',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=GR(/1) NPG=GR(/3) IZF1=KTP(1) SEGACT IZF1*MOD * ON modifie pour avoir les bonnes fonctions de forme. GBM 18/10/99 MP1=FN(/1) * * On fait des packets * NNN=MOD(NBEL,LRV) IF(NNN.EQ.0) THEN NPACK=NBEL/LRV ELSE NPACK=1+(NBEL-NNN)/LRV ENDIF KPACKD=1 KPACKF=NPACK * * boucle sur les packets, on active la structure petrov * SEGINI PETROV DO 7001 KPACK=KPACKD,KPACKF C C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= C C 1. Calcul des limites du paquet courant. KDEB=1+(KPACK-1)*LRV KFIN=MIN(NBEL,KDEB+LRV-1) C DO 7002 K=KDEB,KFIN KP=K-KDEB+1 NK = K + K0 NKR=(1-IK1)*(NK-1)+1 COEFK(KP)=IZTGG1.VPOCHA(NKR,1) + 1.D-30 * WRITE(*,*) IZTGG1.VPOCHA(NKR,1), COEFK(KP) ANUK(KP)=ANUKK/COEFK(KP) 7002 CONTINUE C C Une seule composante de projection car div se projete C sur un scalaire. C c write (*,*) 'IPADU', IPADU.LECT(/1) C CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN DT=0.D0 CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, & NES,NP,NPG,IAXI,XCOOR, & MELEME.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP, & TETAN.VPOCHA,IPADU.LECT,IZTU1.VPOCHA,IPADU.LECT,NPTU,ANUK, & WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL, & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) DO 30 KE=KDEB,KFIN NK=KE-KDEB+1 * DO 12 I=1,NP * J=NUM(I,KE) * DO 12 N=1,IDIM * XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N) * 12 CONTINUE * CALL CALJBR * &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) C write(6,*)' Retour caljbr ',mp1,np,npg DO 324 K=1,IDIM IPM4=KIPM(K) DO M=1,MP1 DO I=1,NP U=0.D0 DO 333 L=1,NPG * ON modifie la formule GBM 18/10/99 U=U+WT(NK,M,L,1)*HK(NK,K,I,L)*PGSK(NK,L)*DEUPI*RPGK(NK,L) * U=U+FN(M,L)*HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L) 333 CONTINUE IF(IAXI.NE.0.AND.K.EQ.1)THEN DO 334 L=1,NPG * modifié par GBM 18/10/99 U=U+WT(NK,M,L,1)*FN(I,L)*PGSK(NK,L)*DEUPI * U=U+FN(M,L)*FN(I,L)*PGSQ(L)*DEUPI 334 CONTINUE ENDIF C C K1=1+(1-IK1)*(NK-1) C attention valable que pour coef au centre U=U*COEFK(NK) * WRITE(*,*) NK,M,I,COEFK(NK) C if(ikas.ne.2)then IPM4.AM(KE,I,M)=IPM4.AM(KE,I,M)+U else IPM4.AM(KE,M,I)=IPM4.AM(KE,M,I)+U endif * write(*,*) KE,PGSQ ENDDO ENDDO 324 CONTINUE 30 CONTINUE 7001 CONTINUE SEGSUP IZHR,IZFFM,PETROV C CAS MACRO CENTRE ELSEIF(KPRE.EQ.2)THEN * initialisation de npg au pif par TC NPG=1 SEGINI PETROV NOM0=NOMS(ITYPEL)//' ' CALL KALPBG(NOM0,'FONFORM ',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=GR(/1) NPG=GR(/3) IZF1=KTP(1) SEGACT IZF1*MOD MPG=IZF1.FN(/2) NP=GR(/2) DO 40 KE=1,NEL NK=K0+KE IX=0 DO I=1,NP J=MELEME.NUM(I,KE) DO N=1,IDIM IX=IX+1 XYZ1(IX)=XCOOR((J-1)*(IDIM+1) +N) ENDDO ENDDO CALL CALJBR(FN,GR,PG,XYZ1,HR,PGSQ,RPG,NES, & IDIM,NP,NPG,IAXI,AIRE,AJ,SGN) DO 424 K=1,IDIM IPM4=KIPM(K) DO 423 I=1,NP U=0.D0 DO 433 L=1,NPG U=U+HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L) 433 CONTINUE IF(IAXI.NE.0.AND.K.EQ.1)THEN DO 434 L=1,NPG U=U+FN(I,L)*PGSQ(L)*DEUPI 434 CONTINUE ENDIF K1=1+(1-IK1)*(NK-1) U=U*IZTGG1.VPOCHA(K1,1) if(ikas.ne.2)then IPM4.AM(KE,I,1)=IPM4.AM(KE,I,1)+U else IPM4.AM(KE,1,I)=IPM4.AM(KE,1,I)+U endif 423 CONTINUE 424 CONTINUE 40 CONTINUE SEGSUP IZHR,IZFFM,PETROV ENDIF RETURN 1002 FORMAT(10(1X,1PE11.4)) 1040 FORMAT(1X,'CALCUL MATRICE AM ',I4/10(1X,1PE11.4)) END