1 : C TESTPG SOURCE GOUNAND 05/12/21 21:37:20 5281 2 : SUBROUTINE TESTPG() 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : TESTPG 7 : C PROJET : Noyau linéaire NLIN 8 : C DESCRIPTION : On vérifie les méthodes d'intégration numériques. On 9 : C vérifie l'exactitude de l'intégration des monomes 10 : C d'ordre inférieur ou égal à celui de la méthode. 11 : C 12 : C LANGAGE : ESOPE 13 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 14 : C mél : gounand@semt2.smts.cea.fr 15 : C*********************************************************************** 16 : C APPELES : 17 : C APPELES (E/S) : 18 : C APPELES (BLAS) : 19 : C APPELES (CALCUL) : 20 : C APPELE PAR : PILOT 21 : C*********************************************************************** 22 : C SYNTAXE GIBIANE : 23 : C ENTREES : 24 : C ENTREES/SORTIES : 25 : C SORTIES : 26 : C*********************************************************************** 27 : C VERSION : v1, 22/07/99, version initiale 28 : C HISTORIQUE : v1, 22/07/99, création 29 : C HISTORIQUE : 30 : C HISTORIQUE : 31 : C*********************************************************************** 32 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 33 : C en cas de modification de ce sous-programme afin de faciliter 34 : C la maintenance ! 35 : C*********************************************************************** 36 : -INC CCOPTIO 37 : CBEGININCLUDE SPOGAU 38 : SEGMENT POGAU 39 : CHARACTER*(LNNPG) NOMPG 40 : CHARACTER*(LNTPG) TYPMPG 41 : CHARACTER*(LNFPG) FORLPG 42 : INTEGER NORDPG 43 : REAL*8 XCOPG(NDLPG,NBPG) 44 : REAL*8 XPOPG(NBPG) 45 : ENDSEGMENT 46 : SEGMENT POGAUS 47 : POINTEUR LISPG(0).POGAU 48 : ENDSEGMENT 49 : CENDINCLUDE SPOGAU 50 : POINTEUR MYPGS.POGAUS 51 : POINTEUR PGCOUR.POGAU 52 : * 53 : * Exposant des monomes 54 : * 55 : SEGMENT XPMONO 56 : INTEGER EXPOSA(NDESP) 57 : ENDSEGMENT 58 : INTEGER NDESP 59 : POINTEUR MYMONO.XPMONO 60 : * 61 : INTEGER IMPR,IRET 62 : * 63 : INTEGER NBMETH,I,J,K 64 : INTEGER NORDM,NESREF,NPOGAU 65 : INTEGER IMETH,IORD,IEXP,IEXP1,IEXP2,IPG,ICOOR 66 : REAL*8 XFACT,XERREL 67 : REAL*8 XPRMON,SOLEX,SOLAPP,SERR 68 : * 69 : * Executable statements 70 : * 71 : * 72 : * Initialisation du segment contenant les informations sur les 73 : * méthodes d'intégration (type Gauss). 74 : * 75 : IMPR=0 75 : CALL INPGS(MYPGS,IMPR,IRET) 76 : IF (IRET.NE.0) GOTO 9999 77 : * 78 : * Puis on teste chacune de méthodes 79 : * 80 : * IMPR=0 81 : SEGACT MYPGS 82 : NBMETH=MYPGS.LISPG(/1) 83 : DO 1 IMETH=1,NBMETH 84 : PGCOUR=MYPGS.LISPG(IMETH) 85 : SEGACT PGCOUR 86 : WRITE(IOIMP,*) 'Méthode d''intégration ',IMETH,' : ' 87 : $ ,PGCOUR.NOMPG 88 : IF (IMPR.GT.1) THEN 88 : CALL PRPG(PGCOUR,IMPR,IRET) 89 : IF (IRET.NE.0) GOTO 9999 90 : ENDIF 91 : NORDM=PGCOUR.NORDPG 92 : NESREF=PGCOUR.XCOPG(/1) 93 : NPOGAU=PGCOUR.XCOPG(/2) 94 : DO 2 IORD=0,NORDM+1 95 : IF (PGCOUR.FORLPG.EQ.'SEGMENT') THEN 96 : NDESP=NESREF 97 : SEGINI MYMONO 98 : MYMONO.EXPOSA(1)=IORD 99 : I=MYMONO.EXPOSA(1) 100 : IF (PGCOUR.TYPMPG.EQ.'GAUSS'.OR. 101 : $ PGCOUR.TYPMPG.EQ.'NEWTON-COTES') THEN 102 : IF (MOD(I,2).EQ.0) THEN 103 : SOLEX=2.D0/(DBLE(I)+1.D0) 104 : ELSE 105 : SOLEX=0.D0 106 : ENDIF 107 : ELSEIF (PGCOUR.TYPMPG.EQ.'GAUSS-JACOBI10') THEN 108 : SOLEX=(1.D0/(DBLE(I+1))) 109 : $ -(1.D0/(DBLE(I+2))) 110 : ELSEIF (PGCOUR.TYPMPG.EQ.'GAUSS-JACOBI20') THEN 111 : SOLEX=(1.D0/(DBLE(I+1))) 112 : $ -(2.D0/(DBLE(I+2))) 113 : $ +(1.D0/(DBLE(I+3))) 114 : ELSE 115 : WRITE(IOIMP,*) 'Type de méthode non reconnue' 116 : GOTO 9999 117 : ENDIF 118 : SOLAPP=0.D0 119 : DO 22 IPG=1,NPOGAU 120 : XPRMON=1.D0 121 : DO 222 ICOOR=1,NDESP 122 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 123 : $ **MYMONO.EXPOSA(ICOOR)) 124 : 222 CONTINUE 125 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 126 : 22 CONTINUE 127 : SEGSUP MYMONO 128 : SERR=XERREL(SOLEX,SOLAPP) 129 : IF (IMPR.GT.3) THEN 130 : WRITE(IOIMP,*) 'Monome : ksi**',I,' :' 131 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 132 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 133 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 134 : ENDIF 135 : IF (SERR.GT.2.D-14.AND.IORD.LE.NORDM) THEN 136 : GOTO 9998 137 : ENDIF 138 : ELSEIF (PGCOUR.FORLPG.EQ.'TRIANGLE') THEN 139 : DO 26 IEXP=IORD,0,-1 140 : NDESP=NESREF 141 : SEGINI MYMONO 142 : MYMONO.EXPOSA(1)=IEXP 143 : MYMONO.EXPOSA(2)=IORD-IEXP 144 : I=MYMONO.EXPOSA(1) 145 : J=MYMONO.EXPOSA(2) 146 : SOLEX=(DFACT(I)*DFACT(J))/DFACT(I+J+2) 147 : SOLAPP=0.D0 148 : DO 262 IPG=1,NPOGAU 149 : XPRMON=1.D0 150 : DO 2622 ICOOR=1,NDESP 151 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 152 : $ **MYMONO.EXPOSA(ICOOR)) 153 : 2622 CONTINUE 154 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 155 : 262 CONTINUE 156 : SEGSUP MYMONO 157 : SERR=XERREL(SOLEX,SOLAPP) 158 : IF (IMPR.GT.3) THEN 159 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J,' :' 160 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 161 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 162 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 163 : ENDIF 164 : IF (SERR.GT.1.D-15.AND.IORD.LE.NORDM) THEN 165 : GOTO 9998 166 : ENDIF 167 : 26 CONTINUE 168 : ELSEIF (PGCOUR.FORLPG.EQ.'CARRE') THEN 169 : DO 24 IEXP=IORD,0,-1 170 : NDESP=NESREF 171 : SEGINI MYMONO 172 : MYMONO.EXPOSA(1)=IEXP 173 : MYMONO.EXPOSA(2)=IORD-IEXP 174 : I=MYMONO.EXPOSA(1) 175 : J=MYMONO.EXPOSA(2) 176 : IF (MOD(I,2).EQ.0.AND.MOD(J,2).EQ.0) THEN 177 : SOLEX=4.D0/((DBLE(I)+1.D0)*(DBLE(J)+1.D0)) 178 : ELSE 179 : SOLEX=0.D0 180 : ENDIF 181 : SOLAPP=0.D0 182 : DO 242 IPG=1,NPOGAU 183 : XPRMON=1.D0 184 : DO 2422 ICOOR=1,NDESP 185 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 186 : $ **MYMONO.EXPOSA(ICOOR)) 187 : 2422 CONTINUE 188 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 189 : 242 CONTINUE 190 : SEGSUP MYMONO 191 : SERR=XERREL(SOLEX,SOLAPP) 192 : IF (IMPR.GT.3) THEN 193 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J,' :' 194 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 195 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 196 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 197 : ENDIF 198 : IF (SERR.GT.1.D-15.AND.IORD.LE.NORDM) THEN 199 : GOTO 9998 200 : ENDIF 201 : 24 CONTINUE 202 : ELSEIF (PGCOUR.FORLPG.EQ.'TETRAEDRE') THEN 203 : DO 28 IEXP1=IORD,0,-1 204 : DO 282 IEXP2=IORD-IEXP1,0,-1 205 : NDESP=NESREF 206 : SEGINI MYMONO 207 : MYMONO.EXPOSA(1)=IEXP1 208 : MYMONO.EXPOSA(2)=IEXP2 209 : MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2 210 : I=MYMONO.EXPOSA(1) 211 : J=MYMONO.EXPOSA(2) 212 : K=MYMONO.EXPOSA(3) 213 : SOLEX=(DFACT(I)*DFACT(J)*DFACT(K))/DFACT(I+J+K+3) 214 : SOLAPP=0.D0 215 : DO 2822 IPG=1,NPOGAU 216 : XPRMON=1.D0 217 : DO 28222 ICOOR=1,NDESP 218 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 219 : $ **MYMONO.EXPOSA(ICOOR)) 220 : 28222 CONTINUE 221 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 222 : 2822 CONTINUE 223 : SEGSUP MYMONO 224 : SERR=XERREL(SOLEX,SOLAPP) 225 : IF (IMPR.GT.3) THEN 226 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J 227 : $ ,' :',' dzeta**',K,' :' 228 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 229 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 230 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 231 : ENDIF 232 : IF (SERR.GT.2.D-15.AND.IORD.LE.NORDM) THEN 233 : GOTO 9998 234 : ENDIF 235 : 282 CONTINUE 236 : 28 CONTINUE 237 : ELSEIF (PGCOUR.FORLPG.EQ.'PYRAMIDE') THEN 238 : NDESP=NESREF 239 : SEGINI MYMONO 240 : MYMONO.EXPOSA(1)=0 241 : MYMONO.EXPOSA(2)=0 242 : MYMONO.EXPOSA(3)=IORD 243 : I=MYMONO.EXPOSA(1) 244 : J=MYMONO.EXPOSA(2) 245 : K=MYMONO.EXPOSA(3) 246 : SOLEX=2.D0*((1.D0/(DBLE(IORD+3))) 247 : $ +(-2.D0/(DBLE(IORD+2))) 248 : $ +(1.D0/(DBLE(IORD+1)))) 249 : SOLAPP=0.D0 250 : DO 2922 IPG=1,NPOGAU 251 : XPRMON=1.D0 252 : DO 29222 ICOOR=1,NDESP 253 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 254 : $ **MYMONO.EXPOSA(ICOOR)) 255 : 29222 CONTINUE 256 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 257 : 2922 CONTINUE 258 : SEGSUP MYMONO 259 : SERR=XERREL(SOLEX,SOLAPP) 260 : IF (IMPR.GT.3) THEN 261 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J 262 : $ ,' :',' dzeta**',K,' :' 263 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 264 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 265 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 266 : ENDIF 267 : IF (SERR.GT.1.D-8.AND.IORD.LE.NORDM) THEN 268 : GOTO 9998 269 : ENDIF 270 : ELSEIF (PGCOUR.FORLPG.EQ.'PRISME') THEN 271 : DO 30 IEXP1=IORD,0,-1 272 : DO 302 IEXP2=IORD-IEXP1,0,-1 273 : NDESP=NESREF 274 : SEGINI MYMONO 275 : MYMONO.EXPOSA(1)=IEXP1 276 : MYMONO.EXPOSA(2)=IEXP2 277 : MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2 278 : I=MYMONO.EXPOSA(1) 279 : J=MYMONO.EXPOSA(2) 280 : K=MYMONO.EXPOSA(3) 281 : IF (MOD(K,2).EQ.0) THEN 282 : SOLEX=((DFACT(I)*DFACT(J))/DFACT(I+J+2))* 283 : $ (2.D0/(DBLE(K)+1.D0)) 284 : ELSE 285 : SOLEX=0.D0 286 : ENDIF 287 : SOLAPP=0.D0 288 : DO 3022 IPG=1,NPOGAU 289 : XPRMON=1.D0 290 : DO 30222 ICOOR=1,NDESP 291 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 292 : $ **MYMONO.EXPOSA(ICOOR)) 293 : 30222 CONTINUE 294 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 295 : 3022 CONTINUE 296 : SEGSUP MYMONO 297 : SERR=XERREL(SOLEX,SOLAPP) 298 : IF (IMPR.GT.3) THEN 299 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J 300 : $ ,' :',' dzeta**',K,' :' 301 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 302 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 303 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 304 : ENDIF 305 : IF (SERR.GT.1.D-15.AND.IORD.LE.NORDM) THEN 306 : GOTO 9998 307 : ENDIF 308 : 302 CONTINUE 309 : 30 CONTINUE 310 : ELSEIF (PGCOUR.FORLPG.EQ.'CUBE') THEN 311 : DO 32 IEXP1=IORD,0,-1 312 : DO 322 IEXP2=IORD-IEXP1,0,-1 313 : NDESP=NESREF 314 : SEGINI MYMONO 315 : MYMONO.EXPOSA(1)=IEXP1 316 : MYMONO.EXPOSA(2)=IEXP2 317 : MYMONO.EXPOSA(3)=IORD-IEXP1-IEXP2 318 : I=MYMONO.EXPOSA(1) 319 : J=MYMONO.EXPOSA(2) 320 : K=MYMONO.EXPOSA(3) 321 : IF (MOD(I,2).EQ.0.AND.MOD(J,2).EQ.0 322 : $ .AND.MOD(K,2).EQ.0) THEN 323 : SOLEX=8.D0/((DBLE(I)+1.D0)*(DBLE(J)+1.D0) 324 : $ *(DBLE(K)+1.D0)) 325 : ELSE 326 : SOLEX=0.D0 327 : ENDIF 328 : SOLAPP=0.D0 329 : DO 3222 IPG=1,NPOGAU 330 : XPRMON=1.D0 331 : DO 32222 ICOOR=1,NDESP 332 : XPRMON=XPRMON*(PGCOUR.XCOPG(ICOOR,IPG) 333 : $ **MYMONO.EXPOSA(ICOOR)) 334 : 32222 CONTINUE 335 : SOLAPP=SOLAPP+(PGCOUR.XPOPG(IPG)*XPRMON) 336 : 3222 CONTINUE 337 : SEGSUP MYMONO 338 : SERR=XERREL(SOLEX,SOLAPP) 339 : IF (IMPR.GT.3) THEN 340 : WRITE(IOIMP,*) 'Monome : ksi**',I,' eta**',J 341 : $ ,' :',' dzeta**',K,' :' 342 : WRITE(IOIMP,*) ' Solution exacte : ',SOLEX 343 : WRITE(IOIMP,*) ' Solution appr. : ',SOLAPP 344 : WRITE(IOIMP,*) ' Erreur relative : ',SERR 345 : ENDIF 346 : IF (SERR.GT.2.D-15.AND.IORD.LE.NORDM) THEN 347 : GOTO 9998 348 : ENDIF 349 : 322 CONTINUE 350 : 32 CONTINUE 351 : ELSE 352 : WRITE(IOIMP,*) PGCOUR.FORLPG,' n''est pas une forme' 353 : WRITE(IOIMP,*) 'd''élément reconnue.' 354 : GOTO 9999 355 : ENDIF 356 : 2 CONTINUE 357 : SEGDES PGCOUR 358 : WRITE(IOIMP,*) 'Ok' 359 : 1 CONTINUE 360 : SEGDES MYPGS 361 : * 362 : * Normal termination 363 : * 364 : IRET=0 364 : CALL ECRLOG(.TRUE.) 365 : RETURN 366 : * 367 : * Format handling 368 : * 369 : * 370 : * Error handling 371 : * 372 : 9998 CONTINUE 373 : WRITE(IOIMP,*) 'SERR=',SERR,' => Pb. ds la méthode !!' 374 : 9999 CONTINUE 375 : IRET=1 376 : WRITE(IOIMP,*) 'An error was detected in subroutine testpg' 376 : CALL ECRLOG(.FALSE.) 377 : RETURN 378 : * 379 : * End of subroutine testpg 380 : * 381 : END 382 : 383 :
© Cast3M 2003 - All rights reserved.
Disclaimer