C PB443 SOURCE CHAT 05/01/13 02:10:40 5004 SUBROUTINE PB443(X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NGG,NPG,NOM2) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C************************************************************************ C C CALCULE LES FONCTIONS DE FORME D'UN : Iso-Q2 (iso P1/P0 nc) TE10 C C C C ^ zeta C | C |10 C | C | C | 9 C | /\ /^ eta C | / \ / C |/____\ / C |7 8 / C | / C | 5 C | / \ C | / \ C | 6_____\4 C | / \ /\ C | / \ / \ C |/_____\/____\ ____________________>ksi C 1 2 3 C C C************************************************************************ CHARACTER*4 NOM2 REAL*8 AL,BE PARAMETER (AL=.5854101966249684D0) PARAMETER (BE=.1381966011250105D0) REAL*8 X(NPG),Y(NPG),Z(NPG) REAL*8 FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) REAL*8 FM(MP,NPG),GM(ND,MP,NPG) INTEGER CONEK(4,8),L,K,I,N,N1,N2,N3,N4,IGAU,IFN INTEGER ND,MP,NGG,NP,NPG REAL*8 KXYZ(3,10),AA(4) REAL*8 T1X,T1Y,T1Z,C1,P1C REAL*8 T2X,T2Y,T2Z,C2,P2C REAL*8 T3X,T3Y,T3Z,C3,P3C REAL*8 T4X,T4Y,T4Z,C4,P4C C C Donnees de la numerotation locale des 4 noeuds de chaque tetraedre C DATA CONEK/1,2,6,7, & 2,8,6,7, & 2,3,4,8, & 2,4,6,8, & 4,5,6,9, & 8,9,4,6, & 6,7,8,9, & 7,8,9,10/ C C Donnees des coordonnees des points du tetraedre TE10 C X Y Z DATA KXYZ/0.0, 0.0, 0.0, & 0.5, 0.0, 0.0, & 1.0, 0.0, 0.0, & 0.5, 0.5, 0.0, & 0.0, 1.0, 0.0, & 0.0, 0.5, 0.0, & 0.0, 0.0, 0.5, & 0.5, 0.0, 0.5, & 0.0, 0.5, 0.5, & 0.0, 0.0, 1.0/ C C On initialise les fonctions tests et les gradients a 0 C DO 10 L=1,NPG C Poids de Gauss PG(L) = 0.1875D0 DO 22 K=1,MP FM(K,L)=0.0D0 DO 22 I=1,ND GM(I,K,L)=0.0D0 22 CONTINUE DO 20 K=1,NP FN(K,L)=0.0D0 DO 20 I=1,ND GR(I,K,L)=0.0D0 20 CONTINUE 10 CONTINUE C C On traite chaque tetraedre elementaire C DO 40 K=1,8 C C On travaille sur chacun des tetraedres composant le macro C N1=CONEK(1,K) N2=CONEK(2,K) N3=CONEK(3,K) N4=CONEK(4,K) C Calcul des equations de plans oppose a chaque point C C Plan 1 T1X=(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(3,N4)-KXYZ(3,N2)) & -(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(2,N4)-KXYZ(2,N2)) T1Y=(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(1,N4)-KXYZ(1,N2)) & -(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(3,N4)-KXYZ(3,N2)) T1Z=(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(2,N4)-KXYZ(2,N2)) & -(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(1,N4)-KXYZ(1,N2)) C1=-(T1X*KXYZ(1,N2)+T1Y*KXYZ(2,N2)+T1Z*KXYZ(3,N2)) P1C=T1X*KXYZ(1,N1)+T1Y*KXYZ(2,N1)+T1Z*KXYZ(3,N1)+C1 T1X=T1X/P1C T1Y=T1Y/P1C T1Z=T1Z/P1C C1=C1/P1C C Plan 2 T2X=(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(3,N1)-KXYZ(3,N3)) & -(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(2,N1)-KXYZ(2,N3)) T2Y=(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(1,N1)-KXYZ(1,N3)) & -(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(3,N1)-KXYZ(3,N3)) T2Z=(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(2,N1)-KXYZ(2,N3)) & -(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(1,N1)-KXYZ(1,N3)) C2=-(T2X*KXYZ(1,N3)+T2Y*KXYZ(2,N3)+T2Z*KXYZ(3,N3)) P2C=T2X*KXYZ(1,N2)+T2Y*KXYZ(2,N2)+T2Z*KXYZ(3,N2)+C2 T2X=T2X/P2C T2Y=T2Y/P2C T2Z=T2Z/P2C C2=C2/P2C C Plan 3 T3X=(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(3,N2)-KXYZ(3,N4)) & -(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(2,N2)-KXYZ(2,N4)) T3Y=(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(1,N2)-KXYZ(1,N4)) & -(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(3,N2)-KXYZ(3,N4)) T3Z=(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(2,N2)-KXYZ(2,N4)) & -(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(1,N2)-KXYZ(1,N4)) C3=-(T3X*KXYZ(1,N4)+T3Y*KXYZ(2,N4)+T3Z*KXYZ(3,N4)) P3C=T3X*KXYZ(1,N3)+T3Y*KXYZ(2,N3)+T3Z*KXYZ(3,N3)+C3 T3X=T3X/P3C T3Y=T3Y/P3C T3Z=T3Z/P3C C3=C3/P3C C Plan 4 T4X=(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(3,N3)-KXYZ(3,N1)) & -(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(2,N3)-KXYZ(2,N1)) T4Y=(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(1,N3)-KXYZ(1,N1)) & -(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(3,N3)-KXYZ(3,N1)) T4Z=(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(2,N3)-KXYZ(2,N1)) & -(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(1,N3)-KXYZ(1,N1)) C4=-(T4X*KXYZ(1,N1)+T4Y*KXYZ(2,N1)+T4Z*KXYZ(3,N1)) P4C=T4X*KXYZ(1,N4)+T4Y*KXYZ(2,N4)+T4Z*KXYZ(3,N4)+C4 T4X=T4X/P4C T4Y=T4Y/P4C T4Z=T4Z/P4C C4=C4/P4C C C Boucle sur les points de Gauss C DO 50 N=1,(NPG/8) IF (NPG.EQ.32) THEN DO 60 L=1,4 AA(L)=BE 60 CONTINUE AA(N)=AL ELSE DO 70 L=1,4 AA(L)=0.5D1 70 CONTINUE ENDIF C C Indice globale du point de Gauss C IGAU=(K-1)*(NPG/8)+N C C Coordonnees du point de Gauss (barycentre des 4 points du tetraedre) C X(IGAU)=AA(1)*KXYZ(1,N1)+AA(2)*KXYZ(1,N2) & +AA(3)*KXYZ(1,N3)+AA(4)*KXYZ(1,N4) Y(IGAU)=AA(1)*KXYZ(2,N1)+AA(2)*KXYZ(2,N2) & +AA(3)*KXYZ(2,N3)+AA(4)*KXYZ(2,N4) Z(IGAU)=AA(1)*KXYZ(3,N1)+AA(2)*KXYZ(3,N2) & +AA(3)*KXYZ(3,N3)+AA(4)*KXYZ(3,N4) C C Boucle sur les fonctions tests C C C Numerotation globale de la fonction test C IFN=(K-1)*4 C C Calcul des fonctions tests et du gradient au point IGAU C FN((IFN+1),IGAU)=T1X*X(IGAU)+T1Y*Y(IGAU) & +T1Z*Z(IGAU)+C1 GR(1,(IFN+1),IGAU)=T1X GR(2,(IFN+1),IGAU)=T1Y GR(3,(IFN+1),IGAU)=T1Z FN((IFN+2),IGAU)=T2X*X(IGAU)+T2Y*Y(IGAU) & +T2Z*Z(IGAU)+C2 GR(1,(IFN+2),IGAU)=T2X GR(2,(IFN+2),IGAU)=T2Y GR(3,(IFN+2),IGAU)=T2Z FN((IFN+3),IGAU)=T3X*X(IGAU)+T3Y*Y(IGAU) & +T3Z*Z(IGAU)+C3 GR(1,(IFN+3),IGAU)=T3X GR(2,(IFN+3),IGAU)=T3Y GR(3,(IFN+3),IGAU)=T3Z FN((IFN+4),IGAU)=T4X*X(IGAU)+T4Y*Y(IGAU) & +T4Z*Z(IGAU)+C4 GR(1,(IFN+4),IGAU)=T4X GR(2,(IFN+4),IGAU)=T4Y GR(3,(IFN+4),IGAU)=T4Z C C Calcul des fonctions tests pour la pression C FM(((K+1)/2),IGAU)=1.0D0 50 CONTINUE 40 CONTINUE IF(NOM2.EQ.'MCF1')THEN CO=6.D0 ** (1.D0/3.D0) UNSCO=1.D0/CO XXXX=-1.D0*UNSCO DO 51 L=1,NPG FM(1,L)=1.D0-(X(L)+Y(L)+Z(L))*UNSCO FM(2,L)=X(L)*UNSCO FM(3,L)=Y(L)*UNSCO FM(4,L)=Z(L)*UNSCO C GM(1,1,L)=XXXX GM(2,1,L)=XXXX GM(3,1,L)=XXXX C GM(1,2,L)=UNSCO GM(2,2,L)=0.D0 GM(3,2,L)=0.D0 C GM(1,3,L)=0.D0 GM(2,3,L)=UNSCO GM(3,3,L)=0.D0 C GM(1,4,L)=0.D0 GM(2,4,L)=0.D0 GM(3,4,L)=UNSCO C 51 CONTINUE ENDIF C write(6,*)'X ' C write(6,1008)X 1008 FORMAT(8(1X,1PE11.4)) RETURN END