C PB902 SOURCE MAGN 10/06/02 21:15:00 6681 SUBROUTINE PB902 &(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NG,NPG,NOM2,ITYPI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C************************************************************************ C C CALCULE LES FONCTIONS DE FORME D'UN : QUA9 C eta C ^ C | C 7 6 5 C 1. x------x------x C | | C | | C alfa 8x 9x 4x C | | C | | C 0. x------x------x--> ksi C 1 2 3 C C 0. alfa 1. C C ITYPI=2 pts d'intégration de type Gauss Lobatto i.e. sommets élément C C REMARQUE : Pour l'élément quadratique. C En coordonnées cylindriques si les points milieux (2,4,6,8,9) sont C pile poil au milieu, la matrice masse diagonale s'annule sur l'axe. C /1 C | 2 pi x N1 dx = 0. C /0 C Pour remédier à cela on décale un peu les points mileux sur l'élément C de référence, alfa=0.50005 au lieu de 0.5 . C Si on rapprochait les points milieux de l'axe (alfa = 0.45 par exemple) C la matrice masse diagonale serait négative sur l'axe, ce que l'on ne C veut pas non plus. En conséquence le choix qu'on a fait doit être C accompagné d'une numérotation dans le sens trigo de l'élément courant. C************************************************************************ REAL*8 XREF(ND,NP),X(NPG),Y(NPG) DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) DIMENSION FM(MP,NPG),GM(ND,MP,NPG) CHARACTER*4 NOM2 REAL*8 A,B,C,D,U(5),H(5),AUX1,AUX2 -INC PPARAM -INC CCOPTIO C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> AUX1=1.D0/2.D0+3.D0**(1.D0/2.D0)/6.D0 AUX2=1.D0/2.D0-3.D0**(1.D0/2.D0)/6.D0 C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF(IFOMOD.EQ.0)THEN alfa=0.50005D0 ELSE alfa=0.5D0 ENDIF alfai=1.D0/alfa alfa2=1.D0/alfa/(alfa - 1.D0) alfa3=1.D0/(alfai - 1.D0) XREF(1,1)=0.D0 XREF(2,1)=0.D0 XREF(1,2)=alfa XREF(2,2)=0.D0 XREF(1,3)=1.D0 XREF(2,3)=0.D0 XREF(1,4)=1.D0 XREF(2,4)=alfa XREF(1,5)=1.D0 XREF(2,5)=1.D0 XREF(1,6)=alfa XREF(2,6)=1.D0 XREF(1,7)=0.D0 XREF(2,7)=1.D0 XREF(1,8)=0.D0 XREF(2,8)=alfa XREF(1,9)=alfa XREF(2,9)=alfa C Verification des coordonnées c IF(.TRUE.)THEN IF(.FALSE.)THEN DO 11 L=1,NP X(L)=XREF(1,L) Y(L)=XREF(2,L) 11 CONTINUE DO 12 L=1,NP FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) FN(2,L)=alfa2*X(L)*(X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) FN(6,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) FN(9,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L) & ,FN(7,L),FN(8,L),FN(9,L) 12 CONTINUE 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4)) ENDIF C Fin Vérification CALL CALUHG(U,H,NG) C A=0.D0 B=1.D0 C=0.D0 D=1.D0 IF(ITYPI.EQ.2)THEN X(1)=0.D0 Y(1)=0.D0 X(2)=alfa Y(2)=0.D0 X(3)=1.D0 Y(3)=0.D0 X(4)=1.D0 Y(4)=alfa X(5)=1.D0 Y(5)=1.D0 X(6)=alfa Y(6)=1.D0 X(7)=0.D0 Y(7)=1.D0 X(8)=0.D0 Y(8)=alfa X(9)=alfa Y(9)=alfa DO 2 L=1,NPG PG(L)=1.D0/9.D0 2 CONTINUE ELSE CALL CALG2(A,B,C,D,NG,H,U,X,Y,PG) ENDIF DO 1 L=1,NPG C IF(NOM2.EQ.'PRP0')THEN FM(1,L)= 1.D0 GM(1,1,L)=0.D0 GM(2,1,L)=0.D0 ELSEIF(NOM2.EQ.'PRQ1')THEN FM(1,L)= 3.D0*(X(L)-AUX1)*(Y(L)-AUX1) FM(2,L)=-3.D0*(X(L)-AUX2)*(Y(L)-AUX1) FM(3,L)= 3.D0*(X(L)-AUX2)*(Y(L)-AUX2) FM(4,L)=-3.D0*(X(L)-AUX1)*(Y(L)-AUX2) C GM(1,1,L)=3.D0*(Y(L)-AUX1) GM(2,1,L)=3.D0*(X(L)-AUX1) GM(1,2,L)=-3.D0*(Y(L)-AUX1) GM(2,2,L)=-3.D0*(X(L)-AUX2) GM(1,3,L)=3.D0*(Y(L)-AUX2) GM(2,3,L)=3.D0*(X(L)-AUX2) GM(1,4,L)=-3.D0*(Y(L)-AUX2) GM(2,4,L)=-3.D0*(X(L)-AUX1) ELSEIF(NOM2.EQ.'PRP1')THEN C? FM(1,L)= (3.D0-4.D0*X(L))/2.D0 C? FM(2,L)= 2.D0*(X(L)-Y(L)) C? FM(3,L)= (4.D0*Y(L)-1.D0)/2.D0 FM(1,L)= (4.D0*X(L)+4.D0*Y(L)-3.D0)/3.D0 FM(2,L)=-1.D0/3.D0*(8.D0*X(L)-4.D0*Y(L)-3.D0) FM(3,L)= 1.D0/3.D0*(-8.D0*Y(L)+4.D0*X(L)+3.D0) C GM(1,1,L)=0.D0 GM(2,1,L)=0.D0 GM(1,2,L)=0.D0 GM(2,2,L)=0.D0 GM(1,3,L)=0.D0 GM(2,3,L)=0.D0 ELSEIF(NOM2.EQ.'PFP1')THEN FM(1,L)= (X(L)-1.D0)*(Y(L)-1.D0) FM(2,L)=-X(L)*(Y(L)-1.D0) FM(3,L)= X(L)*Y(L) FM(4,L)=-Y(L)*(X(L)-1.D0) GM(1,1,L)=Y(L)-1.D0 GM(2,1,L)=X(L)-1.D0 GM(1,2,L)=-(Y(L)-1.D0) GM(2,2,L)=-X(L) GM(1,3,L)=Y(L) GM(2,3,L)=X(L) GM(1,4,L)=-Y(L) GM(2,4,L)=-(X(L)-1.D0) ENDIF C c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0) c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c c N2X=alfa2*X(L)*(X(L)-1.D0) c N2Y=alfa2*Y(L)*(Y(L)-1.D0) c c N3X=alfa3*X(L)*(alfai*X(L)-1.D0) c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0) C c FN(1,L)=N1X*N1Y FN(1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c FN(2,L)=N2X*N1Y FN(2,L)=alfa2*X(L)*(X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c FN(3,L)=N3X*N1Y FN(3,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c FN(4,L)=N3X*N2Y FN(4,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) c FN(5,L)=N3X*N3Y FN(5,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c FN(6,L)=N2X*N3Y FN(6,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c FN(7,L)=N1X*N3Y FN(7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c FN(8,L)=N1X*N2Y FN(8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) c FN(9,L)=N2X*N2Y FN(9,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) C c FN(1,L)=(X(L)-1.D0)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0) c FN(2,L)=-4.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0) c FN(3,L)=X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0) c FN(4,L)=-4.D0*X(L)*Y(L)*(2.D0*X(L)-1.D0)*(Y(L)-1.D0) c FN(5,L)=X(L)*Y(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0) c FN(6,L)=-4.D0*X(L)*Y(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0) c FN(7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0)*(X(L)-1.D0) c FN(8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(X(L)-1.D0)*(2.D0*X(L)-1.D0) c FN(9,L)=16.D0*X(L)*Y(L)*(X(L)-1.D0)*(Y(L)-1.D0) C c N1X=(X(L)-1.D0)*(alfai*X(L)-1.D0) c N1Y=(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c c N2X=alfa2*X(L)*(X(L)-1.D0) c N2Y=alfa2*Y(L)*(Y(L)-1.D0) c c N3X=alfa3*X(L)*(alfai*X(L)-1.D0) c N3Y=alfa3*Y(L)*(alfai*Y(L)-1.D0) c c G1X=2.D0*alfai*X(L)-(1.D0 + alfai) c G1Y=2.D0*alfai*Y(L)-(1.D0 + alfai) c c G2X=alfa2*(2.D0*X(L)-1.D0) c G2Y=alfa2*(2.D0*Y(L)-1.D0) c c G3X=alfa3*(2.D0*alfai*X(L)-1.D0) c G3Y=alfa3*(2.D0*alfai*Y(L)-1.D0) C c GR(1,1,L)=G1X*N1Y GR(1,1,L)=(2.D0*alfai*X(L)-(1.D0 + alfai)) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c GR(2,1,L)=N1X*G1Y GR(2,1,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *(2.D0*alfai*Y(L)-(1.D0 + alfai)) c GR(1,2,L)=G2X*N1Y GR(1,2,L)=(Y(L)-1.D0)*(alfai*Y(L)-1.D0) & *alfa2*(2.D0*X(L)-1.D0) c GR(2,2,L)=N2X*G1Y GR(2,2,L)=alfa2*X(L)*(X(L)-1.D0) & *(2.D0*alfai*Y(L)-(1.D0 + alfai)) c GR(1,3,L)=G3X*N1Y GR(1,3,L)=alfa3*(2.D0*alfai*X(L)-1.D0) & *(Y(L)-1.D0)*(alfai*Y(L)-1.D0) c GR(2,3,L)=N3X*G1Y GR(2,3,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *(2.D0*alfai*Y(L)-(1.D0 + alfai)) c GR(1,4,L)=G3X*N2Y GR(1,4,L)=alfa3*(2.D0*alfai*X(L)-1.D0) & *alfa2*Y(L)*(Y(L)-1.D0) c GR(2,4,L)=N3X*G2Y GR(2,4,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa2*(2.D0*Y(L)-1.D0) c GR(1,5,L)=G3X*N3Y GR(1,5,L)=alfa3*(2.D0*alfai*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c GR(2,5,L)=N3X*G3Y GR(2,5,L)=alfa3*X(L)*(alfai*X(L)-1.D0) & *alfa3*(2.D0*alfai*Y(L)-1.D0) c GR(1,6,L)=G2X*N3Y GR(1,6,L)=alfa2*(2.D0*X(L)-1.D0) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c GR(2,6,L)=N2X*G3Y GR(2,6,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa3*(2.D0*alfai*Y(L)-1.D0) c GR(1,7,L)=G1X*N3Y GR(1,7,L)=(2.D0*alfai*X(L)-(1.D0 + alfai)) & *alfa3*Y(L)*(alfai*Y(L)-1.D0) c GR(2,7,L)=N1X*G3Y GR(2,7,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa3*(2.D0*alfai*Y(L)-1.D0) c GR(1,8,L)=G1X*N2Y GR(1,8,L)=(2.D0*alfai*X(L)-(1.D0 + alfai)) & *alfa2*Y(L)*(Y(L)-1.D0) c GR(2,8,L)=N1X*G2Y GR(2,8,L)=(X(L)-1.D0)*(alfai*X(L)-1.D0) & *alfa2*(2.D0*Y(L)-1.D0) c GR(1,9,L)=G2X*N2Y GR(1,9,L)=alfa2*Y(L)*(Y(L)-1.D0) & *alfa2*(2.D0*X(L)-1.D0) c GR(2,9,L)=N2X*G2Y GR(2,9,L)=alfa2*X(L)*(X(L)-1.D0) & *alfa2*(2.D0*Y(L)-1.D0) C c GR(1,1,L)=(Y(L)-1.D0)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0) c GR(2,1,L)=(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0) c GR(1,2,L)=-4.D0*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0)*(Y(L)-1.D0) c GR(2,2,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-3.D0) c GR(1,3,L)=(2.D0*Y(L)-1.D0)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0) c GR(2,3,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-3.D0) c GR(1,4,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-1.D0) c GR(2,4,L)=-4.D0*X(L)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0) c GR(1,5,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-1.D0) c GR(2,5,L)=X(L)*(2.D0*X(L)-1.D0)*(4.D0*Y(L)-1.D0) c GR(1,6,L)=-4.D0*Y(L)*(2.D0*Y(L)-1.D0)*(2.D0*X(L)-1.D0) c GR(2,6,L)=-4.D0*X(L)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0) c GR(1,7,L)=Y(L)*(2.D0*Y(L)-1.D0)*(4.D0*X(L)-3.D0) c GR(2,7,L)=(2.D0*X(L)-1.D0)*(X(L)-1.D0)*(4.D0*Y(L)-1.D0) c GR(1,8,L)=-4.D0*Y(L)*(Y(L)-1.D0)*(4.D0*X(L)-3.D0) c GR(2,8,L)=-4.D0*(X(L)-1.D0)*(2.D0*X(L)-1.D0)*(2.D0*Y(L)-1.D0) c GR(1,9,L)=16.D0*Y(L)*(Y(L)-1.D0)*(2.D0*X(L)-1.D0) c GR(2,9,L)=16.D0*X(L)*(X(L)-1.D0)*(2.D0*Y(L)-1.D0) C 1 CONTINUE C C C WRITE(6,101) C WRITE(6,1002)FM C WRITE(6,1002)GM C WRITE(6,1002)FN C WRITE(6,1002)GR C WRITE(6,101) C write(6,*)' RET PB902 ' RETURN 1002 FORMAT(10(1X,1PD11.4)) 1001 FORMAT(20(1X,I5)) 101 FORMAT(1X,'... SUBPB902 ... FM,GM,FN,GR ',9(10H..........)/) C END