1 : C INGATR SOURCE GOUNAND 05/12/21 21:32:26 5281 2 : SUBROUTINE INGATR(MYPGS,IMPR,IRET) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : INGATR 7 : C PROJET : Noyau linéaire NLIN 8 : C DESCRIPTION : Remplit le segment des méthodes d'intégration 9 : C avec des méthodes d'intégration numérique de cubature 10 : C type Gauss pour le triangle (ordre 1 à 7). 11 : C 12 : C REFERENCES : Le site de Cools (avec 32 chiffres sign.) 13 : C (essentiellement Stroud et al.) dont on reprend la 14 : C nomenclature... 15 : C LANGAGE : ESOPE 16 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 17 : C mél : gounand@semt2.smts.cea.fr 18 : C*********************************************************************** 19 : C APPELES : INIPG, GTSINO, GTRO3I, FIPG, CPROPG 20 : C APPELE PAR : INPGS 21 : C*********************************************************************** 22 : C ENTREES : - 23 : C ENTREES/SORTIES : MYPGS (actif en *MOD) 24 : C SORTIES : - 25 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 26 : C*********************************************************************** 27 : C VERSION : v1, 20/10/99, version initiale 28 : C HISTORIQUE : v1, 20/10/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 : POINTEUR PGPRO1.POGAU 53 : POINTEUR PGPRO2.POGAU 54 : * 55 : INTEGER IMPR,IRET 56 : * On travaille en coordonnées barycentriques. 57 : * D'après le Dhatt et Touzot 58 : * L1 = 1 - \ksi - \eta 59 : * L2 = \ksi 60 : * L3 = \eta 61 : INTEGER DIMSRF 62 : PARAMETER(DIMSRF=2) 63 : REAL*8 XCOR(DIMSRF+1) 64 : REAL*8 POIDS 65 : * 66 : * Générateurs pour la cubature de degré 1 à 1 point : GAT2-1-1 : 67 : * - [ Fully symmetric ] 68 : REAL*8 X1D1,Y1D1,P1D1 69 : PARAMETER (X1D1=1.D0/3.D0) 70 : PARAMETER (Y1D1=1.D0/3.D0) 71 : PARAMETER (P1D1=0.5D0) 72 : * 73 : * Générateurs pour la cubature de degré 2 à 3 points : GAT2-2-3A : 74 : * - [ Fully symmetric ] 75 : REAL*8 X1D2,Y1D2,P1D2 76 : PARAMETER (X1D2=0.5D0) 77 : PARAMETER (Y1D2=0.5D0) 78 : PARAMETER (P1D2=1.D0/6.D0) 79 : * 80 : * Pour le degré 3, méthode produit conique à 4 points : GPT2-3-4 81 : * 82 : * Générateurs pour la cubature de degré 4 à 6 points : GAT2-4-6A : 83 : * - [ Fully symmetric ] 84 : REAL*8 X1D4,Y1D4,P1D4 85 : PARAMETER (X1D4=0.0915762135097707434595714634022015D0) 86 : PARAMETER (Y1D4=0.0915762135097707434595714634022015D0) 87 : PARAMETER (P1D4=0.0549758718276609338191631624501052D0) 88 : * - [ Fully symmetric ] 89 : REAL*8 X2D4,Y2D4,P2D4 90 : PARAMETER (X2D4=0.445948490915964886318329253883051D0) 91 : PARAMETER (Y2D4=0.445948490915964886318329253883051D0) 92 : PARAMETER (P2D4=0.111690794839005732847503504216561D0) 93 : * 94 : * Générateurs pour la cubature de degré 5 à 7 points : GAT2-5-7 : 95 : * - [ Fully symmetric ] 96 : REAL*8 X1D5,Y1D5,P1D5 97 : PARAMETER (X1D5=1.D0/3.D0) 98 : PARAMETER (Y1D5=1.D0/3.D0) 99 : PARAMETER (P1D5=0.1125D0) 100 : * - [ Fully symmetric ] 101 : REAL*8 X2D5,Y2D5,P2D5 102 : PARAMETER (X2D5=0.101286507323456338800987361915123D0) 103 : PARAMETER (Y2D5=0.101286507323456338800987361915123D0) 104 : PARAMETER (P2D5=0.0629695902724135762978419727500906D0) 105 : * - [ Fully symmetric ] 106 : REAL*8 X3D5,Y3D5,P3D5 107 : PARAMETER (X3D5=0.470142064105115089770441209513447D0) 108 : PARAMETER (Y3D5=0.470142064105115089770441209513447D0) 109 : PARAMETER (P3D5=0.0661970763942530903688246939165759D0) 110 : * 111 : * Pas de degré 6, même nb. de points que degré 7 avec la 112 : * qualité PI (coeffs positifs, points à l'intérieur). 113 : * 114 : * Générateurs pour la cubature de degré 7 à 12 points : GAT2-7-12 : 115 : * - [ Ro3 invariant ] 116 : REAL*8 X1D7,Y1D7,P1D7 117 : PARAMETER (X1D7=0.0623822650944021181736830009963499D0) 118 : PARAMETER (Y1D7=0.0675178670739160854425571310508685D0) 119 : PARAMETER (P1D7=0.0265170281574362514287541804607391D0) 120 : * - [ Ro3 invariant ] 121 : REAL*8 X2D7,Y2D7,P2D7 122 : PARAMETER (X2D7=0.0552254566569266117374791902756449D0) 123 : PARAMETER (Y2D7=0.321502493851981822666307849199202D0) 124 : PARAMETER (P2D7=0.0438814087144460550367699031392875D0) 125 : * - [ Ro3 invariant ] 126 : REAL*8 X3D7,Y3D7,P3D7 127 : PARAMETER (X3D7=0.0343243029450971464696306424839376D0) 128 : PARAMETER (Y3D7=0.660949196186735657611980310197799D0) 129 : PARAMETER (P3D7=0.0287750427849815857384454969002185D0) 130 : * - [ Ro3 invariant ] 131 : REAL*8 X4D7,Y4D7,P4D7 132 : PARAMETER (X4D7=0.515842334353591779257463386826430D0) 133 : PARAMETER (Y4D7=0.277716166976391782569581871393723D0) 134 : PARAMETER (P4D7=0.0674931870098027744626970861664214D0) 135 : * 136 : INTEGER NOPG,IPG 137 : * 138 : * Executable statements 139 : * 140 : IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingatr' 141 : * 142 : * Méthode de nom : NCT2-1-3 143 : * Sur un triangle : cubature d'ordre 1 à 3 points 144 : * espace de référence de dimension 2 145 : * 146 : * In INIPG : SEGINI PGCOUR 146 : CALL INIPG('NCT2-1-3','NEWTON-COTES','TRIANGLE', 147 : $ 1,3,2, 148 : $ PGCOUR, 149 : $ IMPR,IRET) 150 : IF (IRET.NE.0) GOTO 9999 151 : NOPG=0 152 : XCOR(1)=1.D0 153 : XCOR(2)=0.D0 154 : XCOR(3)=0.D0 154 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET) 155 : IF (IRET.NE.0) GOTO 9999 156 : SEGDES PGCOUR 157 : MYPGS.LISPG(**)=PGCOUR 158 : * 159 : * Méthode de nom : NCT2-3-7 160 : * Sur un triangle : cubature d'ordre 3 à 7 points 161 : * espace de référence de dimension 2 162 : * 163 : * In INIPG : SEGINI PGCOUR 163 : CALL INIPG('NCT2-3-7','NEWTON-COTES','TRIANGLE', 164 : $ 3,7,2, 165 : $ PGCOUR, 166 : $ IMPR,IRET) 167 : IF (IRET.NE.0) GOTO 9999 168 : NOPG=0 169 : XCOR(1)=1.D0 170 : XCOR(2)=0.D0 171 : XCOR(3)=0.D0 172 : POIDS=0.025D0 172 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET) 173 : IF (IRET.NE.0) GOTO 9999 174 : XCOR(1)=0.5D0 175 : XCOR(2)=0.5D0 176 : XCOR(3)=0.D0 177 : POIDS=0.2D0/3.D0 177 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET) 178 : IF (IRET.NE.0) GOTO 9999 179 : XCOR(1)=1.D0/3.D0 180 : XCOR(2)=1.D0/3.D0 181 : XCOR(3)=1-XCOR(1)-XCOR(2) 182 : POIDS=0.225D0 182 : CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET) 183 : IF (IRET.NE.0) GOTO 9999 184 : 185 : SEGDES PGCOUR 186 : MYPGS.LISPG(**)=PGCOUR 187 : * 188 : * Méthode de nom : GAT2-1-1 189 : * Sur un triangle : cubature d'ordre 1 à 1 point 190 : * espace de référence de dimension 2 191 : * 192 : * In INIPG : SEGINI PGCOUR 192 : CALL INIPG('GAT2-1-1','GAUSS','TRIANGLE', 193 : $ 1,1,2, 194 : $ PGCOUR, 195 : $ IMPR,IRET) 196 : IF (IRET.NE.0) GOTO 9999 197 : NOPG=0 198 : XCOR(1)=X1D1 199 : XCOR(2)=Y1D1 200 : XCOR(3)=1-XCOR(1)-XCOR(2) 200 : CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D1,IMPR,IRET) 201 : IF (IRET.NE.0) GOTO 9999 202 : SEGDES PGCOUR 203 : MYPGS.LISPG(**)=PGCOUR 204 : * 205 : * Méthode de nom : GAT2-2-3A 206 : * Sur un triangle : cubature d'ordre 2 à 3 points 207 : * espace de référence de dimension 2 208 : * 209 : * In INIPG : SEGINI PGCOUR 209 : CALL INIPG('GAT2-2-3A','GAUSS','TRIANGLE', 210 : $ 2,3,2, 211 : $ PGCOUR, 212 : $ IMPR,IRET) 213 : IF (IRET.NE.0) GOTO 9999 214 : NOPG=0 215 : XCOR(1)=X1D2 216 : XCOR(2)=Y1D2 217 : XCOR(3)=1-XCOR(1)-XCOR(2) 217 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET) 218 : IF (IRET.NE.0) GOTO 9999 219 : SEGDES PGCOUR 220 : MYPGS.LISPG(**)=PGCOUR 221 : * 222 : * Méthode de nom : GPT2-3-4 223 : * Sur un triangle : méthode gauss-produit conique d'ordre 3 à 4 points 224 : * espace de référence de dimension 2 225 : * 226 : * In INIPG : SEGINI PGCOUR 226 : CALL INIPG('GPT2-3-4','GAUSS-PRODUIT','TRIANGLE', 227 : $ 3,4,2, 228 : $ PGCOUR, 229 : $ IMPR,IRET) 230 : IF (IRET.NE.0) GOTO 9999 230 : CALL FIPG('GAC1-3-2',MYPGS,PGPRO1,IMPR,IRET) 231 : IF (IRET.NE.0) GOTO 9999 232 : * 233 : * Il faut remettre les points de Gauss-Legendre sur [0,1] 234 : * à la place de [-1,1] 235 : * 236 : SEGINI,PGPRO2=PGPRO1 237 : DO 1 IPG=1,PGPRO2.XCOPG(/2) 238 : PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0 239 : PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0 240 : 1 CONTINUE 241 : SEGDES PGPRO2 242 : PGPRO1=PGPRO2 242 : CALL FIPG('GJ10-3-2',MYPGS,PGPRO2,IMPR,IRET) 243 : IF (IRET.NE.0) GOTO 9999 243 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 244 : IF (IRET.NE.0) GOTO 9999 245 : SEGDES PGCOUR 246 : MYPGS.LISPG(**)=PGCOUR 247 : * 248 : * Méthode de nom : GAT2-4-6A 249 : * Sur un triangle : cubature d'ordre 4 à 6 points 250 : * espace de référence de dimension 2 251 : * 252 : * In INIPG : SEGINI PGCOUR 252 : CALL INIPG('GAT2-4-6A','GAUSS','TRIANGLE', 253 : $ 4,6,2, 254 : $ PGCOUR, 255 : $ IMPR,IRET) 256 : IF (IRET.NE.0) GOTO 9999 257 : NOPG=0 258 : XCOR(1)=X1D4 259 : XCOR(2)=Y1D4 260 : XCOR(3)=1-XCOR(1)-XCOR(2) 260 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D4,IMPR,IRET) 261 : IF (IRET.NE.0) GOTO 9999 262 : XCOR(1)=X2D4 263 : XCOR(2)=Y2D4 264 : XCOR(3)=1-XCOR(1)-XCOR(2) 264 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D4,IMPR,IRET) 265 : IF (IRET.NE.0) GOTO 9999 266 : SEGDES PGCOUR 267 : MYPGS.LISPG(**)=PGCOUR 268 : * 269 : * Méthode de nom : GAT2-5-7 270 : * Sur un triangle : cubature d'ordre 5 à 7 points 271 : * espace de référence de dimension 2 272 : * 273 : * In INIPG : SEGINI PGCOUR 273 : CALL INIPG('GAT2-5-7','GAUSS','TRIANGLE', 274 : $ 5,7,2, 275 : $ PGCOUR, 276 : $ IMPR,IRET) 277 : IF (IRET.NE.0) GOTO 9999 278 : NOPG=0 279 : XCOR(1)=X1D5 280 : XCOR(2)=Y1D5 281 : XCOR(3)=1-XCOR(1)-XCOR(2) 281 : CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET) 282 : IF (IRET.NE.0) GOTO 9999 283 : XCOR(1)=X2D5 284 : XCOR(2)=Y2D5 285 : XCOR(3)=1-XCOR(1)-XCOR(2) 285 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET) 286 : IF (IRET.NE.0) GOTO 9999 287 : XCOR(1)=X3D5 288 : XCOR(2)=Y3D5 289 : XCOR(3)=1-XCOR(1)-XCOR(2) 289 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P3D5,IMPR,IRET) 290 : IF (IRET.NE.0) GOTO 9999 291 : SEGDES PGCOUR 292 : MYPGS.LISPG(**)=PGCOUR 293 : * 294 : * Méthode de nom : GPT2-5-9 295 : * Sur un triangle : méthode gauss-produit conique d'ordre 5 à 9 points 296 : * espace de référence de dimension 2 297 : * 298 : * In INIPG : SEGINI PGCOUR 298 : CALL INIPG('GPT2-5-9','GAUSS-PRODUIT','TRIANGLE', 299 : $ 5,9,2, 300 : $ PGCOUR, 301 : $ IMPR,IRET) 302 : IF (IRET.NE.0) GOTO 9999 302 : CALL FIPG('GAC1-5-3',MYPGS,PGPRO1,IMPR,IRET) 303 : IF (IRET.NE.0) GOTO 9999 304 : * 305 : * Il faut remettre les points de Gauss-Legendre sur [0,1] 306 : * à la place de [-1,1] 307 : * 308 : SEGINI,PGPRO2=PGPRO1 309 : DO 3 IPG=1,PGPRO2.XCOPG(/2) 310 : PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0 311 : PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0 312 : 3 CONTINUE 313 : SEGDES PGPRO2 314 : PGPRO1=PGPRO2 314 : CALL FIPG('GJ10-5-3',MYPGS,PGPRO2,IMPR,IRET) 315 : IF (IRET.NE.0) GOTO 9999 315 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 316 : IF (IRET.NE.0) GOTO 9999 317 : SEGDES PGCOUR 318 : MYPGS.LISPG(**)=PGCOUR 319 : * 320 : * Méthode de nom : GAT2-7-12 321 : * Sur un triangle : cubature d'ordre 7 à 12 points 322 : * espace de référence de dimension 2 323 : * 324 : * In INIPG : SEGINI PGCOUR 324 : CALL INIPG('GAT2-7-12','GAUSS','TRIANGLE', 325 : $ 7,12,2, 326 : $ PGCOUR, 327 : $ IMPR,IRET) 328 : IF (IRET.NE.0) GOTO 9999 329 : NOPG=0 330 : XCOR(1)=X1D7 331 : XCOR(2)=Y1D7 332 : XCOR(3)=1-XCOR(1)-XCOR(2) 332 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET) 333 : IF (IRET.NE.0) GOTO 9999 334 : XCOR(1)=X2D7 335 : XCOR(2)=Y2D7 336 : XCOR(3)=1-XCOR(1)-XCOR(2) 336 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET) 337 : IF (IRET.NE.0) GOTO 9999 338 : XCOR(1)=X3D7 339 : XCOR(2)=Y3D7 340 : XCOR(3)=1-XCOR(1)-XCOR(2) 340 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P3D7,IMPR,IRET) 341 : IF (IRET.NE.0) GOTO 9999 342 : XCOR(1)=X4D7 343 : XCOR(2)=Y4D7 344 : XCOR(3)=1-XCOR(1)-XCOR(2) 344 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P4D7,IMPR,IRET) 345 : IF (IRET.NE.0) GOTO 9999 346 : SEGDES PGCOUR 347 : MYPGS.LISPG(**)=PGCOUR 348 : * 349 : * Méthode de nom : GPT2-7-16 350 : * Sur un triangle : méthode gauss-produit conique d'ordre 7 à 16 points 351 : * espace de référence de dimension 2 352 : * 353 : * In INIPG : SEGINI PGCOUR 353 : CALL INIPG('GPT2-7-16','GAUSS-PRODUIT','TRIANGLE', 354 : $ 7,16,2, 355 : $ PGCOUR, 356 : $ IMPR,IRET) 357 : IF (IRET.NE.0) GOTO 9999 357 : CALL FIPG('GAC1-7-4',MYPGS,PGPRO1,IMPR,IRET) 358 : IF (IRET.NE.0) GOTO 9999 359 : * 360 : * Il faut remettre les points de Gauss-Legendre sur [0,1] 361 : * à la place de [-1,1] 362 : * 363 : SEGINI,PGPRO2=PGPRO1 364 : DO 5 IPG=1,PGPRO2.XCOPG(/2) 365 : PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0 366 : PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0 367 : 5 CONTINUE 368 : SEGDES PGPRO2 369 : PGPRO1=PGPRO2 369 : CALL FIPG('GJ10-7-4',MYPGS,PGPRO2,IMPR,IRET) 370 : IF (IRET.NE.0) GOTO 9999 370 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 371 : IF (IRET.NE.0) GOTO 9999 372 : SEGDES PGCOUR 373 : MYPGS.LISPG(**)=PGCOUR 374 : * 375 : * Normal termination 376 : * 377 : IRET=0 378 : RETURN 379 : * 380 : * Format handling 381 : * 382 : * 383 : * Error handling 384 : * 385 : 9999 CONTINUE 386 : IRET=1 387 : WRITE(IOIMP,*) 'An error was detected in subroutine ingatr' 388 : RETURN 389 : * 390 : * End of subroutine INGATR 391 : * 392 : END 393 : 394 : 395 :
© Cast3M 2003 - All rights reserved.
Disclaimer