1 : C PB443 SOURCE CHAT 05/01/13 02:10:40 5004 2 : SUBROUTINE PB443(X,Y,Z,PG,FN,GR,FM,GM,ND,NP,MP,NGG,NPG,NOM2) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C************************************************************************ 6 : C 7 : C CALCULE LES FONCTIONS DE FORME D'UN : Iso-Q2 (iso P1/P0 nc) TE10 8 : C 9 : C 10 : C 11 : C ^ zeta 12 : C | 13 : C |10 14 : C | 15 : C | 16 : C | 9 17 : C | /\ /^ eta 18 : C | / \ / 19 : C |/____\ / 20 : C |7 8 / 21 : C | / 22 : C | 5 23 : C | / \ 24 : C | / \ 25 : C | 6_____\4 26 : C | / \ /\ 27 : C | / \ / \ 28 : C |/_____\/____\ ____________________>ksi 29 : C 1 2 3 30 : C 31 : C 32 : C************************************************************************ 33 : CHARACTER*4 NOM2 34 : REAL*8 AL,BE 35 : PARAMETER (AL=.5854101966249684D0) 36 : PARAMETER (BE=.1381966011250105D0) 37 : REAL*8 X(NPG),Y(NPG),Z(NPG) 38 : REAL*8 FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) 39 : 40 : REAL*8 FM(MP,NPG),GM(ND,MP,NPG) 41 : INTEGER CONEK(4,8),L,K,I,N,N1,N2,N3,N4,IGAU,IFN 42 : INTEGER ND,MP,NGG,NP,NPG 43 : REAL*8 KXYZ(3,10),AA(4) 44 : REAL*8 T1X,T1Y,T1Z,C1,P1C 45 : REAL*8 T2X,T2Y,T2Z,C2,P2C 46 : REAL*8 T3X,T3Y,T3Z,C3,P3C 47 : REAL*8 T4X,T4Y,T4Z,C4,P4C 48 : C 49 : C Donnees de la numerotation locale des 4 noeuds de chaque tetraedre 50 : C 51 : DATA CONEK/1,2,6,7, 52 : & 2,8,6,7, 53 : & 2,3,4,8, 54 : & 2,4,6,8, 55 : & 4,5,6,9, 56 : & 8,9,4,6, 57 : & 6,7,8,9, 58 : & 7,8,9,10/ 59 : C 60 : C Donnees des coordonnees des points du tetraedre TE10 61 : C X Y Z 62 : DATA KXYZ/0.0, 0.0, 0.0, 63 : & 0.5, 0.0, 0.0, 64 : & 1.0, 0.0, 0.0, 65 : & 0.5, 0.5, 0.0, 66 : & 0.0, 1.0, 0.0, 67 : & 0.0, 0.5, 0.0, 68 : & 0.0, 0.0, 0.5, 69 : & 0.5, 0.0, 0.5, 70 : & 0.0, 0.5, 0.5, 71 : & 0.0, 0.0, 1.0/ 72 : 73 : C 74 : C On initialise les fonctions tests et les gradients a 0 75 : C 76 : DO 10 L=1,NPG 77 : C Poids de Gauss 78 : PG(L) = 0.1875D0 79 : DO 22 K=1,MP 80 : FM(K,L)=0.0D0 81 : DO 22 I=1,ND 82 : GM(I,K,L)=0.0D0 83 : 22 CONTINUE 84 : 85 : DO 20 K=1,NP 86 : FN(K,L)=0.0D0 87 : DO 20 I=1,ND 88 : GR(I,K,L)=0.0D0 89 : 20 CONTINUE 90 : 91 : 10 CONTINUE 92 : C 93 : C On traite chaque tetraedre elementaire 94 : C 95 : DO 40 K=1,8 96 : C 97 : C On travaille sur chacun des tetraedres composant le macro 98 : C 99 : N1=CONEK(1,K) 100 : N2=CONEK(2,K) 101 : N3=CONEK(3,K) 102 : N4=CONEK(4,K) 103 : C Calcul des equations de plans oppose a chaque point 104 : C 105 : C Plan 1 106 : T1X=(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(3,N4)-KXYZ(3,N2)) 107 : & -(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(2,N4)-KXYZ(2,N2)) 108 : T1Y=(KXYZ(3,N3)-KXYZ(3,N2))*(KXYZ(1,N4)-KXYZ(1,N2)) 109 : & -(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(3,N4)-KXYZ(3,N2)) 110 : T1Z=(KXYZ(1,N3)-KXYZ(1,N2))*(KXYZ(2,N4)-KXYZ(2,N2)) 111 : & -(KXYZ(2,N3)-KXYZ(2,N2))*(KXYZ(1,N4)-KXYZ(1,N2)) 112 : C1=-(T1X*KXYZ(1,N2)+T1Y*KXYZ(2,N2)+T1Z*KXYZ(3,N2)) 113 : P1C=T1X*KXYZ(1,N1)+T1Y*KXYZ(2,N1)+T1Z*KXYZ(3,N1)+C1 114 : T1X=T1X/P1C 115 : T1Y=T1Y/P1C 116 : T1Z=T1Z/P1C 117 : C1=C1/P1C 118 : C Plan 2 119 : T2X=(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(3,N1)-KXYZ(3,N3)) 120 : & -(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(2,N1)-KXYZ(2,N3)) 121 : T2Y=(KXYZ(3,N4)-KXYZ(3,N3))*(KXYZ(1,N1)-KXYZ(1,N3)) 122 : & -(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(3,N1)-KXYZ(3,N3)) 123 : T2Z=(KXYZ(1,N4)-KXYZ(1,N3))*(KXYZ(2,N1)-KXYZ(2,N3)) 124 : & -(KXYZ(2,N4)-KXYZ(2,N3))*(KXYZ(1,N1)-KXYZ(1,N3)) 125 : C2=-(T2X*KXYZ(1,N3)+T2Y*KXYZ(2,N3)+T2Z*KXYZ(3,N3)) 126 : P2C=T2X*KXYZ(1,N2)+T2Y*KXYZ(2,N2)+T2Z*KXYZ(3,N2)+C2 127 : T2X=T2X/P2C 128 : T2Y=T2Y/P2C 129 : T2Z=T2Z/P2C 130 : C2=C2/P2C 131 : C Plan 3 132 : T3X=(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(3,N2)-KXYZ(3,N4)) 133 : & -(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(2,N2)-KXYZ(2,N4)) 134 : T3Y=(KXYZ(3,N1)-KXYZ(3,N4))*(KXYZ(1,N2)-KXYZ(1,N4)) 135 : & -(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(3,N2)-KXYZ(3,N4)) 136 : T3Z=(KXYZ(1,N1)-KXYZ(1,N4))*(KXYZ(2,N2)-KXYZ(2,N4)) 137 : & -(KXYZ(2,N1)-KXYZ(2,N4))*(KXYZ(1,N2)-KXYZ(1,N4)) 138 : C3=-(T3X*KXYZ(1,N4)+T3Y*KXYZ(2,N4)+T3Z*KXYZ(3,N4)) 139 : P3C=T3X*KXYZ(1,N3)+T3Y*KXYZ(2,N3)+T3Z*KXYZ(3,N3)+C3 140 : T3X=T3X/P3C 141 : T3Y=T3Y/P3C 142 : T3Z=T3Z/P3C 143 : C3=C3/P3C 144 : C Plan 4 145 : T4X=(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(3,N3)-KXYZ(3,N1)) 146 : & -(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(2,N3)-KXYZ(2,N1)) 147 : T4Y=(KXYZ(3,N2)-KXYZ(3,N1))*(KXYZ(1,N3)-KXYZ(1,N1)) 148 : & -(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(3,N3)-KXYZ(3,N1)) 149 : T4Z=(KXYZ(1,N2)-KXYZ(1,N1))*(KXYZ(2,N3)-KXYZ(2,N1)) 150 : & -(KXYZ(2,N2)-KXYZ(2,N1))*(KXYZ(1,N3)-KXYZ(1,N1)) 151 : C4=-(T4X*KXYZ(1,N1)+T4Y*KXYZ(2,N1)+T4Z*KXYZ(3,N1)) 152 : P4C=T4X*KXYZ(1,N4)+T4Y*KXYZ(2,N4)+T4Z*KXYZ(3,N4)+C4 153 : T4X=T4X/P4C 154 : T4Y=T4Y/P4C 155 : T4Z=T4Z/P4C 156 : C4=C4/P4C 157 : C 158 : C Boucle sur les points de Gauss 159 : C 160 : DO 50 N=1,(NPG/8) 161 : IF (NPG.EQ.32) THEN 162 : DO 60 L=1,4 163 : AA(L)=BE 164 : 60 CONTINUE 165 : AA(N)=AL 166 : ELSE 167 : DO 70 L=1,4 168 : AA(L)=0.5D1 169 : 70 CONTINUE 170 : ENDIF 171 : C 172 : C Indice globale du point de Gauss 173 : C 174 : IGAU=(K-1)*(NPG/8)+N 175 : C 176 : C Coordonnees du point de Gauss (barycentre des 4 points du tetraedre) 177 : C 178 : X(IGAU)=AA(1)*KXYZ(1,N1)+AA(2)*KXYZ(1,N2) 179 : & +AA(3)*KXYZ(1,N3)+AA(4)*KXYZ(1,N4) 180 : Y(IGAU)=AA(1)*KXYZ(2,N1)+AA(2)*KXYZ(2,N2) 181 : & +AA(3)*KXYZ(2,N3)+AA(4)*KXYZ(2,N4) 182 : Z(IGAU)=AA(1)*KXYZ(3,N1)+AA(2)*KXYZ(3,N2) 183 : & +AA(3)*KXYZ(3,N3)+AA(4)*KXYZ(3,N4) 184 : C 185 : C Boucle sur les fonctions tests 186 : C 187 : C 188 : C Numerotation globale de la fonction test 189 : C 190 : IFN=(K-1)*4 191 : C 192 : C Calcul des fonctions tests et du gradient au point IGAU 193 : C 194 : FN((IFN+1),IGAU)=T1X*X(IGAU)+T1Y*Y(IGAU) 195 : & +T1Z*Z(IGAU)+C1 196 : GR(1,(IFN+1),IGAU)=T1X 197 : GR(2,(IFN+1),IGAU)=T1Y 198 : GR(3,(IFN+1),IGAU)=T1Z 199 : FN((IFN+2),IGAU)=T2X*X(IGAU)+T2Y*Y(IGAU) 200 : & +T2Z*Z(IGAU)+C2 201 : GR(1,(IFN+2),IGAU)=T2X 202 : GR(2,(IFN+2),IGAU)=T2Y 203 : GR(3,(IFN+2),IGAU)=T2Z 204 : FN((IFN+3),IGAU)=T3X*X(IGAU)+T3Y*Y(IGAU) 205 : & +T3Z*Z(IGAU)+C3 206 : GR(1,(IFN+3),IGAU)=T3X 207 : GR(2,(IFN+3),IGAU)=T3Y 208 : GR(3,(IFN+3),IGAU)=T3Z 209 : FN((IFN+4),IGAU)=T4X*X(IGAU)+T4Y*Y(IGAU) 210 : & +T4Z*Z(IGAU)+C4 211 : GR(1,(IFN+4),IGAU)=T4X 212 : GR(2,(IFN+4),IGAU)=T4Y 213 : GR(3,(IFN+4),IGAU)=T4Z 214 : C 215 : C Calcul des fonctions tests pour la pression 216 : C 217 : FM(((K+1)/2),IGAU)=1.0D0 218 : 50 CONTINUE 219 : 40 CONTINUE 220 : 221 : IF(NOM2.EQ.'MCF1')THEN 222 : CO=6.D0 ** (1.D0/3.D0) 223 : UNSCO=1.D0/CO 224 : XXXX=-1.D0*UNSCO 225 : DO 51 L=1,NPG 226 : FM(1,L)=1.D0-(X(L)+Y(L)+Z(L))*UNSCO 227 : FM(2,L)=X(L)*UNSCO 228 : FM(3,L)=Y(L)*UNSCO 229 : FM(4,L)=Z(L)*UNSCO 230 : C 231 : GM(1,1,L)=XXXX 232 : GM(2,1,L)=XXXX 233 : GM(3,1,L)=XXXX 234 : C 235 : GM(1,2,L)=UNSCO 236 : GM(2,2,L)=0.D0 237 : GM(3,2,L)=0.D0 238 : C 239 : GM(1,3,L)=0.D0 240 : GM(2,3,L)=UNSCO 241 : GM(3,3,L)=0.D0 242 : C 243 : GM(1,4,L)=0.D0 244 : GM(2,4,L)=0.D0 245 : GM(3,4,L)=UNSCO 246 : C 247 : 51 CONTINUE 248 : ENDIF 249 : 250 : C write(6,*)'X ' 251 : C write(6,1008)X 252 : 1008 FORMAT(8(1X,1PE11.4)) 253 : 254 : RETURN 255 : END 256 : 257 :
© Cast3M 2003 - All rights reserved.
Disclaimer