1 : C PB702 SOURCE MAGN 10/05/31 21:15:16 6679 2 : SUBROUTINE PB702(XREF,X,Y,PG,FN,GR,FM,GM,ND,NP,MP,NPG,NOM2,ITYPI) 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 : TRI7 8 : C 9 : C Cf Crouzeix Raviard. Conforming and nonconforming finite element 10 : C methods for solving the stationary Stokes equations I 11 : C R.A.I.R.O. (7e année, décembre 1973,R-3,p.33 à 76) 12 : C 13 : C 14 : C ^ eta 15 : C | 16 : C |n3 17 : C r2 5 1 -> A1 18 : C |\ 2 -> A12 - 19 : C | \ r2=sqrt(2) 3 -> A2 20 : C | \ 4 -> A23 21 : C | \ 5 -> A3 22 : C | \ 6 -> A13 23 : C | \ 7 -> A123 24 : C r2/2 6------4 25 : C | |\ D1=1-x/r2-y/r2 26 : C | | \ D2=x/r2 27 : C | 7 | \ D3=y/r2 28 : C | | \ 29 : C | | \ 30 : C 1------2-----3--->ksi 31 : C 0 r2/2 r2 32 : C 33 : C Pi = Di(aDi-c)+FD1xD2xD3 1 <= i <= 3 34 : C Pij = bxDi(aDi-c)-GxD1xD2xD3 1 <= i < j <= 3 35 : C P123= QxD1xD2xD3 36 : C Pour l'élément de référence on a a=2 b=4 c=1 F=3 G=12 37 : C Si on veut perturber un peu la position des points il faut néanmoins 38 : C vérifier les propriétés essentielles des fonctions de forme qui outre 39 : C Pi=1 en X=i et Pi=0 ailleur on doit avoir Somme des Pi = 1 40 : C (intégration exacte d'une fonction constante). Ceci conduit à des 41 : C relations entre a,b,c,F,G et Q 42 : C b=2a c=a-1 Q=3(G-F) (Somme Pi = 1) a ne 2 points Aij non au milieux 43 : C des arètes 44 : C 45 : C ITYPI=2 pts d'integration de type Gauss Lobatto i.e. sommets element 46 : C 47 : C REMARQUE : 48 : C En coordonnées cylindriques,pour une distribution régulière des noeuds 49 : C conforme à l'article de Crouzeix Raviard, la matrice masse diagonale 50 : C s'annule sur l'axe, comme pour le QUA9 (voir pb902.eso). 51 : C /1 52 : C | 2 pi x N1 dx = 0. 53 : C /0 54 : C Pour remédier à cela on se sert de la bulle en la décalant un peu par 55 : C rapport à l'élément de référence en jouant sur les constantes F et G. 56 : C On a pas à toucher "a" position des points 2,4 et 6. 57 : C Le choix que l'on a fait conduit à une matrice masse diagonale >0. En 58 : C conséquence ce choix doit être accompagné d'une numérotation dans le 59 : C sens trigo de l'élément courant. 60 : C 61 : C*********************************************************************** 62 : 63 : REAL*8 XREF(ND,NP),X(NPG),Y(NPG) 64 : DIMENSION FN(NP,NPG),GR(ND,NP,NPG),PG(NPG) 65 : DIMENSION FM(MP,NPG),GM(ND,MP,NPG) 66 : REAL*8 A,B,C,D,U(5),H(5) 67 : CHARACTER*4 NOM2 68 : -INC CCOPTIO 69 : C*** 70 : 71 : IF(IFOMOD.EQ.0)THEN 72 : F=3.0001 73 : G=11.9999 74 : F=3. 75 : G=12. 76 : ELSE 77 : F=3. 78 : G=12. 79 : ENDIF 80 : Q=3.*(G-F) 81 : a=2. 82 : b=2.*a 83 : c=a-1. 84 : 85 : 86 : R2=SQRT(2.D0) 87 : IPOSP=2 88 : 89 : 90 : XREF(1,1)=0.D0 91 : XREF(2,1)=0.D0 92 : 93 : XREF(1,2)=R2/2.D0 94 : XREF(2,2)=0.D0 95 : 96 : XREF(1,3)=R2 97 : XREF(2,3)=0.D0 98 : 99 : XREF(1,4)=R2/2.D0 100 : XREF(2,4)=R2/2.D0 101 : 102 : XREF(1,5)=0.D0 103 : XREF(2,5)=R2 104 : 105 : XREF(1,6)=0. 106 : XREF(2,6)=R2/2.D0 107 : 108 : XREF(1,7)=R2/3.D0 109 : XREF(2,7)=R2/3.D0 110 : 111 : C Verification des coordonnées 112 : c IF(.TRUE.)THEN 113 : IF(.FALSE.)THEN 114 : DO 11 L=1,NP 115 : X(L)=XREF(1,L) 116 : Y(L)=XREF(2,L) 117 : 11 CONTINUE 118 : 119 : DO 12 L=1,NP 120 : 121 : D1=1.D0-X(L)/R2-Y(L)/R2 122 : D2=X(L)/R2 123 : D3=Y(L)/R2 124 : 125 : FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3 126 : FN(2,L)=b*D1*D2-G*D1*D2*D3 127 : FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3 128 : FN(4,L)=b*D2*D3-G*D1*D2*D3 129 : FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3 130 : FN(6,L)=b*D1*D3-G*D1*D2*D3 131 : FN(7,L)=Q*D1*D2*D3 132 : write(6,1033)L,FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L) 133 : & ,FN(7,L) 134 : 12 CONTINUE 135 : 1033 FORMAT(1X,I4,' FN',10(1X,1PD11.4)) 136 : ENDIF 137 : C Fin Vérification 138 : 139 : IF(ITYPI.EQ.2)THEN 140 : X(1)=0.D0 141 : Y(1)=0.D0 142 : X(2)=R2/2.D0 143 : Y(2)=0.D0 144 : X(3)=R2 145 : Y(3)=0.D0 146 : X(4)=R2/2.D0 147 : Y(4)=R2/2.D0 148 : X(5)=0.D0 149 : Y(5)=R2 150 : X(6)=0.D0 151 : Y(6)=R2/2.D0 152 : X(7)=(X(1)+X(3)+X(5))/3.D0 153 : Y(7)=(Y(1)+Y(3)+Y(5))/3.D0 154 : DO 2 L=1,7 155 : PG(L)=1.D0/7.D0 156 : 2 CONTINUE 157 : ELSE 157 : CALL CALUHH(X,Y,PG,NPG) 158 : ENDIF 159 : 160 : DO 1 L=1,NPG 161 : C 162 : IF(NOM2.EQ.'PRP0')THEN 163 : FM(1,L)=1.D0 164 : GM(1,1,L)=0.D0 165 : GM(2,1,L)=0.D0 166 : ELSEIF(NOM2.EQ.'PRP1')THEN 167 : C? FM(1,L)=R2*(X(L)+Y(L)-R2/2.D0) 168 : C? FM(2,L)=-R2*(X(L)-R2/2.D0) 169 : C? FM(3,L)=-R2*(Y(L)-R2/2.D0) 170 : A1=6.D0/5.D0/R2 171 : FM(1,L)=5.D0/3.D0*(1.D0-A1*X(L)-A1*Y(L)) 172 : FM(2,L)=-1.D0/3.D0*(1.D0-6.D0/R2*X(L)) 173 : FM(3,L)=-1.D0/3.D0*(1.D0-6.D0/R2*Y(L)) 174 : C 175 : GM(1,1,L)=0.D0 176 : GM(2,1,L)=0.D0 177 : GM(1,2,L)=0.D0 178 : GM(2,2,L)=0.D0 179 : GM(1,3,L)=0.D0 180 : GM(2,3,L)=0.D0 181 : ELSEIF(NOM2.EQ.'PFP1')THEN 182 : FM(1,L)=-R2/2.D0*(X(L)+Y(L)-R2) 183 : FM(2,L)=R2/2.D0*X(L) 184 : FM(3,L)=R2/2.D0*Y(L) 185 : 186 : GM(1,1,L)=-R2/2.D0 187 : GM(2,1,L)=-R2/2.D0 188 : GM(1,2,L)=R2/2.D0 189 : GM(2,2,L)=0.D0 190 : GM(1,3,L)=0.D0 191 : GM(2,3,L)=R2/2.D0 192 : ENDIF 193 : 194 : C 195 : c IF(.FALSE.)THEN 196 : IF(.TRUE.)THEN 197 : C 198 : 199 : D1=1.D0-X(L)/R2-Y(L)/R2 200 : D2=X(L)/R2 201 : D3=Y(L)/R2 202 : 203 : FN(1,L)=D1*(a*D1-c)+F*D1*D2*D3 204 : FN(2,L)=b*D1*D2-G*D1*D2*D3 205 : FN(3,L)=D2*(a*D2-c)+F*D1*D2*D3 206 : FN(4,L)=b*D2*D3-G*D1*D2*D3 207 : FN(5,L)=D3*(a*D3-c)+F*D1*D2*D3 208 : FN(6,L)=b*D1*D3-G*D1*D2*D3 209 : FN(7,L)=Q*D1*D2*D3 210 : c write(6,*)FN(1,L),FN(2,L),FN(3,L),FN(4,L),FN(5,L),FN(6,L),FN(7,L) 211 : 212 : D1x=-1./R2 213 : D1y=-1./R2 214 : 215 : D2x=1./R2 216 : D2y=0. 217 : 218 : D3x=0. 219 : D3y=1/R2 220 : 221 : GR(1,1,L)=D1x*(a*D1-c)+D1*a*D1x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 222 : GR(2,1,L)=D1y*(a*D1-c)+D1*a*D1y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 223 : 224 : GR(1,2,L)=b*(D1x*D2+D1*D2x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 225 : GR(2,2,L)=b*(D1y*D2+D1*D2y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 226 : 227 : GR(1,3,L)=D2x*(a*D2-c)+D2*a*D2x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 228 : GR(2,3,L)=D2y*(a*D2-c)+D2*a*D2y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 229 : 230 : GR(1,4,L)=b*(D2x*D3+D2*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 231 : GR(2,4,L)=b*(D2y*D3+D2*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 232 : 233 : GR(1,5,L)=D3x*(a*D3-c)+D3*a*D3x+F*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 234 : GR(2,5,L)=D3y*(a*D3-c)+D3*a*D3y+F*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 235 : 236 : GR(1,6,L)=b*(D1x*D3+D1*D3x)-G*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 237 : GR(2,6,L)=b*(D1y*D3+D1*D3y)-G*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 238 : 239 : GR(1,7,L)=Q*(D1x*D2*D3+D1*D2x*D3+D1*D2*D3x) 240 : GR(2,7,L)=Q*(D1y*D2*D3+D1*D2y*D3+D1*D2*D3y) 241 : 242 : ELSE 243 : 244 : c write(6,*)'La base' 245 : 246 : C 247 : C 248 : FN(1,L)=(1.D0-X(L)/R2-Y(L)/R2)*(1.D0-2.D0*X(L)/R2-2.D0*Y(L)/R2+ 249 : & 3.D0/2.D0*X(L)*Y(L)) 250 : FN(2,L)=(1.D0-X(L)/R2-Y(L)/R2)*(4.D0/R2-6.D0*Y(L))*X(L) 251 : FN(3,L)=X(L)/R2*(2.D0*X(L)/R2-1.D0+3.D0*(1.D0-X(L)/R2-Y(L)/R2) 252 : & *Y(L)/R2) 253 : FN(4,L)=2.D0*X(L)*Y(L)*(1.D0-3.D0*(1.D0-X(L)/R2-Y(L)/R2)) 254 : FN(5,L)=Y(L)*Y(L)-Y(L)*R2/2.D0+3.D0/2.D0*X(L)*Y(L)- 255 : & 3.D0*R2/4.D0*X(L)*X(L)*Y(L)-3.D0*R2/4.D0*X(L)*Y(L)*Y(L) 256 : FN(6,L)=2.D0*R2*(-R2/2.D0*Y(L)*Y(L)+3.D0/2.D0*(Y(L)*Y(L)*X(L)+ 257 : & X(L)*X(L)*Y(L))-2.D0*R2*X(L)*Y(L)+Y(L)) 258 : FN(7,L)=27.D0/2.D0*(X(L)*Y(L)-R2/2.D0*X(L)*X(L)*Y(L)- 259 : &R2/2.D0*X(L)*Y(L)*Y(L)) 260 : C 261 : C 262 : C 263 : GR(1,1,L)=-3.D0*R2/4.D0*Y(L)*Y(L)+2.D0*X(L)+7./2.D0*Y(L)- 264 : * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0 265 : GR(2,1,L)=-3.D0*R2/4.D0*X(L)*X(L)+2.D0*Y(L)+7./2.D0*X(L)- 266 : * 3.D0*R2/2.D0*X(L)*Y(L)-3.D0*R2/2.D0 267 : 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)+ 268 : * 2.D0*R2 269 : GR(2,2,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L) 270 : GR(1,3,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L)+ 271 : * 2.D0*X(L)+3.D0/2.D0*Y(L)-R2/2.D0 272 : GR(2,3,L)=R2/2.D0*(-3.D0/2.D0*X(L)*X(L)-3.D0*X(L)*Y(L)+ 273 : * 3.D0*R2/2.D0*X(L)) 274 : GR(1,4,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-4.D0*Y(L) 275 : GR(2,4,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-4.D0*X(L) 276 : GR(1,5,L)=-3.D0*R2/4.D0*Y(L)*Y(L)-3.D0*R2/2.D0*X(L)*Y(L) 277 : & +3.D0/2.D0*Y(L) 278 : GR(2,5,L)=-3.D0*R2/4.D0*X(L)*X(L)-3.D0*R2/2.D0*X(L)*Y(L) 279 : & +2.D0*Y(L)+ 280 : * 3.D0/2.D0*X(L)-R2/2.D0 281 : GR(1,6,L)=3.D0*R2*Y(L)*Y(L)+6.D0*R2*X(L)*Y(L)-8.D0*Y(L) 282 : GR(2,6,L)=3.D0*R2*X(L)*X(L)+6.D0*R2*X(L)*Y(L)-8.D0*X(L) 283 : & -4.D0*Y(L)+2.D0*R2 284 : GR(1,7,L)=27.D0/2.D0*(-R2/2*Y(L)*Y(L)-R2*X(L)*Y(L)+Y(L)) 285 : GR(2,7,L)=27.D0/2.D0*(-R2/2*X(L)*X(L)-R2*X(L)*Y(L)+X(L)) 286 : C 287 : ENDIF 288 : 1 CONTINUE 289 : 290 : IF(.FALSE.)THEN 291 : c IF(.TRUE.)THEN 292 : DO 121 L=1,NPG 293 : F1=0. 294 : DO 131 I=1,NP 295 : F1=F1+FN(I,L) 296 : 131 CONTINUE 297 : Write(6,*)' L=',L,F1 298 : 121 CONTINUE 299 : 300 : DO 122 L=1,NPG 301 : G1=0. 302 : G2=0. 303 : DO 141 I=1,NP 304 : G1=G1+GR(1,I,L) 305 : G2=G2+GR(2,I,L) 306 : 141 CONTINUE 307 : Write(6,*)' L=',L,G1,G2 308 : 122 CONTINUE 309 : ENDIF 310 : 311 : c stop 312 : 313 : 314 : C 315 : C 316 : C WRITE(6,101) 317 : C WRITE(6,1002)FM 318 : C WRITE(6,1002)GM 319 : C WRITE(6,1002)FN 320 : C WRITE(6,1002)GR 321 : C WRITE(6,101) 322 : 323 : RETURN 324 : 1002 FORMAT(10(1X,1PD11.4)) 325 : 1001 FORMAT(20(1X,I5)) 326 : 101 FORMAT(1X,'... SUBPB702 ... FM,GM,FN,GR ',9(10H..........)/) 327 : C 328 : END 329 : 330 : 331 : 332 : 333 : 334 : 335 : 336 :
© Cast3M 2003 - All rights reserved.
Disclaimer