1 : C INGAQU SOURCE PV 07/11/23 21:17:18 5978 2 : SUBROUTINE INGAQU(MYPGS,IMPR,IRET) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : INGAQU 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 carré (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, GCSINO, GCPASY, GCROIN, GCCESY, GCRESY 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 : * 53 : integer PGPRO1 54 : INTEGER IMPR,IRET 55 : * 56 : INTEGER DIMSRF 57 : PARAMETER(DIMSRF=2) 58 : REAL*8 XCOR(DIMSRF) 59 : * 60 : * Générateurs pour la cubature de degré 1 à 1 point : GAC2-1-1 : 61 : * - [ Origin ] 62 : REAL*8 X0D1,Y0D1,P0D1 63 : PARAMETER (X0D1= 0.D0) 64 : PARAMETER (Y0D1= 0.D0) 65 : PARAMETER (P0D1= 4.D0) 66 : * 67 : * Générateurs pour la cubature de degré 2 à 3 points : GAC2-2-3 : 68 : * - [ Partial symmetry ] (\sqrt{2./3.} ; 0.) ; (4./3.) 69 : REAL*8 X1D2,Y1D2,P1D2 70 : PARAMETER (X1D2= 0.81649658092772603273242802490196D0) 71 : PARAMETER (Y1D2= 0.D0) 72 : PARAMETER (P1D2= 4.D0/3.D0) 73 : * - [ Partial symmetry ] (-\sqrt{1./6.} ; \sqrt{1./2.}) ; (4./3.) 74 : REAL*8 X2D2,Y2D2,P2D2 75 : PARAMETER (X2D2=-0.40824829046386301636621401245098D0) 76 : PARAMETER (Y2D2= 0.70710678118654752440084436210485D0) 77 : PARAMETER (P2D2= 4.D0/3.D0) 78 : * 79 : * Générateurs pour la cubature de degré 3 à 4 points : GAC2-3-4A : 80 : * - [ Fully symmetric ] (\sqrt{2./3.} ; 0.) ; (1.) 81 : REAL*8 X1D3,Y1D3,P1D3 82 : PARAMETER (X1D3= 0.816496580927726032732428024901963D0) 83 : PARAMETER (Y1D3= 0.D0) 84 : PARAMETER (P1D3= 1.D0) 85 : * 86 : * Générateurs pour la cubature de degré 4 à 6 points : GAC2-4-6C : 87 : * - [ Origin ] 88 : REAL*8 X0D4,Y0D4,P0D4 89 : PARAMETER (X0D4= 0.D0) 90 : PARAMETER (Y0D4= 0.D0) 91 : PARAMETER (P0D4= 1.14285714285714285714285714285714D0) 92 : * - [ Partial symmetry ] 93 : REAL*8 X1D4,Y1D4,P1D4 94 : PARAMETER (X1D4= 0.966091783079295904914577610477984D0) 95 : PARAMETER (Y1D4= 0.D0) 96 : PARAMETER (P1D4= 0.439560439560439560439560439560439D0) 97 : * - [ Partial symmetry ] 98 : REAL*8 X2D4,Y2D4,P2D4 99 : PARAMETER (X2D4= 0.455603727836192842249584745989760D0) 100 : PARAMETER (Y2D4= 0.851914653304600491301835640366669D0) 101 : PARAMETER (P2D4= 0.566072207007532105264050459451812D0) 102 : * - [ Partial symmetry ] 103 : REAL*8 X3D4,Y3D4,P3D4 104 : PARAMETER (X3D4=-0.731629951573134529368035491840613D0) 105 : PARAMETER (Y3D4= 0.630912788976754026243413202993658D0) 106 : PARAMETER (P3D4= 0.642719001783676685944740749339395D0) 107 : * 108 : * Générateurs pour la cubature de degré 5 à 7 points : GAC2-5-7A : 109 : * - [ Origin ] 110 : REAL*8 X0D5,Y0D5,P0D5 111 : PARAMETER (X0D5= 0.D0) 112 : PARAMETER (Y0D5= 0.D0) 113 : PARAMETER (P0D5= 1.14285714285714285714285714285714D0) 114 : * - [ Central symmetry ] 115 : REAL*8 X1D5,Y1D5,P1D5 116 : PARAMETER (X1D5= 0.966091783079295904914577610477984D0) 117 : PARAMETER (Y1D5= 0.D0) 118 : PARAMETER (P1D5= 0.317460317460317460317460317460317D0) 119 : * - [ Rectangular symmetry ] 120 : REAL*8 X2D5,Y2D5,P2D5 121 : PARAMETER (X2D5= 0.577350269189625764509148780501957D0) 122 : PARAMETER (Y2D5= 0.774596669241483377035853079956479D0) 123 : PARAMETER (P2D5= 0.555555555555555555555555555555555D0) 124 : * 125 : * Générateurs pour la cubature de degré 6 à 10 points : GAC2-6-10C : 126 : * - [ Partial symmetry ] 127 : REAL*8 X1D6,Y1D6,P1D6 128 : PARAMETER (X1D6= 0.836405633697625604537901440708439D0) 129 : PARAMETER (Y1D6= 0.D0) 130 : PARAMETER (P1D6= 0.455343245714173438553541721312115D0) 131 : * - [ Partial symmetry ] 132 : REAL*8 X2D6,Y2D6,P2D6 133 : PARAMETER (X2D6=-0.357460165391307184886174384445965D0) 134 : PARAMETER (Y2D6= 0.D0) 135 : PARAMETER (P2D6= 0.827395973202965492356696938299033D0) 136 : * - [ Partial symmetry ] 137 : REAL*8 X3D6,Y3D6,P3D6 138 : PARAMETER (X3D6= 0.872101531193130590032875340374412D0) 139 : PARAMETER (Y3D6= 0.888764014654764533853170471609380D0) 140 : PARAMETER (P3D6= 0.144000884599645384928603086560623D0) 141 : * - [ Partial symmetry ] 142 : REAL*8 X4D6,Y4D6,P4D6 143 : PARAMETER (X4D6= 0.305985162155426661489382735682569D0) 144 : PARAMETER (Y4D6= 0.604857639464685027698532295335587D0) 145 : PARAMETER (P4D6= 0.668259104262665149929832544411851D0) 146 : * - [ Partial symmetry ] 147 : REAL*8 X5D6,Y5D6,P5D6 148 : PARAMETER (X5D6=-0.410270899466657838725680144555973D0) 149 : PARAMETER (Y5D6= 0.955447506641063747757528360774788D0) 150 : PARAMETER (P5D6= 0.225474004890679357411055647855252D0) 151 : * - [ Partial symmetry ] 152 : REAL*8 X6D6,Y6D6,P6D6 153 : PARAMETER (X6D6=-0.872869311156879339251259208496329D0) 154 : PARAMETER (Y6D6= 0.565459993438754045894375389139462D0) 155 : PARAMETER (P6D6= 0.320896396788440642275389391366698D0) 156 : * 157 : * Générateurs pour la cubature de degré 7 à 12 points : GAC2-7-12A : 158 : * - [ Fully symmetric ] 159 : REAL*8 X1D7,Y1D7,P1D7 160 : PARAMETER (X1D7= 0.925820099772551461566566776583999D0) 161 : PARAMETER (Y1D7= 0.D0) 162 : PARAMETER (P1D7= 0.241975308641975308641975308641975D0) 163 : * - [ Fully symmetric ] 164 : REAL*8 X2D7,Y2D7,P2D7 165 : PARAMETER (X2D7= 0.380554433208315656379106359086394D0) 166 : PARAMETER (Y2D7= 0.380554433208315656379106359086394D0) 167 : PARAMETER (P2D7= 0.520592916667394457139919432046731D0) 168 : * - [ Fully symmetric ] 169 : REAL*8 X3D7,Y3D7,P3D7 170 : PARAMETER (X3D7= 0.805979782918598743707856181350744D0) 171 : PARAMETER (Y3D7= 0.805979782918598743707856181350744D0) 172 : PARAMETER (P3D7= 0.237431774690630234218105259311293D0) 173 : * 174 : INTEGER NOPG 175 : * 176 : * Executable statements 177 : * 178 : IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingaqu' 179 : * 180 : * Méthode de nom : NCC2-1-4 181 : * Sur un carré : cubature d'ordre 1 à 4 points 182 : * espace de référence de dimension 2 183 : * 184 : * In INIPG : SEGINI PGCOUR 184 : CALL INIPG('NCC2-1-4','NEWTON-COTES','CARRE', 185 : $ 1,4,2, 186 : $ PGCOUR, 187 : $ IMPR,IRET) 188 : IF (IRET.NE.0) GOTO 9999 189 : NOPG=0 190 : XCOR(1)=1.D0 191 : XCOR(2)=1.D0 191 : CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 192 : IF (IRET.NE.0) GOTO 9999 193 : SEGDES PGCOUR 194 : MYPGS.LISPG(**)=PGCOUR 195 : * 196 : * Méthode de nom : NCC2-3-9 197 : * Sur un carré : cubature d'ordre 3 à 9 points 198 : * espace de référence de dimension 2 199 : * 200 : * In INIPG : SEGINI PGCOUR 200 : CALL INIPG('NCC2-3-9','NEWTON-COTES','CARRE', 201 : $ 3,9,2, 202 : $ PGCOUR, 203 : $ IMPR,IRET) 204 : IF (IRET.NE.0) GOTO 9999 204 : CALL FIPG('NCC1-3-3',MYPGS,PGPRO1,IMPR,IRET) 205 : IF (IRET.NE.0) GOTO 9999 205 : CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET) 206 : IF (IRET.NE.0) GOTO 9999 207 : SEGDES PGCOUR 208 : MYPGS.LISPG(**)=PGCOUR 209 : * 210 : * Méthode de nom : GAC2-1-1 211 : * Sur un carré : cubature d'ordre 1 à 1 point 212 : * espace de référence de dimension 2 213 : * 214 : * In INIPG : SEGINI PGCOUR 214 : CALL INIPG('GAC2-1-1','GAUSS','CARRE', 215 : $ 1,1,2, 216 : $ PGCOUR, 217 : $ IMPR,IRET) 218 : IF (IRET.NE.0) GOTO 9999 219 : NOPG=0 220 : XCOR(1)=X0D1 221 : XCOR(2)=Y0D1 221 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D1,IMPR,IRET) 222 : IF (IRET.NE.0) GOTO 9999 223 : SEGDES PGCOUR 224 : MYPGS.LISPG(**)=PGCOUR 225 : * 226 : * Méthode de nom : GAC2-2-3 227 : * Sur un carré : cubature d'ordre 2 à 3 points 228 : * espace de référence de dimension 2 229 : * 230 : * In INIPG : SEGINI PGCOUR 230 : CALL INIPG('GAC2-2-3','GAUSS','CARRE', 231 : $ 2,3,2, 232 : $ PGCOUR, 233 : $ IMPR,IRET) 234 : IF (IRET.NE.0) GOTO 9999 235 : NOPG=0 236 : XCOR(1)=X1D2 237 : XCOR(2)=Y1D2 237 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET) 238 : IF (IRET.NE.0) GOTO 9999 239 : XCOR(1)=X2D2 240 : XCOR(2)=Y2D2 240 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P2D2,IMPR,IRET) 241 : IF (IRET.NE.0) GOTO 9999 242 : SEGDES PGCOUR 243 : MYPGS.LISPG(**)=PGCOUR 244 : * 245 : * Méthode de nom : GAC2-3-4A 246 : * Sur un carré : cubature d'ordre 3 à 4 points 247 : * espace de référence de dimension 2 248 : * 249 : * In INIPG : SEGINI PGCOUR 249 : CALL INIPG('GAC2-3-4A','GAUSS','CARRE', 250 : $ 3,4,2, 251 : $ PGCOUR, 252 : $ IMPR,IRET) 253 : IF (IRET.NE.0) GOTO 9999 254 : NOPG=0 255 : XCOR(1)=X1D3 256 : XCOR(2)=Y1D3 256 : CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 257 : IF (IRET.NE.0) GOTO 9999 258 : SEGDES PGCOUR 259 : MYPGS.LISPG(**)=PGCOUR 260 : * 261 : * Méthode de nom : GPC2-3-4 262 : * Sur un carré : méthode gauss-produit d'ordre 3 à 4 points 263 : * espace de référence de dimension 2 264 : * 265 : * In INIPG : SEGINI PGCOUR 265 : CALL INIPG('GPC2-3-4','GAUSS-PRODUIT','CARRE', 266 : $ 3,4,2, 267 : $ PGCOUR, 268 : $ IMPR,IRET) 269 : IF (IRET.NE.0) GOTO 9999 269 : CALL FIPG('GAC1-3-2',MYPGS,PGPRO1,IMPR,IRET) 270 : IF (IRET.NE.0) GOTO 9999 270 : CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET) 271 : IF (IRET.NE.0) GOTO 9999 272 : SEGDES PGCOUR 273 : MYPGS.LISPG(**)=PGCOUR 274 : * 275 : * Méthode de nom : GAC2-4-6C 276 : * Sur un carré : cubature d'ordre 4 à 6 points 277 : * espace de référence de dimension 2 278 : * 279 : * In INIPG : SEGINI PGCOUR 279 : CALL INIPG('GAC2-4-6C','GAUSS','CARRE', 280 : $ 4,6,2, 281 : $ PGCOUR, 282 : $ IMPR,IRET) 283 : IF (IRET.NE.0) GOTO 9999 284 : NOPG=0 285 : XCOR(1)=X0D4 286 : XCOR(2)=Y0D4 286 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D4,IMPR,IRET) 287 : IF (IRET.NE.0) GOTO 9999 288 : XCOR(1)=X1D4 289 : XCOR(2)=Y1D4 289 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D4,IMPR,IRET) 290 : IF (IRET.NE.0) GOTO 9999 291 : XCOR(1)=X2D4 292 : XCOR(2)=Y2D4 292 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P2D4,IMPR,IRET) 293 : IF (IRET.NE.0) GOTO 9999 294 : XCOR(1)=X3D4 295 : XCOR(2)=Y3D4 295 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P3D4,IMPR,IRET) 296 : IF (IRET.NE.0) GOTO 9999 297 : SEGDES PGCOUR 298 : MYPGS.LISPG(**)=PGCOUR 299 : * 300 : * Méthode de nom : GAC2-5-7A 301 : * Sur un carré : cubature d'ordre 5 à 7 points 302 : * espace de référence de dimension 2 303 : * 304 : * In INIPG : SEGINI PGCOUR 304 : CALL INIPG('GAC2-5-7A','GAUSS','CARRE', 305 : $ 5,7,2, 306 : $ PGCOUR, 307 : $ IMPR,IRET) 308 : IF (IRET.NE.0) GOTO 9999 309 : NOPG=0 310 : XCOR(1)=X0D5 311 : XCOR(2)=Y0D5 311 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D5,IMPR,IRET) 312 : IF (IRET.NE.0) GOTO 9999 313 : XCOR(1)=X1D5 314 : XCOR(2)=Y1D5 314 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET) 315 : IF (IRET.NE.0) GOTO 9999 316 : XCOR(1)=X2D5 317 : XCOR(2)=Y2D5 317 : CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET) 318 : IF (IRET.NE.0) GOTO 9999 319 : SEGDES PGCOUR 320 : MYPGS.LISPG(**)=PGCOUR 321 : * 322 : * Méthode de nom : GPC2-5-9 323 : * Sur un carré : méthode gauss-produit d'ordre 5 à 9 points 324 : * espace de référence de dimension 2 325 : * 326 : * In INIPG : SEGINI PGCOUR 326 : CALL INIPG('GPC2-5-9','GAUSS-PRODUIT','CARRE', 327 : $ 5,9,2, 328 : $ PGCOUR, 329 : $ IMPR,IRET) 330 : IF (IRET.NE.0) GOTO 9999 330 : CALL FIPG('GAC1-5-3',MYPGS,PGPRO1,IMPR,IRET) 331 : IF (IRET.NE.0) GOTO 9999 331 : CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET) 332 : IF (IRET.NE.0) GOTO 9999 333 : SEGDES PGCOUR 334 : MYPGS.LISPG(**)=PGCOUR 335 : * 336 : * Méthode de nom : GAC2-6-10C 337 : * Sur un carré : cubature d'ordre 6 à 10 points 338 : * espace de référence de dimension 2 339 : * 340 : * In INIPG : SEGINI PGCOUR 340 : CALL INIPG('GAC2-6-10C','GAUSS','CARRE', 341 : $ 6,10,2, 342 : $ PGCOUR, 343 : $ IMPR,IRET) 344 : IF (IRET.NE.0) GOTO 9999 345 : NOPG=0 346 : XCOR(1)=X1D6 347 : XCOR(2)=Y1D6 347 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D6,IMPR,IRET) 348 : IF (IRET.NE.0) GOTO 9999 349 : XCOR(1)=X2D6 350 : XCOR(2)=Y2D6 350 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D6,IMPR,IRET) 351 : IF (IRET.NE.0) GOTO 9999 352 : XCOR(1)=X3D6 353 : XCOR(2)=Y3D6 353 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P3D6,IMPR,IRET) 354 : IF (IRET.NE.0) GOTO 9999 355 : XCOR(1)=X4D6 356 : XCOR(2)=Y4D6 356 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P4D6,IMPR,IRET) 357 : IF (IRET.NE.0) GOTO 9999 358 : XCOR(1)=X5D6 359 : XCOR(2)=Y5D6 359 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P5D6,IMPR,IRET) 360 : IF (IRET.NE.0) GOTO 9999 361 : XCOR(1)=X6D6 362 : XCOR(2)=Y6D6 362 : CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P6D6,IMPR,IRET) 363 : IF (IRET.NE.0) GOTO 9999 364 : SEGDES PGCOUR 365 : MYPGS.LISPG(**)=PGCOUR 366 : * 367 : * Méthode de nom : GAC2-7-12A 368 : * Sur un carré : cubature d'ordre 7 à 12 points 369 : * espace de référence de dimension 2 370 : * 371 : * In INIPG : SEGINI PGCOUR 371 : CALL INIPG('GAC2-7-12A','GAUSS','CARRE', 372 : $ 7,12,2, 373 : $ PGCOUR, 374 : $ IMPR,IRET) 375 : IF (IRET.NE.0) GOTO 9999 376 : NOPG=0 377 : XCOR(1)=X1D7 378 : XCOR(2)=Y1D7 378 : CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET) 379 : IF (IRET.NE.0) GOTO 9999 380 : XCOR(1)=X2D7 381 : XCOR(2)=Y2D7 381 : CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET) 382 : IF (IRET.NE.0) GOTO 9999 383 : XCOR(1)=X3D7 384 : XCOR(2)=Y3D7 384 : CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P3D7,IMPR,IRET) 385 : IF (IRET.NE.0) GOTO 9999 386 : SEGDES PGCOUR 387 : MYPGS.LISPG(**)=PGCOUR 388 : * 389 : * Méthode de nom : GPC2-7-16 390 : * Sur un carré : méthode gauss-produit d'ordre 7 à 16 points 391 : * espace de référence de dimension 2 392 : * 393 : * In INIPG : SEGINI PGCOUR 393 : CALL INIPG('GPC2-7-16','GAUSS-PRODUIT','CARRE', 394 : $ 7,16,2, 395 : $ PGCOUR, 396 : $ IMPR,IRET) 397 : IF (IRET.NE.0) GOTO 9999 397 : CALL FIPG('GAC1-7-4',MYPGS,PGPRO1,IMPR,IRET) 398 : IF (IRET.NE.0) GOTO 9999 398 : CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET) 399 : IF (IRET.NE.0) GOTO 9999 400 : SEGDES PGCOUR 401 : MYPGS.LISPG(**)=PGCOUR 402 : * 403 : * Normal termination 404 : * 405 : IRET=0 406 : RETURN 407 : * 408 : * Format handling 409 : * 410 : * 411 : * Error handling 412 : * 413 : 9999 CONTINUE 414 : IRET=1 415 : WRITE(IOIMP,*) 'An error was detected in subroutine ingaqu' 416 : RETURN 417 : * 418 : * End of subroutine INGAQU 419 : * 420 : END 421 : 422 : 423 : 424 :
© Cast3M 2003 - All rights reserved.
Disclaimer