1 : C INGATE SOURCE GOUNAND 05/12/21 21:32:21 5281 2 : SUBROUTINE INGATE(MYPGS,IMPR,IRET) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : INGATE 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 tétraèdre (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, GT3FS9, GT3F10, 20 : C FIPG, CPROPG 21 : C APPELE PAR : INPGS 22 : C*********************************************************************** 23 : C ENTREES : - 24 : C ENTREES/SORTIES : MYPGS (actif en *MOD) 25 : C SORTIES : - 26 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 27 : C*********************************************************************** 28 : C VERSION : v1, 22/10/99, version initiale 29 : C HISTORIQUE : v1, 22/10/99, création 30 : C HISTORIQUE : 29/5/00 rajout ordre 6 31 : C HISTORIQUE : 32 : C*********************************************************************** 33 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 34 : C en cas de modification de ce sous-programme afin de faciliter 35 : C la maintenance ! 36 : C*********************************************************************** 37 : -INC CCOPTIO 38 : CBEGININCLUDE SPOGAU 39 : SEGMENT POGAU 40 : CHARACTER*(LNNPG) NOMPG 41 : CHARACTER*(LNTPG) TYPMPG 42 : CHARACTER*(LNFPG) FORLPG 43 : INTEGER NORDPG 44 : REAL*8 XCOPG(NDLPG,NBPG) 45 : REAL*8 XPOPG(NBPG) 46 : ENDSEGMENT 47 : SEGMENT POGAUS 48 : POINTEUR LISPG(0).POGAU 49 : ENDSEGMENT 50 : CENDINCLUDE SPOGAU 51 : POINTEUR MYPGS.POGAUS 52 : POINTEUR PGCOUR.POGAU 53 : POINTEUR PGPRO1.POGAU 54 : POINTEUR PGPRO2.POGAU 55 : * 56 : INTEGER IMPR,IRET 57 : * On travaille en coordonnées barycentriques. 58 : * D'après le Dhatt et Touzot 59 : * L1 = 1 - \ksi - \eta -\zeta 60 : * L2 = \ksi 61 : * L3 = \eta 62 : * L4 = \zeta 63 : INTEGER DIMSRF 64 : PARAMETER(DIMSRF=3) 65 : REAL*8 XCOR(DIMSRF+1) 66 : * 67 : * Générateurs pour la cubature de degré 1 à 1 point : GAT3-1-1 : 68 : * - [ Fully symmetric ] 69 : REAL*8 X1D1,Y1D1,Z1D1,P1D1 70 : PARAMETER (X1D1=1.D0/4.D0) 71 : PARAMETER (Y1D1=1.D0/4.D0) 72 : PARAMETER (Z1D1=1.D0/4.D0) 73 : PARAMETER (P1D1=1.D0/6.D0) 74 : * 75 : * Générateurs pour la cubature de degré 2 à 4 points : GAT3-2-4B : 76 : * - [ Fully symmetric ] 77 : REAL*8 X1D2,Y1D2,Z1D2,P1D2 78 : PARAMETER (X1D2=0.138196601125010515179541316563436D0) 79 : PARAMETER (Y1D2=0.138196601125010515179541316563436D0) 80 : PARAMETER (Z1D2=0.138196601125010515179541316563436D0) 81 : PARAMETER (P1D2=1.D0/24.D0) 82 : * 83 : * Générateurs pour la cubature de degré 3 à 8 points : GAT3-3-8B : 84 : * - [ Fully symmetric ] 85 : REAL*8 X1D3,Y1D3,Z1D3,P1D3 86 : PARAMETER (X1D3=1.D0) 87 : PARAMETER (Y1D3=0.D0) 88 : PARAMETER (Z1D3=0.D0) 89 : PARAMETER (P1D3=4.16666666666666666666666666666666D-3) 90 : * - [ Fully symmetric ] 91 : REAL*8 X2D3,Y2D3,Z2D3,P2D3 92 : PARAMETER (X2D3=1.D0/3.D0) 93 : PARAMETER (Y2D3=1.D0/3.D0) 94 : PARAMETER (Z2D3=1.D0/3.D0) 95 : PARAMETER (P2D3=0.0375D0) 96 : * 97 : * Générateurs pour la cubature de degré 5 à 14 points : GAT3-5-14A : 98 : * - [ Fully symmetric ] 99 : REAL*8 X1D5,Y1D5,Z1D5,P1D5 100 : PARAMETER (X1D5=0.0927352503108912264023239137370306D0) 101 : PARAMETER (Y1D5=0.0927352503108912264023239137370306D0) 102 : PARAMETER (Z1D5=0.0927352503108912264023239137370306D0) 103 : PARAMETER (P1D5=0.0122488405193936582572850342477212D0) 104 : * - [ Fully symmetric ] 105 : REAL*8 X2D5,Y2D5,Z2D5,P2D5 106 : PARAMETER (X2D5=0.310885919263300609797345733763457D0) 107 : PARAMETER (Y2D5=0.310885919263300609797345733763457D0) 108 : PARAMETER (Z2D5=0.310885919263300609797345733763457D0) 109 : PARAMETER (P2D5=0.0187813209530026417998642753888810D0) 110 : * - [ Fully symmetric ] 111 : REAL*8 X3D5,Y3D5,Z3D5,P3D5 112 : PARAMETER (X3D5=0.454496295874350350508119473720660D0) 113 : PARAMETER (Y3D5=0.454496295874350350508119473720660D0) 114 : PARAMETER (Z3D5=0.0455037041256496494918805262793394D0) 115 : PARAMETER (P3D5=7.09100346284691107301157135337624D-3) 116 : * 117 : * Générateurs pour la cubature de degré 6 à 24 points : GAT3-6-24 : 118 : * - [ Fully symmetric ] 119 : REAL*8 X1D6,Y1D6,Z1D6,P1D6 120 : PARAMETER (X1D6=0.214602871259152029288839219386284D0) 121 : PARAMETER (Y1D6=0.214602871259152029288839219386284D0) 122 : PARAMETER (Z1D6=0.214602871259152029288839219386284D0) 123 : PARAMETER (P1D6=6.65379170969458201661510459291332D-3) 124 : * - [ Fully symmetric ] 125 : REAL*8 X2D6,Y2D6,Z2D6,P2D6 126 : PARAMETER (X2D6=0.0406739585346113531155794489564100D0) 127 : PARAMETER (Y2D6=0.0406739585346113531155794489564100D0) 128 : PARAMETER (Z2D6=0.0406739585346113531155794489564100D0) 129 : PARAMETER (P2D6=1.67953517588677382466887290765614D-3) 130 : * - [ Fully symmetric ] 131 : REAL*8 X3D6,Y3D6,Z3D6,P3D6 132 : PARAMETER (X3D6=0.322337890142275510343994470762492D0) 133 : PARAMETER (Y3D6=0.322337890142275510343994470762492D0) 134 : PARAMETER (Z3D6=0.322337890142275510343994470762492D0) 135 : PARAMETER (P3D6=9.22619692394245368252554630895433D-3) 136 : * - [ Fully symmetric ] 137 : REAL*8 X4D6,Y4D6,Z4D6,P4D6 138 : PARAMETER (X4D6=0.0636610018750175252992355276057269D0) 139 : PARAMETER (Y4D6=0.0636610018750175252992355276057269D0) 140 : PARAMETER (Z4D6=0.269672331458315808034097805727606D0) 141 : PARAMETER (P4D6=8.03571428571428571428571428571428D-3) 142 : * 143 : INTEGER NOPG 144 : * 145 : * Executable statements 146 : * 147 : IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingate' 148 : * 149 : * Méthode de nom : NCT3-1-4 150 : * Sur un tétraèdre : cubature d'ordre 1 à 4 points 151 : * espace de référence de dimension 3 152 : * 153 : * In INIPG : SEGINI PGCOUR 153 : CALL INIPG('NCT3-1-4','NEWTON-COTES','TETRAEDRE', 154 : $ 1,4,3, 155 : $ PGCOUR, 156 : $ IMPR,IRET) 157 : IF (IRET.NE.0) GOTO 9999 158 : NOPG=0 159 : XCOR(1)=1.D0 160 : XCOR(2)=0.D0 161 : XCOR(3)=0.D0 162 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 162 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET) 163 : IF (IRET.NE.0) GOTO 9999 164 : SEGDES PGCOUR 165 : MYPGS.LISPG(**)=PGCOUR 166 : * 167 : * Méthode de nom : GAT3-1-1 168 : * Sur un tétraèdre : cubature d'ordre 1 à 1 point 169 : * espace de référence de dimension 3 170 : * 171 : * In INIPG : SEGINI PGCOUR 171 : CALL INIPG('GAT3-1-1','GAUSS','TETRAEDRE', 172 : $ 1,1,3, 173 : $ PGCOUR, 174 : $ IMPR,IRET) 175 : IF (IRET.NE.0) GOTO 9999 176 : NOPG=0 177 : XCOR(1)=X1D1 178 : XCOR(2)=Y1D1 179 : XCOR(3)=Z1D1 180 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 180 : CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D1,IMPR,IRET) 181 : IF (IRET.NE.0) GOTO 9999 182 : SEGDES PGCOUR 183 : MYPGS.LISPG(**)=PGCOUR 184 : * 185 : * Méthode de nom : GAT3-2-4B 186 : * Sur un tétraèdre : cubature d'ordre 2 à 4 points 187 : * espace de référence de dimension 3 188 : * 189 : * In INIPG : SEGINI PGCOUR 189 : CALL INIPG('GAT3-2-4B','GAUSS','TETRAEDRE', 190 : $ 2,4,3, 191 : $ PGCOUR, 192 : $ IMPR,IRET) 193 : IF (IRET.NE.0) GOTO 9999 194 : NOPG=0 195 : XCOR(1)=X1D2 196 : XCOR(2)=Y1D2 197 : XCOR(3)=Z1D2 198 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 198 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET) 199 : IF (IRET.NE.0) GOTO 9999 200 : SEGDES PGCOUR 201 : MYPGS.LISPG(**)=PGCOUR 202 : * 203 : * Méthode de nom : GAT3-3-8B 204 : * Sur un tétraèdre : cubature d'ordre 3 à 8 points 205 : * espace de référence de dimension 3 206 : * 207 : * In INIPG : SEGINI PGCOUR 207 : CALL INIPG('GAT3-3-8B','GAUSS','TETRAEDRE', 208 : $ 3,8,3, 209 : $ PGCOUR, 210 : $ IMPR,IRET) 211 : IF (IRET.NE.0) GOTO 9999 212 : NOPG=0 213 : XCOR(1)=X1D3 214 : XCOR(2)=Y1D3 215 : XCOR(3)=Z1D3 216 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 216 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 217 : IF (IRET.NE.0) GOTO 9999 218 : XCOR(1)=X2D3 219 : XCOR(2)=Y2D3 220 : XCOR(3)=Z2D3 221 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 221 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D3,IMPR,IRET) 222 : IF (IRET.NE.0) GOTO 9999 223 : SEGDES PGCOUR 224 : MYPGS.LISPG(**)=PGCOUR 225 : * 226 : * Méthode de nom : GAT3-5-14 227 : * Sur un tétraèdre : cubature d'ordre 5 à 14 points 228 : * espace de référence de dimension 3 229 : * 230 : * In INIPG : SEGINI PGCOUR 230 : CALL INIPG('GAT3-5-14','GAUSS','TETRAEDRE', 231 : $ 5,14,3, 232 : $ PGCOUR, 233 : $ IMPR,IRET) 234 : IF (IRET.NE.0) GOTO 9999 235 : NOPG=0 236 : XCOR(1)=X1D5 237 : XCOR(2)=Y1D5 238 : XCOR(3)=Z1D5 239 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 239 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET) 240 : IF (IRET.NE.0) GOTO 9999 241 : XCOR(1)=X2D5 242 : XCOR(2)=Y2D5 243 : XCOR(3)=Z2D5 244 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 244 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET) 245 : IF (IRET.NE.0) GOTO 9999 246 : XCOR(1)=X3D5 247 : XCOR(2)=Y3D5 248 : XCOR(3)=Z3D5 249 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 249 : CALL GT3FS9(PGCOUR,NOPG,DIMSRF,XCOR,P3D5,IMPR,IRET) 250 : IF (IRET.NE.0) GOTO 9999 251 : SEGDES PGCOUR 252 : MYPGS.LISPG(**)=PGCOUR 253 : * 254 : * Méthode de nom : GPT3-5-21 255 : * Sur un tétraèdre : méthode gauss-produit conique d'ordre 5 à 21 points 256 : * espace de référence de dimension 3 257 : * 258 : * In INIPG : SEGINI PGCOUR 258 : CALL INIPG('GPT3-5-21','GAUSS-PRODUIT','TETRAEDRE', 259 : $ 5,21,3, 260 : $ PGCOUR, 261 : $ IMPR,IRET) 262 : IF (IRET.NE.0) GOTO 9999 262 : CALL FIPG('GAT2-5-7',MYPGS,PGPRO1,IMPR,IRET) 263 : IF (IRET.NE.0) GOTO 9999 263 : CALL FIPG('GJ20-5-3',MYPGS,PGPRO2,IMPR,IRET) 264 : IF (IRET.NE.0) GOTO 9999 264 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 265 : IF (IRET.NE.0) GOTO 9999 266 : SEGDES PGCOUR 267 : MYPGS.LISPG(**)=PGCOUR 268 : * 269 : * Méthode de nom : GPT3-5-27 270 : * Sur un tétraèdre : méthode gauss-produit conique d'ordre 5 à 27 points 271 : * espace de référence de dimension 3 272 : * 273 : * In INIPG : SEGINI PGCOUR 273 : CALL INIPG('GPT3-5-27','GAUSS-PRODUIT','TETRAEDRE', 274 : $ 5,27,3, 275 : $ PGCOUR, 276 : $ IMPR,IRET) 277 : IF (IRET.NE.0) GOTO 9999 277 : CALL FIPG('GPT2-5-9',MYPGS,PGPRO1,IMPR,IRET) 278 : IF (IRET.NE.0) GOTO 9999 278 : CALL FIPG('GJ20-5-3',MYPGS,PGPRO2,IMPR,IRET) 279 : IF (IRET.NE.0) GOTO 9999 279 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 280 : IF (IRET.NE.0) GOTO 9999 281 : SEGDES PGCOUR 282 : MYPGS.LISPG(**)=PGCOUR 283 : * 284 : * Méthode de nom : GAT3-6-24 285 : * Sur un tétraèdre : cubature d'ordre 6 à 24 points 286 : * espace de référence de dimension 3 287 : * 288 : * In INIPG : SEGINI PGCOUR 288 : CALL INIPG('GAT3-6-24','GAUSS','TETRAEDRE', 289 : $ 6,24,3, 290 : $ PGCOUR, 291 : $ IMPR,IRET) 292 : IF (IRET.NE.0) GOTO 9999 293 : NOPG=0 294 : XCOR(1)=X1D6 295 : XCOR(2)=Y1D6 296 : XCOR(3)=Z1D6 297 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 297 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D6,IMPR,IRET) 298 : IF (IRET.NE.0) GOTO 9999 299 : XCOR(1)=X2D6 300 : XCOR(2)=Y2D6 301 : XCOR(3)=Z2D6 302 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 302 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D6,IMPR,IRET) 303 : IF (IRET.NE.0) GOTO 9999 304 : XCOR(1)=X3D6 305 : XCOR(2)=Y3D6 306 : XCOR(3)=Z3D6 307 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 307 : CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P3D6,IMPR,IRET) 308 : IF (IRET.NE.0) GOTO 9999 309 : XCOR(1)=X4D6 310 : XCOR(2)=Y4D6 311 : XCOR(3)=Z4D6 312 : XCOR(4)=1-XCOR(1)-XCOR(2)-XCOR(3) 312 : CALL GT3F10(PGCOUR,NOPG,DIMSRF,XCOR,P4D6,IMPR,IRET) 313 : IF (IRET.NE.0) GOTO 9999 314 : SEGDES PGCOUR 315 : MYPGS.LISPG(**)=PGCOUR 316 : * 317 : * Méthode de nom : GPT3-7-64 318 : * Sur un tétraèdre : méthode gauss-produit conique d'ordre 7 à 64 points 319 : * espace de référence de dimension 3 320 : * 321 : * In INIPG : SEGINI PGCOUR 321 : CALL INIPG('GPT3-7-64','GAUSS-PRODUIT','TETRAEDRE', 322 : $ 7,64,3, 323 : $ PGCOUR, 324 : $ IMPR,IRET) 325 : IF (IRET.NE.0) GOTO 9999 325 : CALL FIPG('GPT2-7-16',MYPGS,PGPRO1,IMPR,IRET) 326 : IF (IRET.NE.0) GOTO 9999 326 : CALL FIPG('GJ20-7-4',MYPGS,PGPRO2,IMPR,IRET) 327 : IF (IRET.NE.0) GOTO 9999 327 : CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET) 328 : IF (IRET.NE.0) GOTO 9999 329 : SEGDES PGCOUR 330 : MYPGS.LISPG(**)=PGCOUR 331 : * 332 : * Normal termination 333 : * 334 : IRET=0 335 : RETURN 336 : * 337 : * Format handling 338 : * 339 : * 340 : * Error handling 341 : * 342 : 9999 CONTINUE 343 : IRET=1 344 : WRITE(IOIMP,*) 'An error was detected in subroutine ingate' 345 : RETURN 346 : * 347 : * End of subroutine INGATE 348 : * 349 : END 350 : 351 : 352 :
© Cast3M 2003 - All rights reserved.
Disclaimer