C VFSYM2 SOURCE OF166741 24/12/13 21:17:36 12097 SUBROUTINE VFSYM2(MELEFA,MELEFL,MLECEN,MELEFP,MLESOM,MPONOR, & MPOSUR,MELTFA,MLEFA,MPOTEN,MPOCHP,MLENCL, & MLENNE,MLENMI,MPOVCL, & MPOVNE,MPOVMI,ICHTE,ICHCL,ICHNE,IPO2,SCMB,INDLI, & TAB,VAL1,VAL2,IND22,IND2,IND,NBFAC,NSOMM,NBMAX) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : NORV2 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 PPARAM -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,MPOVNE.MPOVAL,MPOVMI.MPOVAL POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLESOM.MLENTI, & MLEFA.MLENTI,MLENNE.MLENTI,MLENMI.MLENTI -INC SMCHAML INTEGER NBNN,NBREF 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 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 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, & QIMPX,QIMPY,QIMPS,XLAMBDA1,XLAMBDA2 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 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 INDICE3 INTEGER IND3(K1,K2,K3) ENDSEGMENT SEGMENT REP INTEGER ID(K3) ENDSEGMENT POINTEUR TAB.REP,INDLI.REP INTEGER K5 SEGMENT NBFAC INTEGER NBFACEL(K5) INTEGER IMELEM(K5) ENDSEGMENT c CALCUL DES DIFFERENTS POINTEURS A ACTIVER DANS POUR PLUSIIEURS c SOUS DOMAINE MAUX = MELTFA NMAI1 = 0 NBSO = MAX(1,MELTFA.LISOUS(/1)) c WRITE(6,*) 'NBSO= ',NBSO 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 C SEGMENT SERVANT A UN PRECALCUL DE NBMAX K3 = NSOMM SEGINI INDLI SEGINI TAB DO I = 1,K3 INDLI.ID(I) = 0 TAB.ID(I) = 0 ENDDO NFAC=MELEFL.NUM(/2) NBMAX = 0 C PRECALCUL DE NBMAX DO NLCF= 1, NFAC, 1 c WRITE(6,*) 'NLCF= ',NLCF 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) INDLI.ID(NLS1) = INDLI.ID(NLS1) + 1 NBMAX = MAX(NBMAX,INDLI.ID(NLS1)) INDLI.ID(NLS2) = INDLI.ID(NLS2) + 1 NBMAX = MAX(NBMAX,INDLI.ID(NLS2)) ENDDO SEGSUP INDLI SEGSUP TAB C ON CONNAIT NBMAX, ON PEUT INITIALISER LES SEGMENTS DE TRAVAIL c INITIALISATION DES MATRICES c NBMAX = 10 c NBMAX = NBMAX + 1 c WRITE(6,*) 'NBMAX= ',NBMAX K3 = NSOMM SEGINI INDLI SEGINI TAB DO I = 1,K3 INDLI.ID(I) = 0 TAB.ID(I) = 0 ENDDO K1 = NBMAX K2 = NSOMM SEGINI IND2 SEGINI IND SEGINI IND22 SEGINI VAL1 SEGINI VAL2 SEGINI SCMB C INITIALISATION DU POINTEUR MATRICE2 K3 = NSOMM SEGINI IPO2 DO I = 1,K3 K1 = NBMAX K2 = NBMAX + 1 SEGINI MATR1 IPO2.POINT(I) = MATR1 SEGDES MATR1 ENDDO NFAC=MELEFL.NUM(/2) DO NLCF= 1, NFAC, 1 c WRITE(6,*) 'NLCF= ',NLCF 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) 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) AG1 = 0.0D0 AD1 = 0.0D0 AG2 = 0.0D0 AD2 = 0.0D0 PSCAG1 = 0.0D0 PSCAG2 = 0.0D0 PSCAD1 = 0.0D0 PSCAD2 = 0.0D0 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 c WRITE(6,*) 'NLCG= ',NLCG c WRITE(6,*) 'NBF= ',NBFA c WRITE(6,*) 'MELTFA= ',MELTFA c WRITE(6,*) 'DIMENSION1 ',MELTFA.NUM(/1) c WRITE(6,*) 'DIMENSION2 ',MELTFA.NUM(/2) c WRITE(6,*) 'NGAUX ',MELTFA.NUM(/2) c SEGACT MELTFA 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*(NGS1 -1))+1 XS1=MCOORD.XCOOR(ICELL) YS1=MCOORD.XCOOR(ICELL+1) ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=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 c WRITE(6,*) 'NLCF= ',NLCF,'VECXG11= ',VECXG1(1) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG11= ',VECYG1(1) c WRITE(6,*) 'NLCF= ',NLCF,'VECXG12= ',VECXG1(2) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG12= ',VECYG1(2) c WRITE(6,*) 'NGCF= ',NGCF c WRITE(6,*) 'N1= ',N1,'XF= ',XF,'YF= ',YF c WRITE(6,*) 'N1= ',N1,'XG= ',XG,'YG= ',YG ENDIF IF ((NSOM1.EQ.NGS2).OR.(NSOM2.EQ.NGS2)) THEN ICELL=(3*(NGS2 -1))+1 XS2=MCOORD.XCOOR(ICELL) YS2=MCOORD.XCOOR(ICELL+1) ICELL=(3*(N1 -1))+1 XF=MCOORD.XCOOR(ICELL) YF=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 c SEGDES MELTFA MELTFA = IMELEM(NLCD) NBF = NBFACEL(NLCD) c WRITE(6,*) 'NLCD= ',NLCD c WRITE(6,*) 'NBF= ',NBF c WRITE(6,*) 'MELTFA= ',MELTFA c WRITE(6,*) 'DIMENSION1 ',MELTFA.NUM(/1) c WRITE(6,*) 'DIMENSION2 ',MELTFA.NUM(/2) IF (NLCD.GT.NMAI1) THEN NDAUX = NLCD -NMAI1 ELSE NDAUX = NLCD ENDIF c SEGACT MELTFA 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 c SEGDES MELTFA AG1=0.5D0*ABS( ( (VECXG1(1)*VECYG1(2)) - & (VECYG1(1))*VECXG1(2)) ) AG2=0.5D0*ABS( ( (VECXG2(1)*VECYG2(2)) - & (VECYG2(1))*VECXG2(2)) ) AD1=0.5D0*ABS( ( (VECXD1(1)*VECYD1(2)) - & (VECYD1(1))*VECXD1(2)) ) AD2=0.5D0*ABS( ( (VECXD2(1)*VECYD2(2)) - & (VECYD2(1))*VECXD2(2)) ) c WRITE(6,*) 'NLCF=',NLCF c WRITE(6,*) 'NLCD=',NLCD c WRITE(6,*) 'NLCG=',NLCG c WRite(6,*) 'AG1=',AG1 c WRite(6,*) 'AG2=',AG2 c WRite(6,*) 'AD1=',AD1 c WRite(6,*) 'AD2=',AD2 c WRITE(6,*) 'NLCF= ',NLCF,'VECXG11= ',VECXG1(1) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG11= ',VECYG1(1) c WRITE(6,*) 'NLCF= ',NLCF,'VECXG12= ',VECXG1(2) c WRITE(6,*) 'NLCF= ',NLCF,'VECYG12= ',VECYG1(2) c WRite(6,*) 'PSCAG1=',PSCAG1 c WRite(6,*) 'PSCAG2=',PSCAG2 c WRite(6,*) 'PSCAD1=',PSCAD1 c WRite(6,*) 'PSCAD2=',PSCAD2 c WRite(6,*) 'COEF1D=',COEF1D c WRite(6,*) 'COEF2D=',COEF2D c WRite(6,*) 'BETA1GD=',BETA1GD c WRite(6,*) 'BETA2GD=',BETA2GD c WRite(6,*) 'INDD2=',INDD2 c CALCUL DE MATRICE POUR LE NOEUD D INDICE NS1 IF (ICHTE.EQ.0) THEN COEF1 = ( (VECXG1(INDG1)*SCN1X) + (VECYG1(INDG1)*SCN1Y) ) & / AG1 IAUX = 3 - INDG1 COEF2 = ( (VECXG1(IAUX)*SCN1X) + (VECYG1(IAUX)*SCN1Y) ) & / AG1 COEF3 = ( (VECXD1(INDD1)*SCN1X) + (VECYD1(INDD1)*SCN1Y) ) & / AD1 IAUX = 3 - INDD1 COEF4 = ( (VECXD1(IAUX)*SCN1X) + (VECYD1(IAUX)*SCN1Y) ) & / AD1 ELSE c WRITE(6,*) 'NLCG= ',NLCG,'NLCG2= ',NLCG2 c WRITE(6,*) 'NLCD= ',NLCD,'NLCD2= ',NLCD2 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 c LE TENSEUR EST DIAGONAL 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' STOP ENDIF c xmink11 = min(K11G,xmink11) c xmink11 = min(K11D,xmink11) c xmaxk11 = max(K11G,xmaxk11) c xmaxk11 = max(K11D,xmaxk11) c xmink22 = min(K22G,xmink22) c xmink22 = min(K22D,xmink22) c xmaxk22 = max(K22G,xmaxk22) c xmaxk22 = max(K22D,xmaxk22) c WRITE(6,*) 'NLCF= ',NLCF c WRITE(6,*) 'NLCG= ',NLCG, 'NLCD= ',NLCD c WRite(6,*) 'K11G=',K11G,'K22G= ',K22G,'K21G=',K21G c WRite(6,*) 'K11D=',K11D,'K22D= ',K22D,'K21D=',K21D c ON EST ICI c PRODUIT TENSEUR VECTEUR IAUX = 3 - INDD1 XLONGD = (VECXD1(IAUX)*VECXD1(IAUX)) + & (VECYD1(IAUX)*VECYD1(IAUX)) XLONGD = XLONGD**0.5 IAUX = 3 - INDG1 XLONGG = (VECXG1(IAUX)*VECXG1(IAUX)) + & (VECYG1(IAUX)*VECYG1(IAUX)) XLONGG = XLONGG**0.5 PSCAGX = (K11G*(VECXG1(INDG1)/AG1)) + (K21G*(VECYG1(INDG1)/AG1)) PSCAGY = (K21G*(VECXG1(INDG1))/AG1) + (K22G*(VECYG1(INDG1)/AG1)) COEF1 = ( (PSCAGX*SCN1X) + (PSCAGY*SCN1Y) ) IAUX = 3 - INDG1 PSCAGX = (K11G*(VECXG1(IAUX)/AG1)) + (K21G*(VECYG1(IAUX)/AG1)) PSCAGY = (K21G*(VECXG1(IAUX)/AG1)) + (K22G*(VECYG1(IAUX)/AG1)) COEF2 = ( (PSCAGX*SCN1X) + (PSCAGY*SCN1Y) ) PSCADX = (K11D*(VECXD1(INDD1)/AD1)) + (K21D*(VECYD1(INDD1)/AD1)) PSCADY = (K21D*(VECXD1(INDD1)/AD1)) + (K22D*(VECYD1(INDD1)/AD1)) COEF3 = ( (PSCADX*SCN1X) + (PSCADY*SCN1Y) ) IAUX = 3 - INDD1 PSCADX = (K11D*(VECXD1(IAUX)/AD1)) + (K21D*(VECYD1(IAUX)/AD1)) PSCADY = (K21D*(VECXD1(IAUX)/AD1)) + (K22D*(VECYD1(IAUX)/AD1)) COEF4 = ( (PSCADX*SCN1X) + (PSCADY*SCN1Y) ) ENDIF c WRite(6,*) 'COEF1=',COEF1 c WRite(6,*) 'COEF2=',COEF2 c WRite(6,*) 'COEF3=',COEF3 c WRite(6,*) 'COEF4=',COEF4 MARQ = 0 DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.NGCF) THEN MARQ = 1 IAFF = I5 GOTO 4 ENDIF ENDDO 4 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS1) = INDLI.ID(NLS1) + 1 ICOU = INDLI.ID(NLS1) IND2.NUME(ICOU,NLS1) = NGCF ELSE ICOU = IAFF ENDIF COEF = (COEF1 - COEF3) MATR1 = IPO2.POINT(NLS1) SEGACT MATR1 *MOD MATR1.MAT2(ICOU,ICOU) = COEF MARQ = 0 DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.NG1) THEN MARQ = 1 IAFF = I5 GOTO 5 ENDIF ENDDO 5 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS1) = INDLI.ID(NLS1) + 1 ICOUCO = INDLI.ID(NLS1) IND2.NUME(ICOUCO,NLS1) = NG1 ELSE ICOUCO = IAFF ENDIF ICOUG = ICOUCO MATR1.MAT2(ICOU,ICOUCO) = COEF2 MARQ = 0 DO I5 = 1,INDLI.ID(NLS1) INDAUX = IND2.NUME(I5,NLS1) IF (INDAUX.EQ.ND1) THEN MARQ = 1 IAFF = I5 GOTO 6 ENDIF ENDDO 6 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS1) = INDLI.ID(NLS1) + 1 ICOUCO = INDLI.ID(NLS1) IND2.NUME(ICOUCO,NLS1) = ND1 ELSE ICOUCO = IAFF ENDIF ICOUD = ICOUCO MATR1.MAT2(ICOU,ICOUCO) = -COEF4 SCMB.MAT(ICOU,NLS1) = & (((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1)) - & ((COEF3+COEF4)*MPOCHP.VPOCHA(NLCD,1))) * COEF POUR INVERSER LA MATRICE * ON CORRIGE ICI VAL1.MAT(ICOU,NLS1) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS1) = - ((COEF3 + COEF4)) IND.NUME(ICOU,NLS1) = NGCG IND22.NUME(ICOU,NLS1) = NGCD * CONDITION AUX LIMITE DE DIRICICHLET IF (NGCG.EQ.NGCD) THEN NLFCL=MLENCL.LECT(NGCF) NL1=MLENCL.LECT(NGS1) NL2=MLENCL.LECT(NGS2) IF ((NL1.GT.0).AND.(NL2.GT.0)) THEN COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = 0.0D0 MATR1.MAT2(ICOU,ICOUD) = 0.0D0 IF (NBF.EQ.3) THEN SCMB.MAT(ICOU,NLS1) = c (COEF*2.D0*MPOVCL.VPOCHA(NL1,1)/3.D0) + c (COEF*MPOVCL.VPOCHA(NL2,1)/3.D0) ELSE SCMB.MAT(ICOU,NLS1) = c (COEF*MPOVCL.VPOCHA(NL1,1)/2.D0) + c (COEF*MPOVCL.VPOCHA(NL2,1)/2.D0) ENDIF VAL1.MAT(ICOU,NLS1) = COEF VAL2.MAT(ICOU,NLS1) = 0.0D0 c ON AJOUTE ICI UN POINT FACE POUR COMPATIBILITE AVEC LAPN IND.NUME(ICOU,NLS1) = NGCF IND22.NUME(ICOU,NLS1) = NGCD ELSE NLFNE=MLENNE.LECT(NGCF) c CONDITION DE FLUX IF (NLFNE.GT.0) THEN QIMPX = MPOVNE.VPOCHA(NLFNE,1) C PRODUIT SCALAIRE DU FLUX IMPOSE AVEC LA NORMALE QIMPS = (QIMPX) COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = COEF2 SCMB.MAT(ICOU,NLS1) = ( & ((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) + (2.D0*QIMPS) VAL1.MAT(ICOU,NLS1) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS1) = 2.D0 IND.NUME(ICOU,NLS1) = NGCG IND22.NUME(ICOU,NLS1) = NGCF ELSE c CONDITION MIXTE NLFMI=MLENMI.LECT(NGCF) IF (NLFMI.GT.0) THEN XLAMBDA1 = MPOVMI.VPOCHA(NLFMI,1) XLAMBDA2 = MPOVMI.VPOCHA(NLFMI,2) QIMPX = MPOVMI.VPOCHA(NLFMI,3) C PRODUIT SCALAIRE DU FLUX IMPOSE AVEC LA NORMALE QIMPS = (QIMPX) c WRITE(6,*) 'NLCF= ',NLCF c WRITE(6,*) 'NGCF= ',NGCF c WRITE(6,*) 'XLAMBDA1= ',XLAMBDA1,'XLAMBDA2= ',XLAMBDA2 c WRITE(6,*) 'QIMPX= ',QIMPX,'QIMPY= ',QIMPY COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = (XLAMBDA1*COEF) - & (2.D0*XLAMBDA2) MATR1.MAT2(ICOU,ICOUG) = (XLAMBDA1*COEF2) SCMB.MAT(ICOU,NLS1) = & (XLAMBDA1*((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) & + (2.D0*QIMPS) VAL1.MAT(ICOU,NLS1) = XLAMBDA1*(COEF1 + COEF2) VAL2.MAT(ICOU,NLS1) = 2.D0 IND.NUME(ICOU,NLS1) = NGCG IND22.NUME(ICOU,NLS1) = NGCF ELSE C PAR DEFAUT FLUX NUL QIMPS = 0 COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = COEF2 SCMB.MAT(ICOU,NLS1) = & (((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) VAL1.MAT(ICOU,NLS1) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS1) = 0.D0 IND.NUME(ICOU,NLS1) = NGCG IND22.NUME(ICOU,NLS1) = NGCD ENDIF ENDIF ENDIF ENDIF c WRITE(6,*) 'NLS1= ',NLS1,'ICOU=',ICOU, c & 'IND= ',IND.NUME(ICOU,NLS1), c & 'IND22= ',IND22.NUME(ICOU,NLS1), c & 'SCMB', SCMB.MAT(ICOU,NLS1), c & 'ICOU =',ICOU,'NOUED2= ',MATR1.MAT2(ICOU,ICOU,NLS1), c & 'ICOUG= ',ICOUG,'NOUED2= ',MATR1.MAT2(ICOU,ICOUG,NLS1), c & 'ICOUD= ',ICOUD,'NOUED2= ',MATR1.MAT2(ICOU,ICOUD,NLS1) c WRITE(6,*) 'COEF1 = ',COEF1,'COEF2= ',COEF2,'COEF3= ', c & COEF3,'COEF4=',COEF4,'HG=',MPOCHP.VPOCHA(NLCG,1), c & 'HD= ',MPOCHP.VPOCHA(NLCD,1) SEGDES MATR1 *MOD c CALCUL DE MATRICE POUR LE NOEUD D INDICE NS2 IF (ICHTE.EQ.0) THEN COEF1 = ( (VECXG2(INDG2)*SCN1X) + (VECYG2(INDG2)*SCN1Y) ) & / AG2 IAUX = 3 - INDG2 COEF2 = ( (VECXG2(IAUX)*SCN1X) + (VECYG2(IAUX)*SCN1Y) ) & / AG2 COEF3 = ( (VECXD2(INDD2)*SCN1X) + (VECYD2(INDD2)*SCN1Y) ) & / AD2 IAUX = 3 - INDD2 COEF4 = ( (VECXD2(IAUX)*SCN1X) + (VECYD2(IAUX)*SCN1Y) ) & / AD2 ELSE 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 c LE TENSEUR EST DIAGONAL 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' STOP ENDIF c PRODUIT TENSEUR VECTEUR IAUX = 3 - INDD1 XLONGD = (VECXD1(IAUX)*VECXD1(IAUX)) + & (VECYD1(IAUX)*VECYD1(IAUX)) XLONGD = XLONGD**0.5 IAUX = 3 - INDG1 XLONGG = (VECXG1(IAUX)*VECXG1(IAUX)) + & (VECYG1(IAUX)*VECYG1(IAUX)) XLONGG = XLONGG**0.5 PSCAGX = (K11G*(VECXG2(INDG2)/AG2)) + (K21G*(VECYG2(INDG2)/AG2)) PSCAGY = (K21G*(VECXG2(INDG2))/AG2) + (K22G*(VECYG2(INDG2)/AG2)) COEF1 = ( (PSCAGX*SCN1X) + (PSCAGY*SCN1Y) ) IAUX = 3 - INDG2 PSCAGX = (K11G*(VECXG2(IAUX)/AG2)) + (K21G*(VECYG2(IAUX)/AG2)) PSCAGY = (K21G*(VECXG2(IAUX)/AG2)) + (K22G*(VECYG2(IAUX)/AG2)) COEF2 = ( (PSCAGX*SCN1X) + (PSCAGY*SCN1Y) ) PSCADX = (K11D*(VECXD2(INDD2)/AD2)) + (K21D*(VECYD2(INDD2)/AD2)) PSCADY = (K21D*(VECXD2(INDD2)/AD2)) + (K22D*(VECYD2(INDD2)/AD2)) COEF3 = ( (PSCADX*SCN1X) + (PSCADY*SCN1Y) ) IAUX = 3 - INDD2 PSCADX = (K11D*(VECXD2(IAUX)/AD2)) + (K21D*(VECYD2(IAUX)/AD2)) PSCADY = (K21D*(VECXD2(IAUX)/AD2)) + (K22D*(VECYD2(IAUX)/AD2)) COEF4 = ( (PSCADX*SCN1X) + (PSCADY*SCN1Y) ) ENDIF MARQ = 0 DO I5 = 1,INDLI.ID(NLS2) INDAUX = IND2.NUME(I5,NLS2) IF (INDAUX.EQ.NGCF) THEN MARQ = 1 IAFF = I5 GOTO 41 ENDIF ENDDO 41 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS2) = INDLI.ID(NLS2) + 1 ICOU = INDLI.ID(NLS2) IND2.NUME(ICOU,NLS2) = NGCF ELSE ICOU = IAFF ENDIF COEF = (COEF1 - COEF3) MATR1 = IPO2.POINT(NLS2) SEGACT MATR1 *MOD MATR1.MAT2(ICOU,ICOU) = COEF MARQ = 0 DO I5 = 1,INDLI.ID(NLS2) INDAUX = IND2.NUME(I5,NLS2) IF (INDAUX.EQ.NG2) THEN MARQ = 1 IAFF = I5 GOTO 51 ENDIF ENDDO 51 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS2) = INDLI.ID(NLS2) + 1 ICOUCO = INDLI.ID(NLS2) IND2.NUME(ICOUCO,NLS2) = NG2 ELSE ICOUCO = IAFF ENDIF ICOUG = ICOUCO MATR1.MAT2(ICOU,ICOUCO) = COEF2 MARQ = 0 DO I5 = 1,INDLI.ID(NLS2) INDAUX = IND2.NUME(I5,NLS2) IF (INDAUX.EQ.ND2) THEN MARQ = 1 IAFF = I5 GOTO 61 ENDIF ENDDO 61 CONTINUE IF (MARQ.EQ.0) THEN INDLI.ID(NLS2) = INDLI.ID(NLS2) + 1 ICOUCO = INDLI.ID(NLS2) IND2.NUME(ICOUCO,NLS2) = ND2 ELSE ICOUCO = IAFF ENDIF ICOUD = ICOUCO MATR1.MAT2(ICOU,ICOUCO) = -COEF4 SCMB.MAT(ICOU,NLS2) =( & ((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1)) - & ((COEF3+COEF4)*MPOCHP.VPOCHA(NLCD,1))) VAL1.MAT(ICOU,NLS2) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS2) = - ((COEF3 + COEF4)) IND.NUME(ICOU,NLS2) = NGCG IND22.NUME(ICOU,NLS2) = NGCD * CONDITION AUX LIMITE DE DIRICICHLET IF (NGCG.EQ.NGCD) THEN NLFCL=MLENCL.LECT(NGCF) NL1=MLENCL.LECT(NGS1) NL2=MLENCL.LECT(NGS2) IF ((NL1.GT.0).AND.(NL2.GT.0)) THEN c WRITE(6,*) 'NLCF= ',NLCF,'NGCF= ',NGCF c WRITE(6,*) 'CLIM= ', MPOVCL.VPOCHA(NLFCL,1) COEF = MAX(ABS(COEF1),ABS(COEF2)) c WRITE(6,*) 'COEF= ',COEF c WRITE(6,*) 'COEF1= ',COEF1 c WRITE(6,*) 'COEF2= ',COEF2 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = 0.0D0 MATR1.MAT2(ICOU,ICOUD) = 0.0D0 c SCMB.MAT(ICOU,NLS2) = (COEF*MPOVCL.VPOCHA(NLFCL,1)) IF (NBF.EQ.3) THEN SCMB.MAT(ICOU,NLS2) = c (COEF*2.D0*MPOVCL.VPOCHA(NL2,1)/3.D0) + c (COEF*MPOVCL.VPOCHA(NL1,1)/3.D0) ELSE SCMB.MAT(ICOU,NLS2) = c (COEF*MPOVCL.VPOCHA(NL2,1)/2.D0) + c (COEF*MPOVCL.VPOCHA(NL1,1)/2.D0) ENDIF VAL1.MAT(ICOU,NLS2) = COEF VAL2.MAT(ICOU,NLS2) = 0.0D0 c ON AJOUTE ICI UN POINT FACE POUR COMPATIBILITE AVEC LAPN IND.NUME(ICOU,NLS2) = NGCF IND22.NUME(ICOU,NLS2) = NGCD ELSE c CONDITION DE FLUX NLFNE=MLENNE.LECT(NGCF) IF (NLFNE.GT.0) THEN QIMPX = MPOVNE.VPOCHA(NLFNE,1) C PRODUIT SCALAIRE DU FLUX IMPOSE AVEC LA NORMALE QIMPS = (QIMPX) COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = COEF2 SCMB.MAT(ICOU,NLS2) = & (((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) & + (2.D0*QIMPS) VAL1.MAT(ICOU,NLS2) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS2) = 2.D0 IND.NUME(ICOU,NLS2) = NGCG IND22.NUME(ICOU,NLS2) = NGCF ELSE c CONDITION MIXTE NLFMI=MLENMI.LECT(NGCF) IF (NLFMI.GT.0) THEN XLAMBDA1 = MPOVMI.VPOCHA(NLFMI,1) XLAMBDA2 = MPOVMI.VPOCHA(NLFMI,2) QIMPX = MPOVMI.VPOCHA(NLFMI,3) C PRODUIT SCALAIRE DU FLUX IMPOSE AVEC LA NORMALE QIMPS = (QIMPX) c WRITE(6,*) 'QIMPX= ',QIMPX,'QIMPY= ',QIMPY c COEF = COEF1 c WRITE(6,*) 'NGCF= ',NGCF c WRITE(6,*) 'XLAMBDA1= ',XLAMBDA1,'XLAMBDA2= ',XLAMBDA2 c WRITE(6,*) 'COEF= ',COEF MATR1.MAT2(ICOU,ICOU) = (XLAMBDA1*COEF) - & (2.D0*XLAMBDA2) MATR1.MAT2(ICOU,ICOUG) = (XLAMBDA1*COEF2) SCMB.MAT(ICOU,NLS2) = & (XLAMBDA1*((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) & + (2.D0*QIMPS) VAL1.MAT(ICOU,NLS2) = XLAMBDA1*(COEF1 + COEF2) VAL2.MAT(ICOU,NLS2) = 2.D0 IND.NUME(ICOU,NLS2) = NGCG IND22.NUME(ICOU,NLS2) = NGCF ELSE C PAR DEFAUT FLUX NUL QIMPS = 0 COEF = COEF1 MATR1.MAT2(ICOU,ICOU) = COEF MATR1.MAT2(ICOU,ICOUG) = COEF2 SCMB.MAT(ICOU,NLS2) = & (((COEF1+COEF2)*MPOCHP.VPOCHA(NLCG,1))) VAL1.MAT(ICOU,NLS2) = (COEF1 + COEF2) VAL2.MAT(ICOU,NLS2) = 0.D0 IND.NUME(ICOU,NLS2) = NGCG IND22.NUME(ICOU,NLS2) = NGCD ENDIF ENDIF ENDIF ENDIF SEGDES MATR1 *MOD c WRITE(6,*) 'NLS2= ',NLS2,'ICOU=',ICOU, c & 'IND= ',IND.NUME(ICOU,NLS2), c & 'IND22= ',IND22.NUME(ICOU,NLS2) c WRITE(6,*) 'NLS2= ',NLS2,'ICOU=',ICOU,'SCMB', SCMB.MAT(NLS2,ICOU) c WRITE(6,*) 'COEF1 = ',COEF1,'COEF2= ',COEF2,'COEF3= ', c & COEF3,'COEF4=',COEF4,'HG=',MPOCHP.VPOCHA(NLCG,1), c & 'HD= ',MPOCHP.VPOCHA(NLCD,1) c WRITE(6,*) 'NLS2= ',NLS2,'ICOU=',ICOU, c & 'IND= ',IND.NUME(ICOU,NLS2), c & 'IND22= ',IND22.NUME(ICOU,NLS2), c & 'SCMB', SCMB.MAT(ICOU,NLS2), c & 'ICOU =',ICOU,'NOUED2= ',MATR1.MAT2(ICOU,ICOU,NLS2), c & 'ICOUG= ',ICOUG,'NOUED2= ',MATR1.MAT2(ICOU,ICOUG,NLS2), c & 'ICOUD= ',ICOUD,'NOUED2= ',MATR1.MAT2(ICOU,ICOUD,NLS2) c WRITE(6,*) 'COEF1 = ',COEF1,'COEF2= ',COEF2,'COEF3= ', c DO I=1,INDLI.ID(NLS1) c DO J = 1,INDLI.ID(NLS1) c WRITE(6,*) 'NLS1= ',NLS1,'I=',I,'J=',J,MATR1.MAT2(I,J,NLS1) c WRITE(6,*) 'NLS2= ',NLS2,'I=',I,'J=',J,MATR1.MAT2(I,J,NLS2) c ENDDO c ENDDO ENDDO c DO J= 1,INDLI.ID(NLS1) c WRITE(6,*) 'MELVA1=',MELVA1.VELCHE(J,NLCF) c WRITE(6,*) 'MELVA2=',MELVA1.VELCHE(J,NLCF) c WRITE(6,*) 'MELEME=',MELEME.NUM(J,NLCF) c ENDDO MELTFA = MAUX IF (NBSO.EQ.2) THEN SEGDES IPT1 SEGDES IPT2 ENDIF 9999 CONTINUE RETURN END