C PB663 SOURCE CHAT 05/01/13 02:10:51 5004 SUBROUTINE PB663(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) PR18 C C P4 equation du plan zeta=0 C P5 equation du plan zeta=1 C ^ zeta etc C | 14 ^ eta |X3 3 C | / \ / |Y3 o C | / \ / |\ C | 15_____\13 / | \ C | / \ /\ / P3 | \ P2 C | / \ / \ / | \ C |/_____\/____\ / | \ C 10 11 12 / | \ C | / 1 o______o 2 C | 9 / |X1 |X2 C | / \ / |Y1 P1 |Y2 C | / \ / C | 18_____17 C | / \ /\ XI abcisse intersection P2 et ksi C | / \ / \ YI ordonnée intersection P2 et eta C |/_____\/____\ C 7 16 8 C | / fct de formes N1=P2*P5/F1C C | 5 tel que N1=1 au pt 1 C | / \ C | / \ C | 6_____\4 C | / \ /\ C | / \ / \ C |/_____\/____\ ____________________>ksi C 1 2 3 C C C************************************************************************ CHARACTER*4 NOM2 REAL*8 X(NPG),Y(NPG),Z(NPG) DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) DIMENSION FM(MP,NPG),GM(ND,MP,NPG) DIMENSION U(5),H(5),XA(3),XB(3),XC(3),XD(3) SAVE XA,XB,XC,XD DATA XA/3*0.25D0/,XB/0.75D0,0.25D0,0.75D0/ DATA XC/2*0.75D0,0.25D0/,XD/0.25D0,2*0.75D0/ C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C R2=SQRT(2.D0) R22=R2*0.5D0 DO 1 K=1,8 C IF(K.EQ.1)THEN X1=0.D0 Y1=0.D0 X2=R22 Y2=0.D0 X3=0.D0 Y3=R22 Z1=0.D0 Z2=0.5D0 XI=R22 YI=R22 ELSEIF(K.EQ.2)THEN X1=R22 Y1=0.D0 X2=R2 Y2=0.D0 X3=R22 Y3=R22 Z1=0.D0 Z2=0.5D0 XI=R2 YI=R2 ELSEIF(K.EQ.3)THEN X1=0.D0 Y1=R22 X2=R22 Y2=R22 X3=0.D0 Y3=R2 Z1=0.D0 Z2=0.5D0 XI=R2 YI=R2 ELSEIF(K.EQ.4)THEN X1=R22 Y1=R22 X2=0.D0 Y2=R22 X3=R22 Y3=0.D0 Z1=0.D0 Z2=0.5D0 XI=R22 YI=R22 ELSEIF(K.EQ.5)THEN X1=0.D0 Y1=0.D0 X2=R22 Y2=0.D0 X3=0.D0 Y3=R22 Z1=0.5D0 Z2=1.D0 XI=R22 YI=R22 ELSEIF(K.EQ.6)THEN X1=R22 Y1=0.D0 X2=R2 Y2=0.D0 X3=R22 Y3=R22 Z1=0.5D0 Z2=1.D0 XI=R2 YI=R2 ELSEIF(K.EQ.7)THEN X1=0.D0 Y1=R22 X2=R22 Y2=R22 X3=0.D0 Y3=R2 Z1=0.5D0 Z2=1.D0 XI=R2 YI=R2 ELSEIF(K.EQ.8)THEN X1=R22 Y1=R22 X2=0.D0 Y2=R22 X3=R22 Y3=0.D0 Z1=0.5D0 Z2=1.D0 XI=R22 YI=R22 ENDIF K0=(K-1)*NGG CALL KHPRM(X(K0+1),Y(K0+1),Z(K0+1),PG(K0+1),NGG, &X1,Y1,X2,Y3,Z1,Z2) C write(6,*)' KHPRM NGG=',ngg C write(6,1002)(x(ii),ii=1,ngg) C write(6,1002)(y(ii),ii=1,ngg) C write(6,1002)(z(ii),ii=1,ngg) P1C=Y3-Y1 P2C=X1/XI+Y1/YI-1.D0 P3C=X2-X1 P4C=Z2-Z1 P5C=Z1-Z2 C write(6,*)'PiC ' C write(6,1002)P1C,P2C,P3C,P4C,P5C DO 1 L=1,NGG P1=Y(L+K0)-Y1 P2=X(L+K0)/XI+Y(L+K0)/YI-1.D0 P3=X(L+K0)-X1 P4=Z(L+K0)-Z1 P5=Z(L+K0)-Z2 F1C=P2C*P5C F1=P2*P5/F1C F2C=P3C*P5C F2=P3*P5/F2C F3C=P1C*P5C F3=P1*P5/F3C F4C=P2C*P4C F4=P2*P4/F4C F5C=P3C*P4C F5=P3*P4/F5C F6C=P1C*P4C F6=P1*P4/F6C C write(6,*)'FiC' C write(6,1002)F1C,F2C,F3C,F4C,F5C,F6C GX1=P5/XI/F1C GY1=P5/YI/F1C GZ1=P2/F1C GX2=P5/F2C GY2=0.D0 GZ2=P3/F2C GX3=0.D0 GY3=P5/F3C GZ3=P1/F3C GX4=P4/XI/F4C GY4=P4/YI/F4C GZ4=P2/F4C GX5=P4/F5C GY5=0.D0 GZ5=P3/F5C GX6=0.D0 GY6=P4/F6C GZ6=P1/F6C C write(6,1002)F1,f2,f3,f4,f5,f6 C C write(6,1002)gx1,gx2,gx3,gx4,gx5,gx6 C write(6,1002)gy1,gy2,gy3,gy4,gy5,gy6 C write(6,1002)gz1,gz2,gz3,gz4,gz5,gz6 LL=K0+L IF(K.EQ.1)THEN FN( 1,LL)=F1 FN( 2,LL)=F2 FN( 6,LL)=F3 FN( 7,LL)=F4 FN(16,LL)=F5 FN(18,LL)=F6 GR(1, 1,LL)=GX1 GR(1, 2,LL)=GX2 GR(1, 6,LL)=GX3 GR(1, 7,LL)=GX4 GR(1,16,LL)=GX5 GR(1,18,LL)=GX6 GR(2, 1,LL)=GY1 GR(2, 2,LL)=GY2 GR(2, 6,LL)=GY3 GR(2, 7,LL)=GY4 GR(2,16,LL)=GY5 GR(2,18,LL)=GY6 GR(3, 1,LL)=GZ1 GR(3, 2,LL)=GZ2 GR(3, 6,LL)=GZ3 GR(3, 7,LL)=GZ4 GR(3,16,LL)=GZ5 GR(3,18,LL)=GZ6 ELSEIF(K.EQ.2)THEN FN( 2,LL)=F1 FN( 3,LL)=F2 FN( 4,LL)=F3 FN(16,LL)=F4 FN( 8,LL)=F5 FN(17,LL)=F6 GR(1, 2,LL)=GX1 GR(1, 3,LL)=GX2 GR(1, 4,LL)=GX3 GR(1,16,LL)=GX4 GR(1, 8,LL)=GX5 GR(1,17,LL)=GX6 GR(2, 2,LL)=GY1 GR(2, 3,LL)=GY2 GR(2, 4,LL)=GY3 GR(2,16,LL)=GY4 GR(2, 8,LL)=GY5 GR(2,17,LL)=GY6 GR(3, 2,LL)=GZ1 GR(3, 3,LL)=GZ2 GR(3, 4,LL)=GZ3 GR(3,16,LL)=GZ4 GR(3, 8,LL)=GZ5 GR(3,17,LL)=GZ6 ELSEIF(K.EQ.3)THEN FN( 6,LL)=F1 FN( 4,LL)=F2 FN( 5,LL)=F3 FN(18,LL)=F4 FN(17,LL)=F5 FN( 9,LL)=F6 GR(1, 6,LL)=GX1 GR(1, 4,LL)=GX2 GR(1, 5,LL)=GX3 GR(1,18,LL)=GX4 GR(1,17,LL)=GX5 GR(1, 9,LL)=GX6 GR(2, 6,LL)=GY1 GR(2, 4,LL)=GY2 GR(2, 5,LL)=GY3 GR(2,18,LL)=GY4 GR(2,17,LL)=GY5 GR(2, 9,LL)=GY6 GR(3, 6,LL)=GZ1 GR(3, 4,LL)=GZ2 GR(3, 5,LL)=GZ3 GR(3,18,LL)=GZ4 GR(3,17,LL)=GZ5 GR(3, 9,LL)=GZ6 ELSEIF(K.EQ.4)THEN FN( 4,LL)=F1 FN( 6,LL)=F2 FN( 2,LL)=F3 FN(17,LL)=F4 FN(18,LL)=F5 FN(16,LL)=F6 GR(1, 4,LL)=GX1 GR(1, 6,LL)=GX2 GR(1, 2,LL)=GX3 GR(1,17,LL)=GX4 GR(1,18,LL)=GX5 GR(1,16,LL)=GX6 GR(2, 4,LL)=GY1 GR(2, 6,LL)=GY2 GR(2, 2,LL)=GY3 GR(2,17,LL)=GY4 GR(2,18,LL)=GY5 GR(2,16,LL)=GY6 GR(3, 4,LL)=GZ1 GR(3, 6,LL)=GZ2 GR(3, 2,LL)=GZ3 GR(3,17,LL)=GZ4 GR(3,18,LL)=GZ5 GR(3,16,LL)=GZ6 ELSEIF(K.EQ.5)THEN FN( 7,LL)=F1 FN(16,LL)=F2 FN(18,LL)=F3 FN(10,LL)=F4 FN(11,LL)=F5 FN(15,LL)=F6 GR(1, 7,LL)=GX1 GR(1,16,LL)=GX2 GR(1,18,LL)=GX3 GR(1,10,LL)=GX4 GR(1,11,LL)=GX5 GR(1,15,LL)=GX6 GR(2, 7,LL)=GY1 GR(2,16,LL)=GY2 GR(2,18,LL)=GY3 GR(2,10,LL)=GY4 GR(2,11,LL)=GY5 GR(2,15,LL)=GY6 GR(3, 7,LL)=GZ1 GR(3,16,LL)=GZ2 GR(3,18,LL)=GZ3 GR(3,10,LL)=GZ4 GR(3,11,LL)=GZ5 GR(3,15,LL)=GZ6 ELSEIF(K.EQ.6)THEN FN(16,LL)=F1 FN( 8,LL)=F2 FN(17,LL)=F3 FN(11,LL)=F4 FN(12,LL)=F5 FN(13,LL)=F6 GR(1,16,LL)=GX1 GR(1, 8,LL)=GX2 GR(1,17,LL)=GX3 GR(1,11,LL)=GX4 GR(1,12,LL)=GX5 GR(1,13,LL)=GX6 GR(2,16,LL)=GY1 GR(2, 8,LL)=GY2 GR(2,17,LL)=GY3 GR(2,11,LL)=GY4 GR(2,12,LL)=GY5 GR(2,13,LL)=GY6 GR(3,16,LL)=GZ1 GR(3, 8,LL)=GZ2 GR(3,17,LL)=GZ3 GR(3,11,LL)=GZ4 GR(3,12,LL)=GZ5 GR(3,13,LL)=GZ6 ELSEIF(K.EQ.7)THEN FN(18,LL)=F1 FN(17,LL)=F2 FN( 9,LL)=F3 FN(15,LL)=F4 FN(13,LL)=F5 FN(14,LL)=F6 GR(1,18,LL)=GX1 GR(1,17,LL)=GX2 GR(1, 9,LL)=GX3 GR(1,15,LL)=GX4 GR(1,13,LL)=GX5 GR(1,14,LL)=GX6 GR(2,18,LL)=GY1 GR(2,17,LL)=GY2 GR(2, 9,LL)=GY3 GR(2,15,LL)=GY4 GR(2,13,LL)=GY5 GR(2,14,LL)=GY6 GR(3,18,LL)=GZ1 GR(3,17,LL)=GZ2 GR(3, 9,LL)=GZ3 GR(3,15,LL)=GZ4 GR(3,13,LL)=GZ5 GR(3,14,LL)=GZ6 ELSEIF(K.EQ.8)THEN FN(17,LL)=F1 FN(18,LL)=F2 FN(16,LL)=F3 FN(13,LL)=F4 FN(15,LL)=F5 FN(11,LL)=F6 GR(1,17,LL)=GX1 GR(1,18,LL)=GX2 GR(1,16,LL)=GX3 GR(1,13,LL)=GX4 GR(1,15,LL)=GX5 GR(1,11,LL)=GX6 GR(2,17,LL)=GY1 GR(2,18,LL)=GY2 GR(2,16,LL)=GY3 GR(2,13,LL)=GY4 GR(2,15,LL)=GY5 GR(2,11,LL)=GY6 GR(3,17,LL)=GZ1 GR(3,18,LL)=GZ2 GR(3,16,LL)=GZ3 GR(3,13,LL)=GZ4 GR(3,15,LL)=GZ5 GR(3,11,LL)=GZ6 ENDIF 1 CONTINUE IF(NOM2.EQ.'MCP0')THEN DO 2 L=1,NPG FM(1,L)=1.D0 GM(1,1,L)=0.D0 GM(2,1,L)=0.D0 GM(3,1,L)=0.D0 2 CONTINUE ELSEIF(NOM2.EQ.'MCP1')THEN DO 3 LL=1,(2*NGG) FM(1,LL)=0.25D0 FM(2,LL+2*NGG)=0.25D0 FM(3,LL+4*NGG)=0.25D0 FM(4,LL+6*NGG)=0.25D0 3 CONTINUE CALL INITD(GM,(12*NPG),0.D0) ELSEIF(NOM2.EQ.'MCF1')THEN DO 4 L=1,NPG FM(1,L)=(X(L)+Y(L)-R2)*(Z(L)-1.D0)/R2 FM(2,L)=-X(L)*(Z(L)-1.D0)/R2 FM(3,L)=-Y(L)*(Z(L)-1.D0)/R2 FM(4,L)=-(X(L)+Y(L)-R2)*Z(L)/R2 FM(5,L)=X(L)*Z(L)/R2 FM(6,L)=Y(L)*Z(L)/R2 GM(1,1,L)=(Z(L)-1.D0)/R2 GM(2,1,L)=(Z(L)-1.D0)/R2 GM(3,1,L)=(X(L)+Y(L)-R2)/R2 GM(1,2,L)=-(Z(L)-1.D0)/R2 GM(2,2,L)=0.D0 GM(3,2,L)=-X(L)/R2 GM(1,3,L)=0.D0 GM(2,3,L)=-(Z(L)-1.D0)/R2 GM(3,3,L)=-Y(L)/R2 GM(1,4,L)=-Z(L)/R2 GM(2,4,L)=-Z(L)/R2 GM(3,4,L)=-(X(L)+Y(L)-R2)/R2 GM(1,5,L)=Z(L)/R2 GM(2,5,L)=0.D0 GM(3,5,L)=X(L)/R2 GM(1,6,L)=0.D0 GM(2,6,L)=Z(L)/R2 GM(3,6,L)=Y(L)/R2 C 4 CONTINUE ENDIF C write(6,*)' FM' C write(6,1002)((FM(II,JJ),II=1,MP),JJ=1,NPG) C write(6,*)' VERIF ,PG=' C write(6,1002)(pg(ii),ii=1,npg) C do 75 ll=1,npg C write(6,*)' VERIF ,fn,gr ll=',ll C write(6,1002)(fn(ii,ll),ii=1,np) C write(6,1002)(gr(1,ii,ll),ii=1,np) C write(6,1002)(gr(2,ii,ll),ii=1,np) C write(6,1002)(gr(3,ii,ll),ii=1,np) C75 continue C write(6,*)' VERIF ,NPG,NP,NGG=',NPG,NP,NGG C UPG=0.D0 C do 72 L=1,NPG C UPG=UPG+PG(L) C UF=0.D0 C UG1=0.D0 C UG2=0.D0 C UG3=0.D0 C DO 71 I=1,NP C UF=UF+FN(I,L) C UG1=UG1+GR(1,I,L) C UG2=UG2+GR(2,I,L) C UG3=UG3+GR(3,I,L) C71 CONTINUE C WRITE(6,*)' VERIF L=',L,UF,UG1,UG2,UG3 C72 CONTINUE C WRITE(6,*)' VERIF PG=',UPG C WRITE(6,101) C IF(MP.EQ.1)THEN C DO 4 L=1,NPG C FM(1,L)=1.D0 C GM(1,1,L)=0.D0 C GM(2,1,L)=0.D0 C GM(3,1,L)=0.D0 C4 CONTINUE C ELSEIF(MP.EQ.4)THEN C XX(1)=X(L) C XX(2)=Y(L) C XX(3)=Z(L) C FM(1,L)=EQPL3P(XX,XB,XC,XD)/ C & EQPL3P(XA,XB,XC,XD) C FM(2,L)=EQPL3P(XX,XA,XC,XD)/ C & EQPL3P(XB,XA,XC,XD) C FM(3,L)=EQPL3P(XX,XA,XB,XD)/ C & EQPL3P(XC,XA,XB,XD) C C FM(4,L)=EQPL3P(XX,XA,XB,XC)/ C & EQPL3P(XD,XA,XB,XC) C CALL INITD(GM,(12*NPG),0.D0) C3 CONTINUE C ENDIF C WRITE(6,101) RETURN 1002 FORMAT(10(1X,1PD11.4)) 1001 FORMAT(20(1X,I5)) 101 FORMAT(1X,'... SUBPB663 ... FM,GM,FN,GR ',9(10H..........)/) C END