C VFSYM4 SOURCE LEPOTIER 09/10/12 21:15:15 6527 C NORV4 SOURCE PV 09/03/12 21:29:43 6325 SUBROUTINE VFSYM4(MELEFA,MELEFL,MLECEN,MELEFP,MLESOM,MPONOR, & MPOSUR,MELTFA,MLEFA,MPOTEN,MPOCHP,MLENCL, & MPOVCL,ICHTE,ICHCL,ICHCO,MPOVCO,IOP, & IPO2,SCMB,INDLI,VAL1,VAL2,IND22,IND2,IND, & IPO3,VAUX,TAB,MELEME,MPOGRA,MELVA1,MELVA2, & NBNN,NBFAC,MCHEL2,MCHAM2) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : NORV4 C C DESCRIPTION : Appelle par NORV1 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) C C AUTEUR : C. LE POTIER, DM2S/SFME/MTMS C C************************************************************************ C IMPLICIT INTEGER(I-N) IMPLICIT real*8 (a-h,o-z) -INC SMLENTI -INC SMELEME -INC SMCHPOI -INC CCOPTIO -INC SMCOORD -INC SMLREEL POINTEUR MELEFL.MELEME, MELEFP.MELEME, MELEFA.MELEME, & MELTFA.MELEME POINTEUR MPOSUR.MPOVAL, MPONOR.MPOVAL, & MPOCHP.MPOVAL, MPOVCL.MPOVAL, MPGSOM.MPOVAL, MPVOSO.MPOVAL, & MPOGRA.MPOVAL,MPOTEN.MPOVAL,MPOVCO.MPOVAL POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLESOM.MLENTI, & MLEFA.MLENTI -INC SMCHAML INTEGER NBNN,NBREF C C**** Variables de COOPTIO C C C**** Variables de COOPTIO C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ & ,IOPER, IOSGB, IOGRA, IOSAU, IORES & ,IECHO, IIMPI, IOSPI & ,IDIM C & ,MCOORD & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU & ,NORINC,NORVAL,NORIND,NORVAD & ,NUCROU, IPSAUV C**** Variable de SMLENTI, SMCHPOI C INTEGER JG, N, NC, NSOUPO, NAT, NBSOUS, NBNO,NBELEM C C**** Les includes C INTEGER I1,ICOMP,ICOMGR,IGEOM & ,IOP1,ICEN,ISOMM,IFAC,IFACEL,IFACEP,INORM & ,ISURF,IMAIL,ICHPO,ICHCL,ICHGRA,ICOEFF & ,NTOT,NSOMM,NCOMP,NFAC,NCEN & ,NLCF,NGCF,NGCF1,NGCF2,NGCG,NGCD,NLCG,NLCD,NGS1,NGS2 & ,NLS1,NLS2,NLFCL & ,ISOUS,IELEM,INOEUD,ICELL,NCON,NFIN,INDGA,INDDR INTEGER ICEN2 REAL*8 SCNX,SCNY,SURF,VOL,VAL,VALX,VALY,XG,XD,XF,XS1,XS2 & ,YG,YD,YF,YS1,YS2,PSCA,XNORM,VECX,VECY,PSCAGX,PSCAGY, & PSCADX,PSCADY,K11G,K22G,K21G,K11D,K22D,K21D,VXG1,VXG2, & VXAU,VYAU,VXD1,VXD2,VYG1,VYG2,TRG1,TRG2, & TRD1,TRD2,TRG,TRD,AUX,AUY,UN,TRGNS1,TRGNS2,EXPR1,EXPR2,COEFD REAL*8 XLONG,AG1,AG2,AD1,AD2,PSCAG1,PSCAG2,PSCAD1,PSCAD2, & COEF1,COEF2,COEF3,COEF4,SCN1X,SCN1Y,VX,VY,COEF1X,COEF2X, & COEF1Y,COEF2Y,CX,CY,ANCX,ANCY,DIFFX,DIFFY,XLONGG,XLONGD & VALD,VALG,COEF,GX,GY,XMINK11,XMAXK11,XMINK22,XMAXK22,THETA, & AUXRE,AUXMA REAL*8 VECXG1(2),VECYG1(2) REAL*8 VECXG2(2),VECYG2(2) REAL*8 VECXD1(2),VECYD1(2) REAL*8 VECXD2(2),VECYD2(2) REAL*8 EPS,TRACE INTEGER ICRIT CHARACTER*(4) NOMCOM(18) CHARACTER*8 TYPE C DATA NOMCOM /'P1DX','P1DY', & 'P2DX','P2DY', & 'P3DX','P3DY', & 'P4DX','P4DY', & 'P5DX','P5DY', & 'P6DX','P6DY', & 'P7DX','P7DY', & 'P8DX','P8DY', & 'P9DX','P9DY'/ INTEGER NDIM SEGMENT MMAT1 REAL*8 PM(NDIM,NDIM),PM1(NDIM,NDIM),XSOL(NDIM) INTEGER IC(NDIM) ENDSEGMENT INTEGER K1,K2 SEGMENT INDICE INTEGER NUME(K1,K2) ENDSEGMENT POINTEUR IND.INDICE,IND2.INDICE,IND22.INDICE SEGMENT MATRICE REAL*8 MAT(K1,K2) ENDSEGMENT POINTEUR VAL1.MATRICE,VAL2.MATRICE,SCMB.MATRICE INTEGER K3 SEGMENT POINT2 INTEGER POINT(K3) ENDSEGMENT POINTEUR IPO2.POINT2 SEGMENT MATRICE2 REAL*8 MAT2(K1,K2) ENDSEGMENT POINTEUR MATR1.MATRICE2,MATR2.MATRICE2 SEGMENT POINT3 INTEGER POINT33(K3) ENDSEGMENT POINTEUR IPO3.POINT3 SEGMENT INDICE3 INTEGER NU(K1,K2) ENDSEGMENT POINTEUR INDIC.INDICE3,INDIC2.INDICE3 SEGMENT REP INTEGER ID(K3) ENDSEGMENT POINTEUR TAB.REP,INDLI.REP INTEGER K5 SEGMENT NBFAC INTEGER NBFACEL(K5) INTEGER IMELEM(K5) ENDSEGMENT INTEGER K7,K8 SEGMENT MISZERO INTEGER TABL(K7) INTEGER TABL2(K7) INTEGER TABL1(K8),IPOS(K8),ICOURANT(K8) REAL*8 XMAX(K7) ENDSEGMENT POINTEUR ITAB.MISZERO NMAI1 = 0 NBSO = MAX(1,MELTFA.LISOUS(/1)) IELTFA = MELTFA IF (NBSO.EQ.1) THEN K5 = MELTFA.NUM(/2) ELSEIF (NBSO.EQ.2) THEN IPT1 = MELTFA.LISOUS(1) SEGACT IPT1 N1 = IPT1.NUM(/2) NMAI1 = N1 SEGDES IPT1 IPT2 = MELTFA.LISOUS(2) SEGACT IPT2 N2 = IPT2.NUM(/2) NMAI2 = N2 SEGDES IPT2 K5 = N1 + N2 ENDIF IF (NBSO.EQ.1) THEN DO I = 1,K5 NTYPE = MELTFA.ITYPEL IF (NTYPE .EQ. 4) THEN NBFACEL(I) = 3 IMELEM(I) = MELTFA ELSE NBFACEL(I) = 4 IMELEM(I) = MELTFA ENDIF c SEGDES MELTFA ENDDO ELSEIF (NBSO.EQ.2) THEN IPT1 = MELTFA.LISOUS(1) SEGACT IPT1 IPT2 = MELTFA.LISOUS(2) SEGACT IPT2 DO I = 1,K5 N1 = IPT1.NUM(/2) IF (I.LE.N1) THEN IF (IPT1.ITYPEL .EQ. 4) THEN NBFACEL(I) = 3 IMELEM(I) = IPT1 ELSE NBFACEL(I) = 4 IMELEM(I) = IPT1 ENDIF c SEGDES IPT1 ELSE IF (IPT2.ITYPEL .EQ. 4) THEN NBFACEL(I) = 3 IMELEM(I) = IPT2 ELSE NBFACEL(I) = 4 IMELEM(I) = IPT2 ENDIF c SEGDES IPT2 ENDIF ENDDO ENDIF NFAC=MELEFL.NUM(/2) C EVALUATION DE NBNN NCON = 0 DO NLCF= 1, NFAC, 1 NGCF=MELEFL.NUM(2,NLCF) NGCF1=MELEFA.NUM(1,NLCF) NGCF2=MELEFP.NUM(3,NLCF) IF((NGCF.NE.NGCF1) .OR. (NGCF.NE.NGCF2))THEN WRITE(IOIMP,*) & 'Il ne faut pas jouer avec la table domaine!' CALL ERREUR(5) GOTO 9999 ENDIF NGCG=MELEFL.NUM(1,NLCF) NGCD=MELEFL.NUM(3,NLCF) NLCG=MLECEN.LECT(NGCG) NLCD=MLECEN.LECT(NGCD) NGS1=MELEFP.NUM(1,NLCF) NGS2=MELEFP.NUM(2,NLCF) NLS1=MLESOM.LECT(NGS1) NLS2=MLESOM.LECT(NGS2) DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.NGCF) THEN IAFF = I5 IG = IAFF GOTO 6399 ENDIF ENDDO 6399 CONTINUE INDIC = IPO3.POINT33(NLS1) SEGACT INDIC *MOD INDIC2 = IPO3.POINT33(NLS2) SEGACT INDIC2 *MOD ICON = TAB.ID(NLS1) DO JAUX = 1,TAB.ID(NLS2) J1 = INDIC2.NU(IG,JAUX) c WRITE(6,*) 'IG= ',IG,'JAUX= ',JAUX,'NLS2= ',NLS2,'J1= ',J1 DO IAUX2 = 1,TAB.ID(NLS1) J2 = INDIC.NU(IG,IAUX2) c WRITE(6,*) 'IG= ',IG,'IAUX2= ',IAUX2,'NLS1= ',NLS1,'J2= ',J2 IF ((J1.EQ.J2).OR.(J1.EQ.0).OR.(J2.EQ.0)) THEN GOTO 5999 ENDIF ENDDO ICON = ICON +1 5999 CONTINUE ENDDO SEGDES INDIC *MOD SEGDES INDIC2 *MOD NCON = MAX(NCON,ICON) c WRITE(6,*) 'NCON= ',NCON ENDDO NCON = NCON + 1 NBNN = NCON C DEFINITION DES PARAMETRES DU CHAMELEM DES COEFFICIENTS c INITIALISATION DU CHAMELEM N1=1 N2=1 N3=6 L1=8 SEGINI MCHELM ICOEFF = MCHELM MCHELM.TITCHE='Gradient' MCHELM.IFOCHE=IFOUR C ISOUS=0 NBSOUS=0 NBREF=0 NBELEM = NFAC ISOUS=ISOUS+1 SEGINI MELEME C ITYPEL=32 -> 'POLY' ITYPEL=32 MCHELM.IMACHE(ISOUS)=MELEME SEGINI MCHAML MCHELM.ICHAML(ISOUS)=MCHAML MCHAML.NOMCHE(1)='SCAL' MCHAML.TYPCHE(1)='REAL*8 ' N1PTEL=NCON N1EL=NBELEM N2PTEL=0 N2EL=0 SEGINI MELVA1 MCHAML.IELVAL(1)=MELVA1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCUL DE LA VITESSE EN CHAQUE FACE AA = 0.0 BB = 0.0 C INDICE QUI COMPTE LES COEFFICIENTS POUR CHAQUE FACE NAUX2 = 0 NMOY = 0 DO NLCF= 1, NFAC, 1 NGCF=MELEFL.NUM(2,NLCF) NGCF1=MELEFA.NUM(1,NLCF) NGCF2=MELEFP.NUM(3,NLCF) IF((NGCF.NE.NGCF1) .OR. (NGCF.NE.NGCF2))THEN WRITE(IOIMP,*) & 'Il ne faut pas jouer avec la table domaine!' CALL ERREUR(5) GOTO 9999 ENDIF MELEME.NUM(1,NLCF)=NGCF MELVA1.VELCHE(1,NLCF)=0.0D0 NGCG=MELEFL.NUM(1,NLCF) NGCD=MELEFL.NUM(3,NLCF) NLCG=MLECEN.LECT(NGCG) NLCD=MLECEN.LECT(NGCD) NGS1=MELEFP.NUM(1,NLCF) NGS2=MELEFP.NUM(2,NLCF) NLS1=MLESOM.LECT(NGS1) NLS2=MLESOM.LECT(NGS2) SCNX=MPONOR.VPOCHA(NLCF,1) SCNY=MPONOR.VPOCHA(NLCF,2) SCN1X = SCNX SCN1Y = SCNY SURF=0.5D0*MPOSUR.VPOCHA(NLCF,1) SCNX=SCNX*SURF SCNY=SCNY*SURF C 3=IDIM+1 ICELL=(3*(NGCG -1))+1 XG=MCOORD.XCOOR(ICELL) YG=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGCD -1))+1 XD=MCOORD.XCOOR(ICELL) YD=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGCF -1))+1 XF=MCOORD.XCOOR(ICELL) YF=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS1 -1))+1 XS1=MCOORD.XCOOR(ICELL) YS1=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS2 -1))+1 XS2=MCOORD.XCOOR(ICELL) YS2=MCOORD.XCOOR(ICELL+1) XLONG = (((XS1-XS2)**2) + ((YS1-YS2)**2)) XLONG = SQRT(XLONG) IG1 = 1 ID1 = 1 IG2 = 1 ID2 = 1 MELTFA = IMELEM(NLCG) NBF = NBFACEL(NLCG) IF (NLCG.GT.NMAI1) THEN NGAUX = NLCG - NMAI1 ELSE NGAUX = NLCG ENDIF DO J = 1,NBF N1 = MELTFA.NUM(J,NGAUX) NL1 = MLEFA.LECT(N1) NSOM1 = MELEFP.NUM(1,NL1) NSOM2 = MELEFP.NUM(2,NL1) IF ((NSOM1.EQ.NGS1).OR.(NSOM2.EQ.NGS1)) THEN ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS1 -1))+1 XS1=MCOORD.XCOOR(ICELL) YS1=MCOORD.XCOOR(ICELL+1) c on corrige pour VFSYM IF (NBF.EQ.3) THEN XF = ((2.D0*XF/3.D0) + (XS1/3.D0)) YF = ((2.D0*YF/3.D0) + (YS1/3.D0)) ENDIF VECXG1(IG1) = -(YF - YG) VECYG1(IG1) = (XF - XG) VX = (XG - XS1) VY = (YG - YS1) PSCA = (VX*VECXG1(IG1)) + (VY*VECYG1(IG1)) IF (PSCA.LT.0.0D0) THEN VECXG1(IG1) = +(YF - YG) VECYG1(IG1) = -(XF - XG) ENDIF c ON REPERE l'INDICE IF ((NSOM2.NE.NGS2).AND.(NSOM1.NE.NGS2)) THEN INDG1 = IG1 NG1 = N1 ENDIF IG1 = IG1 + 1 ENDIF IF ((NSOM1.EQ.NGS2).OR.(NSOM2.EQ.NGS2)) THEN ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS2 -1))+1 XS2=MCOORD.XCOOR(ICELL) YS2=MCOORD.XCOOR(ICELL+1) c on corrige pour VFSYM IF (NBF.EQ.3) THEN XF = ((2.D0*XF/3.D0) + (XS2/3.D0)) YF = ((2.D0*YF/3.D0) + (YS2/3.D0)) ENDIF VECXG2(IG2) = -(YF - YG) VECYG2(IG2) = (XF - XG) VX = (XG - XS2) VY = (YG - YS2) PSCA = (VX*VECXG2(IG2)) + (VY*VECYG2(IG2)) IF (PSCA.LT.0.0D0) THEN VECXG2(IG2) = +(YF - YG) VECYG2(IG2) = -(XF - XG) ENDIF IF ((NSOM2.NE.NGS1).AND.(NSOM1.NE.NGS1)) THEN INDG2 = IG2 NG2 = N1 ENDIF IG2 = IG2 + 1 ENDIF ENDDO MELTFA = IMELEM(NLCD) NBF = NBFACEL(NLCD) IF (NLCD.GT.NMAI1) THEN NDAUX = NLCD -NMAI1 ELSE NDAUX = NLCD ENDIF DO J = 1,NBF N1 = MELTFA.NUM(J,NDAUX) NL1 = MLEFA.LECT(N1) NSOM1 = MELEFP.NUM(1,NL1) NSOM2 = MELEFP.NUM(2,NL1) IF ((NSOM1.EQ.NGS1).OR.(NSOM2.EQ.NGS1)) THEN ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS1 -1))+1 XS1=MCOORD.XCOOR(ICELL) YS1=MCOORD.XCOOR(ICELL+1) c on corrige pour VFSYM IF (NBF.EQ.3) THEN XF = ((2.D0*XF/3.D0) + (XS1/3.D0)) YF = ((2.D0*YF/3.D0) + (YS1/3.D0)) ENDIF VECXD1(ID1) = - (YF - YD) VECYD1(ID1) = (XF - XD) VX = (XD - XS1) VY = (YD - YS1) PSCA = (VX*VECXD1(ID1)) + (VY*VECYD1(ID1)) IF (PSCA.LT.0.0D0) THEN VECXD1(ID1) = +(YF - YD) VECYD1(ID1) = -(XF - XD) ENDIF IF ((NSOM2.NE.NGS2).AND.(NSOM1.NE.NGS2)) THEN INDD1 = ID1 ND1 = N1 ENDIF ID1 = ID1 + 1 ENDIF IF ((NSOM1.EQ.NGS2).OR.(NSOM2.EQ.NGS2)) THEN ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=MCOORD.XCOOR(ICELL+1) ICELL=(3*(NGS2 -1))+1 XS2=MCOORD.XCOOR(ICELL) YS2=MCOORD.XCOOR(ICELL+1) c on corrige pour VFSYM IF (NBF.EQ.3) THEN XF = ((2.D0*XF/3.D0) + (XS2/3.D0)) YF = ((2.D0*YF/3.D0) + (YS2/3.D0)) ENDIF VECXD2(ID2) = - (YF - YD) VECYD2(ID2) = (XF - XD) VX = (XD - XS2) VY = (YD - YS2) PSCA = (VX*VECXD2(ID2)) + (VY*VECYD2(ID2)) IF (PSCA.LT.0.0D0) THEN VECXD2(ID2) = +(YF - YD) VECYD2(ID2) = -(XF - XD) ENDIF IF ((NSOM2.NE.NGS1).AND.(NSOM1.NE.NGS1)) THEN INDD2 = ID2 ND2 = N1 ENDIF ID2 = ID2 + 1 ENDIF ENDDO AG1=ABS( ( (VECXG1(1)*VECYG1(2)) - & (VECYG1(1))*VECXG1(2)) ) AG2=ABS( ( (VECXG2(1)*VECYG2(2)) - & (VECYG2(1))*VECXG2(2)) ) AD1=ABS( ( (VECXD1(1)*VECYD1(2)) - & (VECYD1(1))*VECXD1(2)) ) AD2=ABS( ( (VECXD2(1)*VECYD2(2)) - & (VECYD2(1))*VECXD2(2)) ) c CALCUL DE MATRICE POUR LE NOEUD D INDICE NS1 COEF1 = ( (VECXG1(INDG1)*SCNX) + (VECYG1(INDG1)*SCNY) ) & / AG1 IAUX = 3 - INDG1 COEF2 = ( (VECXG1(IAUX)*SCNX) + (VECYG1(IAUX)*SCNY) ) & / AG1 COEF3 = ( (VECXD1(INDD1)*SCNX) + (VECYD1(INDD1)*SCNY) ) & / AD1 IAUX = 3 - INDD1 COEF4 = ( (VECXD1(IAUX)*SCNX) + (VECYD1(IAUX)*SCNY) ) & / AD1 MARQ = 0 DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.NG1) THEN IAFF = I5 IG1 = IAFF GOTO 62 ENDIF ENDDO 62 CONTINUE TRG1 = SCMB.MAT(IAFF,NLS1) MARQ = 0 DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.NGCF) THEN IAFF = I5 IG = IAFF IGNS1 = IAFF GOTO 63 ENDIF ENDDO 63 CONTINUE TRG = SCMB.MAT(IAFF,NLS1) TRGAUX = TRG TRGNS1 = TRG IAUX = 3 - INDG1 VXG1 = VECXG1(INDG1)/AG1 VXAU = VECXG1(IAUX)/AG1 VYG1 = VECYG1(INDG1)/AG1 VYAU = VECYG1(IAUX)/AG1 IF (ICHTE.GT.0) THEN IF (MPOTEN.VPOCHA(/2) .EQ.3) THEN c LE TENSEUR EST ANISOTROPE K11G = MPOTEN.VPOCHA(NLCG,1) K22G = MPOTEN.VPOCHA(NLCG,2) K21G = MPOTEN.VPOCHA(NLCG,3) K11D = MPOTEN.VPOCHA(NLCD,1) K22D = MPOTEN.VPOCHA(NLCD,2) K21D = MPOTEN.VPOCHA(NLCD,3) ELSEIF (MPOTEN.VPOCHA(/2) .EQ.1) THEN K11G = MPOTEN.VPOCHA(NLCG,1) K22G = K11G K21G = 0.0D0 K11D = MPOTEN.VPOCHA(NLCD,1) K22D = K11D K21D = 0.0D0 ELSE WRITE(6,*) 'TENSEUR NON PREVU' ENDIF c NLCG2 = MLECEN2.LECT(NGCG) c NLCD2 = MLECEN2.LECT(NGCD) c K11G = MPOTEN.VPOCHA(NLCG2,1) c K22G = MPOTEN.VPOCHA(NLCG2,2) c K21G = MPOTEN.VPOCHA(NLCG2,3) c K11D = MPOTEN.VPOCHA(NLCD2,1) c K22D = MPOTEN.VPOCHA(NLCD2,2) c K21D = MPOTEN.VPOCHA(NLCD2,3) c PRODUIT TENSEUR VECTEUR VXG1 = ((VECXG1(INDG1)/AG1)*K11G) + ((VECYG1(INDG1)/AG1)*K21G) VYG1 = ((VECXG1(INDG1)/AG1)*K21G) + ((VECYG1(INDG1)/AG1)*K22G) IAUX = 3 - INDG1 VXAU = ((VECXG1(IAUX)/AG1)*K11G) + ((VECYG1(IAUX)/AG1)*K21G) VYAU = ((VECXG1(IAUX)/AG1)*K21G) + ((VECYG1(IAUX)/AG1)*K22G) c WRite(6,*) 'K11G=',K11G,'K22G= ',K22G,'K21G=',K21G c WRite(6,*) 'K11D=',K11D,'K22D= ',K22D,'K21D=',K21D c WRITE(6,*) 'NLCF= ',NLCF,'VECXG1= ',VXG1 c IAUX = 3 - INDG1 c WRITE(6,*) 'NLCF= ',NLCF,'VXAU= ',VXAU c WRITE(6,*) 'NLCF= ',NLCF,'VECYG1= ',VYG1 c IAUX = 3 - INDG1 c WRITE(6,*) 'NLCF= ',NLCF,'VYAU= ',VYAU c WRITE(6,*) 'NLCF= ',NLCF,'VALD= ',VALD c WRITE(6,*) 'NLCF= ',NLCF,'VECXG1= ',VECXG1(INDG1) c IAUX = 3 - INDG1 c WRITE(6,*) 'NLCF= ',NLCF,'VXAU= ',VECXG1(IAUX) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG1= ',VECYG1(INDG1) c IAUX = 3 - INDG1 c WRITE(6,*) 'NLCF= ',NLCF,'VYAU= ',VECYG1(IAUX) c WRITE(6,*) 'NLCF= ',NLCF,'VALD= ',VALD ENDIF VAL = MPOCHP.VPOCHA(NLCG,1) VALG = VAL COEF1X = VXG1 COEF2X = VXAU COEF1Y = VYG1 COEF2Y = VYAU AUX = SCN1X* ( (COEF1X * (TRG - VAL)) + & (COEF2X * (TRG1 - VAL)) ) AUY = SCN1Y* ( (COEF1Y * (TRG - VAL)) + & (COEF2Y * (TRG1 - VAL)) ) MPOGRA.VPOCHA(NLCF,1) = AUX + AUY c WRITE(6,*) 'NLCF= ',NLCF c WRITE(6,*) 'COEF1X= ',COEF1X c WRITE(6,*) 'COEF2X= ',COEF2X c WRITE(6,*) 'COEF1Y= ',COEF1Y c WRITE(6,*) 'COEF2Y= ',COEF2Y INDIC = IPO3.POINT33(NLS1) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS1) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS1) CX = MATR1.MAT2(IG,JAUX) CY = MATR1.MAT2(IG,JAUX) INDAUX = INDIC.NU(IG,JAUX) MELEME.NUM(JAUX+1,NLCF) = INDAUX IF (INDAUX.EQ.NGCG) THEN CX = CX - 1.D0 CY = CY - 1.D0 INDGA = JAUX +1 ENDIF IF (INDAUX.EQ.NGCD) THEN INDDR = JAUX + 1 ENDIF MELVA1.VELCHE(JAUX+1,NLCF) = & (COEF1X*CX*SCN1X) + (COEF1Y*CY*SCN1Y) ENDDO ITROUVE = 0 DO ICOUR = 1,TAB.ID(NLS1) IA = ICOUR J1 = INDIC.NU(IG1,IA) DO IAUX2 = 2,TAB.ID(NLS1)+1 J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 ITROUVE = 1 GOTO 511 ENDIF ENDDO 511 CONTINUE IF (ITROUVE.EQ.0) THEN WRITE(6,*) 'PROBLEME ALGO1' WRITE(6,*) 'INDIC= ',INDIC DO IAUX2 = 1,TAB.ID(NLS1)+1 J2 = MELEME.NUM(IAUX2,NLCF) WRITE(6,*) 'NLCF= ',NLCF,'IAUX2=',IAUX2,'MELEME= ',J2 WRITE(6,*) 'IG= ',IG,'IAUX2=',IAUX2, & 'NLS1= ',NLS1,'INDIC=',INDIC.NU(IG,IAUX2) WRITE(6,*) 'IG1= ',IG1,'IAUX2=',IAUX2, & 'NLS1= ',NLS1,'INDIC= ',INDIC.NU(IG1,IAUX2) ENDDO CALL ERREUR(5) ENDIF CX = MATR1.MAT2(IG1,IA) CY = MATR1.MAT2(IG1,IA) IF (J1.EQ.NGCG) THEN CX = CX - 1.D0 CY = CY - 1.D0 ENDIF MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) + & (COEF2X*CX*SCN1X) + (COEF2Y*CY*SCN1Y) ENDDO c DO IAUX2 = 2,TAB.ID(NLS1)+1 c J1 = MELEME.NUM(IAUX2,NLCF) c WRITE(6,*) 'NLCF= ',NLCF,'IAUX2= ',IAUX2 c WRITE(6,*) 'MELEME2= ',J1 c COEF1 = MELVA1.VELCHE(IAUX2,NLCF) c COEF2 = MELVA2.VELCHE(IAUX2,NLCF) c WRITE(6,*) 'COEF1= ',COEF1,'COEF2= ',COEF2 c ENDDO * POUR NLS1 c DO IAUX2 = 1,INDLI.ID(NLS1) c WRITE(6,*) 'NLS1= ',NLS1,'INDIC= ',INDIC(ID,IAUX2) c WRITE(6,*) 'NLS1= ',NLS1,'INDIC= ',INDIC(ID1,IAUX2) c ENDDO c DO IAUX2 = 2,TAB.ID(NLS1)+1 c J1 = MELEME.NUM(IAUX2,NLCF) c WRITE(6,*) 'NLCF= ',NLCF,'IAUX2= ',IAUX2 c WRITE(6,*) 'MELEME3= ',J1 c COEF1 = MELVA1.VELCHE(IAUX2,NLCF) c COEF2 = MELVA2.VELCHE(IAUX2,NLCF) c WRITE(6,*) 'COEF1= ',COEF1,'COEF2= ',COEF2 c ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD INDIC = IPO3.POINT33(NLS2) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS2) SEGACT MATR1 *MOD c CALCUL DE MATRICE POUR LE NOEUD D INDICE NS2 COEF1 = ( (VECXG2(INDG2)*SCNX) + (VECYG2(INDG2)*SCNY) ) & / AG2 IAUX = 3 - INDG2 COEF2 = ( (VECXG2(IAUX)*SCNX) + (VECYG2(IAUX)*SCNY) ) & / AG2 COEF3 = ( (VECXD2(INDD2)*SCNX) + (VECYD2(INDD2)*SCNY) ) & / AD2 IAUX = 3 - INDD2 COEF4 = ( (VECXD2(IAUX)*SCNX) + (VECYD2(IAUX)*SCNY) ) & / AD2 DO I5 = 1,INDLI.ID(NLS2) INDAUX = IND2.NUME(I5,NLS2) IF (INDAUX.EQ.NG2) THEN IAFF = I5 IG2 = IAFF GOTO 411 ENDIF ENDDO 411 CONTINUE TRG2 = SCMB.MAT(IAFF,NLS2) MARQ = 0 DO I5 = 1,INDLI.ID(NLS2) INDAUX = IND2.NUME(I5,NLS2) IF (INDAUX.EQ.NGCF) THEN IAFF = I5 IG = IAFF IGNS2 = IAFF GOTO 635 ENDIF ENDDO 635 CONTINUE TRG = SCMB.MAT(IAFF,NLS2) TRGNS2 = TRG IAUX = 3 - INDG2 VXG2 = VECXG2(INDG2)/AG2 VXAU = VECXG2(IAUX)/AG2 VYG2 = VECYG2(INDG2)/AG2 VYAU = VECYG2(IAUX)/AG2 c WRITE(6,*) 'NLCF= ',NLCF,'VECXG2= ',VECXG2(INDG2) c IAUX = 3 - INDG2 c WRITE(6,*) 'NLCF= ',NLCF,'VXAU= ',VECXG2(IAUX) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG2= ',VECYG2(INDG2) c IAUX = 3 - INDG2 c WRITE(6,*) 'NLCF= ',NLCF,'VYAU= ',VECYG2(IAUX) IF (ICHTE.GT.0) THEN IF (MPOTEN.VPOCHA(/2) .EQ.3) THEN c LE TENSEUR EST ANISOTROPE K11G = MPOTEN.VPOCHA(NLCG,1) K22G = MPOTEN.VPOCHA(NLCG,2) K21G = MPOTEN.VPOCHA(NLCG,3) K11D = MPOTEN.VPOCHA(NLCD,1) K22D = MPOTEN.VPOCHA(NLCD,2) K21D = MPOTEN.VPOCHA(NLCD,3) ELSEIF (MPOTEN.VPOCHA(/2) .EQ.1) THEN K11G = MPOTEN.VPOCHA(NLCG,1) K22G = K11G K21G = 0.0D0 K11D = MPOTEN.VPOCHA(NLCD,1) K22D = K11D K21D = 0.0D0 ELSE WRITE(6,*) 'TENSEUR NON PREVU' ENDIF c NLCG2 = MLECEN2.LECT(NGCG) c NLCD2 = MLECEN2.LECT(NGCD) c K11G = MPOTEN.VPOCHA(NLCG2,1) c K22G = MPOTEN.VPOCHA(NLCG2,2) c K21G = MPOTEN.VPOCHA(NLCG2,3) c c K11D = MPOTEN.VPOCHA(NLCD2,1) c K22D = MPOTEN.VPOCHA(NLCD2,2) c K21D = MPOTEN.VPOCHA(NLCD2,3) c PRODUIT TENSEUR VECTEUR VXG2 = ((VECXG2(INDG2)/AG2)*K11G) + ((VECYG2(INDG2)/AG2)*K21G) VYG2 = ((VECXG2(INDG2)/AG2)*K21G) + ((VECYG2(INDG2)/AG2)*K22G) IAUX = 3 - INDG2 VXAU = ((VECXG2(IAUX)/AG2)*K11G) + ((VECYG2(IAUX)/AG2)*K21G) VYAU = ((VECXG2(IAUX)/AG2)*K21G) + ((VECYG2(IAUX)/AG2)*K22G) c WRite(6,*) 'K11G=',K11G,'K22G= ',K22G,'K21G=',K21G c WRite(6,*) 'K11D=',K11D,'K22D= ',K22D,'K21D=',K21D c WRITE(6,*) 'NLCF= ',NLCF,'VECXG2= ',VXG2 c IAUX = 3 - INDG2 c WRITE(6,*) 'NLCF= ',NLCF,'VXAU= ',VXAU c WRITE(6,*) 'NLCF= ',NLCF,'VECYG2= ',VYG2 c IAUX = 3 - INDG2 c WRITE(6,*) 'NLCF= ',NLCF,'VYAU= ',VYAU ENDIF COEF1X = VXG2 COEF2X = VXAU COEF1Y = VYG2 COEF2Y = VYAU AUX = SCN1X* ( (COEF1X * (TRG - VAL)) + & (COEF2X * (TRG2 - VAL)) ) AUY = SCN1Y* ( (COEF1Y * (TRG - VAL)) + & (COEF2Y * (TRG2 - VAL)) ) MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) + & AUX + AUY ANCX = MPOGRA.VPOCHA(NLCF,1) * POUR NLS2 c WRITE(6,*) 'COEF1X= ',COEF1X c WRITE(6,*) 'COEF2X= ',COEF2X c WRITE(6,*) 'COEF1Y= ',COEF1Y c WRITE(6,*) 'COEF2Y= ',COEF2Y ICON = TAB.ID(NLS1)+1 DO JAUX = 1,TAB.ID(NLS2) J1 = INDIC.NU(IG,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 59 ENDIF ENDDO ICON = ICON +1 IAUX = ICON IF (IAUX.GT.NBNN) THEN WRITE(6,*) 'NBNN TROP PETIT' CALL ERREUR(5) ENDIF 59 CONTINUE CX = MATR1.MAT2(IG,JAUX) CY = MATR1.MAT2(IG,JAUX) INDAUX = INDIC.NU(IG,JAUX) MELEME.NUM(IAUX,NLCF) = INDAUX IF (INDAUX.EQ.NGCG) THEN CX = CX - 1.0D0 CY = CY - 1.0D0 ENDIF MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) + & (COEF1X*CX*SCN1X) + (COEF1Y*CY*SCN1Y) ENDDO c DO IAUX2 = 2,ICON c J1 = MELEME.NUM(IAUX2,NLCF) c WRITE(6,*) 'NLCF= ',NLCF,'IAUX2= ',IAUX2 c WRITE(6,*) 'MELEME4= ',J1 c COEF1 = MELVA1.VELCHE(IAUX2,NLCF) c COEF2 = MELVA2.VELCHE(IAUX2,NLCF) c WRITE(6,*) 'COEF1= ',COEF1,'COEF2= ',COEF2 c ENDDO ITROUVE = 0 DO ICOUR = 1,TAB.ID(NLS2) IA = ICOUR J1 = INDIC.NU(IG2,IA) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 ITROUVE = 1 GOTO 577 ENDIF ENDDO 577 CONTINUE IF (ITROUVE.EQ.0) THEN WRITE(6,*) 'PROBLEME ALGO3' WRITE(6,*) 'INDIC= ',INDIC DO IAUX2 = 2,ICON J1 = MELEME.NUM(IAUX2,NLCF) WRITE(6,*) 'NLCF= ',NLCF,'IAUX2=',IAUX2,'MELEME= ',J1 WRITE(6,*) 'IG= ',IG,'IAUX2=',IAUX2,'INDIC=', & 'NLS2= ',NLS2,INDIC.NU(IG,IAUX2) WRITE(6,*) 'IG2= ',IG2,'IAUX2=',IAUX2,'INDIC= ', & 'NLS2= ',NLS2,INDIC.NU(IG2,IAUX2) ENDDO CALL ERREUR(5) ENDIF CX = MATR1.MAT2(IG2,IA) CY = MATR1.MAT2(IG2,IA) IF (J1.EQ.NGCG) THEN CX = CX - 1.0D0 CY = CY - 1.D0 ENDIF MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) + & (COEF2X*CX*SCN1X) + (COEF2Y*CY*SCN1Y) ENDDO MPOGRA.VPOCHA(NLCF,1) = (-0.5D0)* MPOGRA.VPOCHA(NLCF,1) ISUP = ICON DO J= ISUP+1,NBNN MELVA1.VELCHE(J,NLCF) = 0.0D0 MELEME.NUM(J,NLCF) = MELEME.NUM(ISUP,NLCF) ENDDO DO J= 1,NBNN MELVA1.VELCHE(J,NLCF) = (-0.5D0*MELVA1.VELCHE(J,NLCF)) ENDDO C ON RAJOUTE LE CONVECTIF IF (ICHCO.GT.0) THEN C BOUCLE POUR CALCUER INDGA,INDDR INDFR = 0 DO J= 1,ISUP IF (MELEME.NUM(J,NLCF).EQ.NGCG) INDGA = J IF (MELEME.NUM(J,NLCF).EQ.NGCD) INDDR = J IF (MELEME.NUM(J,NLCF).EQ.NGCF) INDFR = J ENDDO UN = MPOVCO.VPOCHA(NLCF,1) c WRITE(6,*) 'UN= ',UN C OPTION CENTRE IF (IOP.EQ.2) THEN IF (NLCD.NE.NLCG) THEN VAL = 0.5D0*(MPOCHP.VPOCHA(NLCG,1) + & MPOCHP.VPOCHA(NLCD,1))*UN MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ELSE C CONDITIONS AUX LIMITES : TRACE CALCULEE PAR LA DIFFUSION TRACE = 0.5*(TRGNS1+TRGNS2) c WRITE(6,*) 'NLCF= ',NLCF,'TRACE= ',TRACE VAL = TRACE*UN MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ENDIF C ON COMPLETE MELVA1 POUR LE CONVECTIF IF (NLCD.NE.NLCG) THEN MELVA1.VELCHE(INDGA,NLCF) = MELVA1.VELCHE(INDGA,NLCF) - & (0.5D0*UN) MELVA1.VELCHE(INDDR,NLCF) = MELVA1.VELCHE(INDDR,NLCF) - & (0.5D0*UN) C CONDITION AUX LIMITES : ON RAJOUTE LES DEPENDANCES DES TRACES c POUR LES CONDITIONS MIXTES OU DE NEUMAN ELSE NLFCL = MLENCL.LECT(NGCF) C ON RAJOUTE CECI POUR L OPTION GRADGEO IF (NLFCL.NE.0) THEN MELVA1.VELCHE(ICON+1,NLCF) = - UN MELEME.NUM(ICON+1,NLCF) = NGCF ENDIF IF (NLFCL.EQ.0) THEN INDIC = IPO3.POINT33(NLS1) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS1) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS1) J1 = INDIC.NU(IGNS1,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5459 ENDIF ENDDO 5459 CONTINUE CX = MATR1.MAT2(IGNS1,JAUX) MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - (UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD INDIC = IPO3.POINT33(NLS2) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS2) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS2) J1 = INDIC.NU(IGNS2,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5988 ENDIF ENDDO 5988 CONTINUE CX = MATR1.MAT2(IGNS2,JAUX) c WRITE(6,*) 'NLCF= ',NLCF,'CX= ',CX MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - (UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD ENDIF ENDIF C OPTION UPWIND ELSEIF (IOP.EQ.1) THEN IF (NLCD.NE.NLCG) THEN IF (UN.GT.0.0D0) THEN VAL = (MPOCHP.VPOCHA(NLCG,1)*UN) ELSE VAL = (MPOCHP.VPOCHA(NLCD,1)*UN) ENDIF MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ELSE C CONDITIONS AUX LIMITES : TRACE CALCULEE PAR LA DIFFUSION C CONDITIONS AUX LIMITES : TRACE CALCULEE PAR LA DIFFUSION IF (UN.GT.0.0D0) THEN VAL = (MPOCHP.VPOCHA(NLCG,1)*UN) ELSE TRACE = 0.5*(TRGNS1+TRGNS2) VAL = TRACE*UN ENDIF MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ENDIF C ON COMPLETE MELVA1 POUR LE CONVECTIF IF (NLCD.NE.NLCG) THEN IF (UN.GT.0.0D0) THEN MELVA1.VELCHE(INDGA,NLCF) = MELVA1.VELCHE(INDGA,NLCF) - & (UN) ELSE MELVA1.VELCHE(INDDR,NLCF) = MELVA1.VELCHE(INDDR,NLCF) - & (UN) ENDIF C CONDITION AUX LIMITES : ON RAJOUTE LES DEPENDANCES DES TRACES c POUR LES CONDITIONS MIXTES OU DE NEUMAN ELSE IF (UN.GT.0.0D0) THEN MELVA1.VELCHE(INDGA,NLCF) = MELVA1.VELCHE(INDGA,NLCF) - & (UN) ELSE NLFCL = MLENCL.LECT(NGCF) C ON RAJOUTE CECI POUR L OPTION GRADGEO IF (NLFCL.NE.0) THEN MELVA1.VELCHE(ICON+1,NLCF) = - UN MELEME.NUM(ICON+1,NLCF) = NGCF ENDIF c IF (1.EQ.0) THEN IF (NLFCL.EQ.0) THEN INDIC = IPO3.POINT33(NLS1) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS1) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS1) J1 = INDIC.NU(IGNS1,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5449 ENDIF ENDDO 5449 CONTINUE CX = MATR1.MAT2(IGNS1,JAUX) MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - (UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD INDIC = IPO3.POINT33(NLS2) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS2) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS2) J1 = INDIC.NU(IGNS2,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5989 ENDIF ENDDO 5989 CONTINUE CX = MATR1.MAT2(IGNS2,JAUX) MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - (UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD c ENDIF ENDIF ENDIF ENDIF C OPTION UPWICENT ELSEIF (IOP.EQ.3) THEN IF (NLCD.NE.NLCG) THEN c CALCUL DE THETA c DANS MELVA1 : INDDR EST LA NUM de NLCD C COEFD = MELVA1(INDDR,NLCF) COEFD = MELVA1.VELCHE(INDDR,NLCF) THETA = 0.5D0 EXPR1 = (THETA*UN) + COEFD EXPR2 = ((1.D0 - THETA)*UN) - COEFD IF (EXPR1.LT.0.0D0) THEN IF (ABS(UN) .GT. 1e-20) THEN THETA = - (0.5*COEFD / UN) ENDIF ENDIF IF (EXPR2.GT.0.0D0) THEN IF (ABS(UN) .GT. 1e-20) THEN THETA = 1.D0 - (0.5*COEFD / UN) ENDIF ENDIF THETA = MIN(1.D0,THETA) THETA = MAX(0.D0,THETA) c WRITE(6,*) 'THETA= ',THETA c WRITE(6,*) 'COEFD= ',COEFD c WRITE(6,*) 'UN= ',UN c WRITE(6,*) 'EXPR1= ',EXPR1 c WRITE(6,*) 'EXPR2= ',EXPR2 VAL = (MPOCHP.VPOCHA(NLCG,1)*UN*THETA) + & (MPOCHP.VPOCHA(NLCD,1)*(1.D0-THETA)*UN) MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ELSE C CONDITIONS AUX LIMITES : TRACE CALCULEE PAR LA DIFFUSION IF (INDFR.NE.0) THEN COEFD = MELVA1.VELCHE(INDFR,NLCF) ELSE COEFD = 0.0D0 ENDIF C C'EST LE THETA LE MIEUX ADPATE POUR LES CONDITIONS AUX LIMITES THETA = 0.0D0 EXPR1 = (THETA*UN) + COEFD EXPR2 = ((1.D0 - THETA)*UN) - COEFD IF (EXPR1.LT.0.0D0) THEN IF (ABS(UN) .GT. 1e-20) THEN THETA = - (0.5*COEFD / UN) ENDIF ENDIF IF (EXPR2.GT.0.0D0) THEN IF (ABS(UN) .GT. 1e-20) THEN THETA = 1.D0 - (0.5*COEFD / UN) ENDIF ENDIF THETA = MIN(1.D0,THETA) THETA = MAX(0.D0,THETA) TRACE = 0.5*(TRGNS1+TRGNS2) VAL = (MPOCHP.VPOCHA(NLCG,1)*UN*THETA) + & (TRACE*(1.D0-THETA)*UN) MPOGRA.VPOCHA(NLCF,1) = MPOGRA.VPOCHA(NLCF,1) - VAL ENDIF C ON COMPLETE MELVA1 POUR LE CONVECTIF IF (NLCD.NE.NLCG) THEN MELVA1.VELCHE(INDGA,NLCF) = MELVA1.VELCHE(INDGA,NLCF) - & (THETA*UN) MELVA1.VELCHE(INDDR,NLCF) = MELVA1.VELCHE(INDDR,NLCF) - & ((1.D0-THETA)*UN) C CONDITION AUX LIMITES : ON RAJOUTE LES DEPENDANCES DES TRACES c POUR LES CONDITIONS MIXTES OU DE NEUMAN ELSE MELVA1.VELCHE(INDGA,NLCF) = MELVA1.VELCHE(INDGA,NLCF) - & (THETA*UN) NLFCL = MLENCL.LECT(NGCF) C ON RAJOUTE CECI POUR L OPTION GRADGEO c IF (NLFCL.NE.0) THEN c MELVA1.VELCHE(ICON+1,NLCF) = - UN c MELEME.NUM(ICON+1,NLCF) = NGCF c ENDIF c IF (1.EQ.0) THEN c IF (NLFCL.EQ.0) THEN INDIC = IPO3.POINT33(NLS1) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS1) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS1) J1 = INDIC.NU(IGNS1,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5549 ENDIF ENDDO 5549 CONTINUE CX = MATR1.MAT2(IGNS1,JAUX) MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - ((1.D0-THETA)*UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD INDIC = IPO3.POINT33(NLS2) SEGACT INDIC *MOD MATR1 = IPO2.POINT(NLS2) SEGACT MATR1 *MOD DO JAUX = 1,TAB.ID(NLS2) J1 = INDIC.NU(IGNS2,JAUX) DO IAUX2 = 2,ICON J2 = MELEME.NUM(IAUX2,NLCF) IF (J1.EQ.J2) THEN IAUX = IAUX2 GOTO 5289 ENDIF ENDDO 5289 CONTINUE CX = MATR1.MAT2(IGNS2,JAUX) MELVA1.VELCHE(IAUX,NLCF) = MELVA1.VELCHE(IAUX,NLCF) & - ((1.D0-THETA)*UN*CX*0.5D0) ENDDO SEGDES INDIC *MOD SEGDES MATR1 *MOD c ENDIF c ENDIF ENDIF ENDIF ENDIF c DO IAUX2 = 2,ICON c J1 = MELEME.NUM(IAUX2,NLCF) c WRITE(6,*) 'NLCF= ',NLCF,'IAUX2= ',IAUX2 c WRITE(6,*) 'MELEME6= ',J1 c COEF1 = MELVA1.VELCHE(IAUX2,NLCF) c WRITE(6,*) 'COEF1= ',COEF1 c ENDDO c WRITE(6,*) 'MPOGRA= ',MPOGRA.VPOCHA(NLCF,1) c IF (ICHTE.GT.0) THEN c IAUX = 3 - INDD1 c XLONGD = (VECXD1(IAUX)*VECXD1(IAUX)) + c & (VECYD1(IAUX)*VECYD1(IAUX)) c XLONGD = XLONGD**0.5 c IAUX = 3 - INDG1 c XLONGG = (VECXG1(IAUX)*VECXG1(IAUX)) + c & (VECYG1(IAUX)*VECYG1(IAUX)) c XLONGG = XLONGG**0.5 c ANCX = ((K11G/XLONGG)*VALG) + ((K11D/XLONGD)*VALD) c ANCX = ANCX/((K11G/XLONGG) + (K11D/XLONGD)) c ANCX = (VALD - VALG) / c & ( (XLONGG/K11G) + (XLONGD/K11D)) c ANCX = (VALD - VALG) / c & ( (XLONGG) + (XLONGD)) c MPOGRA.VPOCHA(NLCF,1) = SCN1X*ANCX c MPOGRA.VPOCHA(NLCF,2) = SCN1Y*ANCX c ENDIF c IF (NGCG.NE.NGCD) THEN c AA = MPOGRA.VPOCHA(NLCF,1) + AA c BB = MPOGRA.VPOCHA(NLCF,2) + BB c ENDIF IF (1.EQ.0) THEN WRITE(6,*) 'NLCF= ',NLCF,'GR1= ',MPOGRA.VPOCHA(NLCF,1) WRITE(6,*) 'NLCF= ',NLCF,'GR2= ',MPOGRA.VPOCHA(NLCF,2) WRITE(6,*) 'NLCF= ',NLCF,'TRG= ',TRGAUX WRITE(6,*) 'NLCF= ',NLCF,'TRG1= ',TRG1 WRITE(6,*) 'NLCF= ',NLCF,'TRG= ',TRG WRITE(6,*) 'NLCF= ',NLCF,'TRG2= ',TRG2 WRITE(6,*) 'NLCF= ',NLCF,'VAL= ',VAL VALD = MPOCHP.VPOCHA(NLCD,1) WRITE(6,*) 'NLCF= ',NLCF,'VALD= ',VALD WRITE(6,*) 'NLCF= ',NLCF,'VECXG1= ',VECXG1(INDG1) IAUX = 3 - INDG1 WRITE(6,*) 'NLCF= ',NLCF,'VECXG1= ',VECXG1(IAUX) WRITE(6,*) 'NLCF= ',NLCF,'VECYG1= ',VECYG1(INDG1) IAUX = 3 - INDG1 WRITE(6,*) 'NLCF= ',NLCF,'VECYG1= ',VECYG1(IAUX) WRITE(6,*) 'NLCF= ',NLCF,'VECXG2= ',VECXG2(INDG2) IAUX = 3 - INDG2 WRITE(6,*) 'NLCF= ',NLCF,'VECXG2= ',VECXG2(IAUX) WRITE(6,*) 'NLCF= ',NLCF,'VECYG2= ',VECYG2(INDG2) IAUX = 3 - INDG2 WRITE(6,*) 'NLCF= ',NLCF,'VECYG2= ',VECYG2(IAUX) WRite(6,*) 'AG1=',AG1 WRite(6,*) 'AG2=',AG2 WRite(6,*) 'AD1=',AD1 WRite(6,*) 'AD2=',AD2 ENDIF c CALCUL DE MATRICE POUR LE NOEUD D INDICE NS2 NAUX2 = MAX(NAUX2,ICON) NMOY = NMOY + ICON SEGDES INDIC *MOD SEGDES MATR1 *MOD ENDDO NMOY = NMOY/(NFAC*1.D0) c WRITE(6,*) 'NAUX2= ',NAUX2 c WRITE(6,*) 'NMOY2= ',NMOY IF (NBSO.EQ.2) THEN SEGDES IPT1 SEGDES IPT2 ENDIF K7 = NFAC K8 = NAUX2 SEGINI ITAB c ON SUPPRIME LES ZEROS INTERIEURS IREP = 2 DO NLCF = 1,NFAC IREP = 2 NFIN = NAUX2 C ON CALCULE D4ABORD LE MAXIMUM DE LA LIGNE AUXMA = 0.0 DO J=2,NFIN AUXRE = ABS(MELVA1.VELCHE(J,NLCF)) AUXMA = MAX(AUXRE,AUXMA) ENDDO ITAB.XMAX(NLCF) = AUXMA DO J=IREP, NFIN IF (ABS(MELVA1.VELCHE(J,NLCF)).Le.(1e-14*AUXMA)) THEN DO K=1,NFIN-J AUX = MELVA1.VELCHE(J+K,NLCF) C ON DECALE DE K CRAN IF (ABS(AUX).gt.(1e-14*AUXMA)) THEN DO I=J,NFIN - K MELVA1.VELCHE(I,NLCF) = MELVA1.VELCHE(I+K,NLCF) MELEME.NUM(I,NLCF) = MELEME.NUM(I+K,NLCF) ENDDO C MISE A ZERO DES TERMES DU BOUT DO I=NFIN-K+1,NFIN MELVA1.VELCHE(I,NLCF) = 0.0 ENDDO GOTO 2000 ENDIF ENDDO 2000 CONTINUE ENDIF ENDDO ENDDO c TABLEAU CALCULANT LE NOMBRE DE VOISIN NON NUL POUR CHAQUE FACE NMOY = 0 INF = NAUX2 DO NLCF = 1,NFAC IREC = 3 DO J=NAUX2,1,-1 IF (ABS(MELVA1.VELCHE(J,NLCF)).gt. & (1e-14*ITAB.XMAX(NLCF))) THEN IREC = J GOTO 1111 ENDIF ENDDO 1111 CONTINUE ITAB.TABL(NLCF) = IREC NMOY = NMOY + IREC ENDDO NMOY = NMOY/(NFAC*1.D0) c WRITE(6,*) 'NEWMOY2= ',NMOY c WRITE(6,*) 'NFAC= ',NFAC C TAB1(U) TABLEAU QUI CONTIENT LE NOMBRE DE FACE AYANT U VOISIN NMAX = 0 DO ICOUR = 1,NAUX2 ITAB.TABL1(ICOUR) = 0 ENDDO DO NLCF = 1,NFAC ICOUR = ITAB.TABL(NLCF) ITAB.TABL1(ICOUR)=ITAB.TABL1(ICOUR) + 1 ENDDO C ON COMPTE LE NOMVRE DE SOUS DOMAINE NTSOUS = 0 DO ICOUR = 1,NAUX2 IF (ITAB.TABL1(ICOUR) .NE.0) NTSOUS = NTSOUS +1 ENDDO C IPOS INDICE DE LA PREMIERE FACE AYANT I VOISIN C ICOUR INDICE COURANT INITIALISE A IPOS ITAB.IPOS(1) = 1 DO I = 2,NAUX2 ITAB.IPOS(I) = ITAB.IPOS(I-1) + ITAB.TABL1(I-1) ITAB.ICOURANT(I) = ITAB.IPOS(I) ENDDO c TABL2 TABLEAU QUI RANGE DANS L4ODRES DES SOUS DOMAINES LES FACES NLCF DO I =1,NFAC IHELP = ITAB.TABL(I) IAUX = ITAB.ICOURANT(IHELP) ITAB.TABL2(IAUX) = I IAUX2 = ITAB.TABL(I) ITAB.ICOURANT(IAUX2) = ITAB.ICOURANT(IAUX2) + 1 ENDDO c WRITE(6,*) 'NTSOUS= ',NTSOUS C**** Initialisation du MCHELM C N1=NTSOUS N2=1 N3=6 L1=8 SEGINI MCHEL2 MCHEL2.TITCHE='Gradient' MCHEL2.IFOCHE=IFOUR C ISOUS=0 NBSOUS=0 NBREF=0 DO I1 = 1, NAUX2, 1 NBELEM=ITAB.TABL1(I1) IF(NBELEM .GT. 0)THEN ISOUS=ISOUS+1 NBNN=I1 SEGINI IPT8 C ITYPEL=32 -> 'POLY' ITYPEL=32 MCHEL2.IMACHE(ISOUS)=IPT8 SEGINI MCHAM2 MCHEL2.ICHAML(ISOUS)=MCHAM2 MCHAM2.NOMCHE(1)='SCAL' MCHAM2.TYPCHE(1)='REAL*8 ' N1PTEL=NBNN N1EL=NBELEM N2PTEL=0 N2EL=0 SEGINI MELVA2 MCHAM2.IELVAL(1)=MELVA2 C WRITE(6,*) 'I1= ',I1 C WRITE(6,*) 'NBELEM= ',NBELEM C WRITE(6,*) 'NBNN= ',NBNN DO I2=1,NBELEM,1 DO I3=1,NBNN,1 IAUX = ITAB.IPOS(I1) IP= ITAB.TABL2(I2 + IAUX - 1) c IAUX = MELEME.NUM(I3,1) c IAUX2 = MLEFA.LECT c WRITE(6,*) 'IP= ',IP,'I1= ',I1,'I2=',I2,'I3=',I3 c WRITE(6,*) 'ISOUS= ',ISOUS,'MELEME= ',MELEME.NUM(I3,IP) IPT8.NUM(I3,I2)= MELEME.NUM(I3,IP) MELVA2.VELCHE(I3,I2)=MELVA1.VELCHE(I3,IP) ENDDO ENDDO C SEGDES MCHAM2 SEGDES IPT8 SEGDES MELVA2 ENDIF ENDDO SEGDES MCHEL2 SEGDES ITAB SEGDES MCHAML SEGDES MELEME SEGDES MELVA1 SEGDES MCHELM SEGSUP IPO3 SEGSUP IPO2 SEGSUP INDLI SEGSUP TAB SEGSUP IND2 SEGSUP IND SEGSUP IND22 SEGSUP VAL1 SEGSUP VAL2 SEGSUP SCMB 9999 CONTINUE RETURN END