C PB702 SOURCE MAGN 10/05/31 21:15:16 6679 SUBROUTINE PB702(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NPG,NOM2,ITYPI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C CALCULE LES FONCTIONS DE FORME D'UN : TRI7 C C Cf Crouzeix Raviard. Conforming and nonconforming finite element C methods for solving the stationary Stokes equations I C R.A.I.R.O. (7e année, décembre 1973,R-3,p.33 à 76) C C C ^ eta C | C |n3 C r2 5 1 -> A1 C |\ 2 -> A12 - C | \ r2=sqrt(2) 3 -> A2 C | \ 4 -> A23 C | \ 5 -> A3 C | \ 6 -> A13 C | \ 7 -> A123 C r2/2 6------4 C | |\ D1=1-x/r2-y/r2 C | | \ D2=x/r2 C | 7 | \ D3=y/r2 C | | \ C | | \ C 1------2-----3--->ksi C 0 r2/2 r2 C C Pi = Di(aDi-c)+FD1xD2xD3 1 <= i <= 3 C Pij = bxDi(aDi-c)-GxD1xD2xD3 1 <= i < j <= 3 C P123= QxD1xD2xD3 C Pour l'élément de référence on a a=2 b=4 c=1 F=3 G=12 C Si on veut perturber un peu la position des points il faut néanmoins C vérifier les propriétés essentielles des fonctions de forme qui outre C Pi=1 en X=i et Pi=0 ailleur on doit avoir Somme des Pi = 1 C (intégration exacte d'une fonction constante). Ceci conduit à des C relations entre a,b,c,F,G et Q C b=2a c=a-1 Q=3(G-F) (Somme Pi = 1) a ne 2 points Aij non au milieux C des arètes C C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element C C REMARQUE : C En coordonnées cylindriques,pour une distribution régulière des noeuds C conforme à l'article de Crouzeix Raviard, la matrice masse diagonale C s'annule sur l'axe, comme pour le QUA9 (voir pb902.eso). C /1 C | 2 pi x N1 dx = 0. C /0 C Pour remédier à cela on se sert de la bulle en la décalant un peu par C rapport à l'élément de référence en jouant sur les constantes F et G. C On a pas à toucher "a" position des points 2,4 et 6. C Le choix que l'on a fait conduit à une matrice masse diagonale >0. En C conséquence ce choix doit être accompagné d'une numérotation dans le C sens trigo de l'élément courant. C 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) REAL*8 A,B,C,D,U(5),H(5) CHARACTER*4 NOM2 -INC PPARAM -INC CCOPTIO C*** IF(IFOMOD.EQ.0)THEN F=3.0001 G=11.9999 F=3. G=12. ELSE F=3. G=12. ENDIF Q=3.*(G-F) a=2. b=2.*a c=a-1. R2=SQRT(2.D0) IPOSP=2 XREF(1,1)=0.D0 XREF(2,1)=0.D0 XREF(1,2)=R2/2.D0 XREF(2,2)=0.D0 XREF(1,3)=R2 XREF(2,3)=0.D0 XREF(1,4)=R2/2.D0 XREF(2,4)=R2/2.D0 XREF(1,5)=0.D0 XREF(2,5)=R2 XREF(1,6)=0. XREF(2,6)=R2/2.D0 XREF(1,7)=R2/3.D0 XREF(2,7)=R2/3.D0 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 D1=1.D0-X(L)/R2-Y(L)/R2 D2=X(L)/R2 D3=Y(L)/R2 FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3 FN(2,L)=b*D1*D2-G*D1*D2*D3 FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3 FN(4,L)=b*D2*D3-G*D1*D2*D3 FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3 FN(6,L)=b*D1*D3-G*D1*D2*D3 FN(7,L)=Q*D1*D2*D3 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) 12 CONTINUE 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4)) ENDIF C Fin Vérification IF(ITYPI.EQ.2)THEN X(1)=0.D0 Y(1)=0.D0 X(2)=R2/2.D0 Y(2)=0.D0 X(3)=R2 Y(3)=0.D0 X(4)=R2/2.D0 Y(4)=R2/2.D0 X(5)=0.D0 Y(5)=R2 X(6)=0.D0 Y(6)=R2/2.D0 X(7)=(X(1)+X(3)+X(5))/3.D0 Y(7)=(Y(1)+Y(3)+Y(5))/3.D0 DO 2 L=1,7 PG(L)=1.D0/7.D0 2 CONTINUE ELSE CALL CALUHH(X,Y,PG,NPG) 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.'PRP1')THEN C? FM(1,L)=R2*(X(L)+Y(L)-R2/2.D0) C? FM(2,L)=-R2*(X(L)-R2/2.D0) C? FM(3,L)=-R2*(Y(L)-R2/2.D0) A1=6.D0/5.D0/R2 FM(1,L)=5.D0/3.D0*(1.D0-A1*X(L)-A1*Y(L)) FM(2,L)=-1.D0/3.D0*(1.D0-6.D0/R2*X(L)) FM(3,L)=-1.D0/3.D0*(1.D0-6.D0/R2*Y(L)) 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)=-R2/2.D0*(X(L)+Y(L)-R2) FM(2,L)=R2/2.D0*X(L) FM(3,L)=R2/2.D0*Y(L) GM(1,1,L)=-R2/2.D0 GM(2,1,L)=-R2/2.D0 GM(1,2,L)=R2/2.D0 GM(2,2,L)=0.D0 GM(1,3,L)=0.D0 GM(2,3,L)=R2/2.D0 ENDIF C c IF(.FALSE.)THEN IF(.TRUE.)THEN C D1=1.D0-X(L)/R2-Y(L)/R2 D2=X(L)/R2 D3=Y(L)/R2 FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3 FN(2,L)=b*D1*D2-G*D1*D2*D3 FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3 FN(4,L)=b*D2*D3-G*D1*D2*D3 FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3 FN(6,L)=b*D1*D3-G*D1*D2*D3 FN(7,L)=Q*D1*D2*D3 c write(6,*)FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L),FN(7,L) D1x=-1./R2 D1y=-1./R2 D2x=1./R2 D2y=0. D3x=0. D3y=1/R2 GR(1,1,L)=D1x*(a*D1-c)+D1*a*D1x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,1,L)=D1y*(a*D1-c)+D1*a*D1y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,2,L)=b*(D1x*D2+D1*D2x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,2,L)=b*(D1y*D2+D1*D2y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,3,L)=D2x*(a*D2-c)+D2*a*D2x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,3,L)=D2y*(a*D2-c)+D2*a*D2y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,4,L)=b*(D2x*D3+D2*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,4,L)=b*(D2y*D3+D2*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,5,L)=D3x*(a*D3-c)+D3*a*D3x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,5,L)=D3y*(a*D3-c)+D3*a*D3y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,6,L)=b*(D1x*D3+D1*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,6,L)=b*(D1y*D3+D1*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) GR(1,7,L)=Q*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) GR(2,7,L)=Q*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) ELSE c write(6,*)'La base' C C FN(1,L)=(1.D0-X(L)/R2-Y(L)/R2)*(1.D0-2.D0*X(L)/R2-2.D0*Y(L)/R2+ & 3.D0/2.D0*X(L)*Y(L)) FN(2,L)=(1.D0-X(L)/R2-Y(L)/R2)*(4.D0/R2-6.D0*Y(L))*X(L) FN(3,L)=X(L)/R2*(2.D0*X(L)/R2-1.D0+3.D0*(1.D0-X(L)/R2-Y(L)/R2) & *Y(L)/R2) FN(4,L)=2.D0*X(L)*Y(L)*(1.D0-3.D0*(1.D0-X(L)/R2-Y(L)/R2)) FN(5,L)=Y(L)*Y(L)-Y(L)*R2/2.D0+3.D0/2.D0*X(L)*Y(L)- & 3.D0*R2/4.D0*X(L)*X(L)*Y(L)-3.D0*R2/4.D0*X(L)*Y(L)*Y(L) FN(6,L)=2.D0*R2*(-R2/2.D0*Y(L)*Y(L)+3.D0/2.D0*(Y(L)*Y(L)*X(L)+ & X(L)*X(L)*Y(L))-2.D0*R2*X(L)*Y(L)+Y(L)) FN(7,L)=27.D0/2.D0*(X(L)*Y(L)-R2/2.D0*X(L)*X(L)*Y(L)- &R2/2.D0*X(L)*Y(L)*Y(L)) C C C GR(1,1,L)=-3.D0*R2/4.D0*Y(L)*Y(L)+2.D0*X(L)+7./2.D0*Y(L)- * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0 GR(2,1,L)=-3.D0*R2/4.D0*X(L)*X(L)+2.D0*Y(L)+7./2.D0*X(L)- * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0 GR(1,2,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L)-8.D0*Y(L)+ * 2.D0*R2 GR(2,2,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L) GR(1,3,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)+ * 2.D0*X(L)+3.D0/2.D0*Y(L)-R2/2.D0 GR(2,3,L)=R2/2.D0*(-3.D0/2.D0*X(L)*X(L)-3.D0*X(L)*Y(L)+ * 3.D0*R2/2.D0*X(L)) GR(1,4,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*Y(L) GR(2,4,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L) GR(1,5,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L) & +3.D0/2.D0*Y(L) GR(2,5,L)=-3.D0*R2/4.D0*X(L)*X(L)-3.D0*R2/2.D0*X(L)*Y(L) & +2.D0*Y(L)+ * 3.D0/2.D0*X(L)-R2/2.D0 GR(1,6,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-8.D0*Y(L) GR(2,6,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L) & -4.D0*Y(L)+2.D0*R2 GR(1,7,L)=27.D0/2.D0*(-R2/2*Y(L)*Y(L)-R2*X(L)*Y(L)+Y(L)) GR(2,7,L)=27.D0/2.D0*(-R2/2*X(L)*X(L)-R2*X(L)*Y(L)+X(L)) C ENDIF 1 CONTINUE IF(.FALSE.)THEN c IF(.TRUE.)THEN DO 121 L=1,NPG F1=0. DO 131 I=1,NP F1=F1+FN(I,L) 131 CONTINUE Write(6,*)' L=',L,F1 121 CONTINUE DO 122 L=1,NPG G1=0. G2=0. DO 141 I=1,NP G1=G1+GR(1,I,L) G2=G2+GR(2,I,L) 141 CONTINUE Write(6,*)' L=',L,G1,G2 122 CONTINUE ENDIF c stop 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) RETURN 1002 FORMAT(10(1X,1PD11.4)) 1001 FORMAT(20(1X,I5)) 101 FORMAT(1X,'... SUBPB702 ... FM,GM,FN,GR ',9(10H..........)/) C END