1 : C INGASE SOURCE GOUNAND 05/12/21 21:32:14 5281 2 : SUBROUTINE INGASE(MYPGS,IMPR,IRET) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : INGASE 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 Gauss 10 : C à une dimension (ordre 1 à 11). 11 : C 12 : C REFERENCES : Numerical recipes (sous-programme gauleg modifié) 13 : C on a recalculé les poids et points de Gauss en REAL*16 14 : C donc avec environ 32 (plutôt 31) chiffres significatifs 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, GCCESY 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, 19/10/99, version initiale 28 : C HISTORIQUE : v1, 19/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 IMPR,IRET 54 : * 55 : INTEGER DIMSRF 56 : PARAMETER(DIMSRF=1) 57 : REAL*8 XCOR(DIMSRF) 58 : REAL*8 POIDS 59 : * 60 : * Pour les méthodes sur un segment, tous les générateurs sont 61 : * à symétrie centrales 62 : * 63 : * Générateurs pour la méthode de Gauss de degré 1 : GAC1-1-1 : 64 : * - l'origine 65 : REAL*8 X0D1,P0D1 66 : PARAMETER (X0D1=0.D0) 67 : PARAMETER (P0D1=2.D0) 68 : * 69 : * Générateurs pour la méthode de Gauss de degré 3 : GAC1-3-2 : 70 : * - le générateur symétrique (1./\sqrt{3.});(1.) 71 : REAL*8 X1D3,P1D3 72 : PARAMETER (X1D3=0.57735026918962576450914878050195D0) 73 : PARAMETER (P1D3=1.D0) 74 : * 75 : * Générateurs pour la méthode de Gauss de degré 5 : GAC1-5-3 : 76 : * - l'origine 77 : REAL*8 X0D5,P0D5 78 : PARAMETER (X0D5=0.D0) 79 : PARAMETER (P0D5=8.D0/9.D0) 80 : * - le générateur symétrique (\sqrt{3./5.});(5./9.) 81 : REAL*8 X1D5,P1D5 82 : PARAMETER (X1D5=0.77459666924148337703585307995648D0) 83 : PARAMETER (P1D5=5.D0/9.D0) 84 : * 85 : * Générateurs pour la méthode de Gauss de degré 7 : GAC1-7-4 : 86 : * - les générateurs symétriques : 87 : * (\sqrt{\frac{3.-2\sqrt{6./5.}}{7.}});(1./2.+1./(6\sqrt{6./5.})) 88 : * (\sqrt{\frac{3.+2\sqrt{6./5.}}{7.}});(1./2.-1./(6\sqrt{6./5.})) 89 : REAL*8 X1D7,P1D7,X2D7,P2D7 90 : PARAMETER (X1D7=0.33998104358485626480266575910324D0) 91 : PARAMETER (P1D7=0.65214515486254614262693605077800D0) 92 : PARAMETER (X2D7=0.86113631159405257522394648889281D0) 93 : PARAMETER (P2D7=0.34785484513745385737306394922200D0) 94 : * 95 : * Générateurs pour la méthode de Gauss de degré 9 : GAC1-9-5 : 96 : * - l'origine 97 : REAL*8 X0D9,P0D9 98 : PARAMETER (X0D9=0.D0) 99 : PARAMETER (P0D9=128.D0/225.D0) 100 : * - les générateurs symétriques 101 : * ((1./3.)\sqrt{5.-4.\sqrt{5./14.}});(161./450.+13./(180.\sqrt{5./14.})) 102 : * ((1./3.)\sqrt{5.+4.\sqrt{5./14.}});(161./450.-13./(180.\sqrt{5./14.})) 103 : REAL*8 X1D9,P1D9,X2D9,P2D9 104 : PARAMETER (X1D9=0.53846931010568309103631442070021D0) 105 : PARAMETER (P1D9=0.47862867049936646804129151483564D0) 106 : PARAMETER (X2D9=0.90617984593866399279762687829939D0) 107 : PARAMETER (P2D9=0.23692688505618908751426404071992D0) 108 : * 109 : * Générateurs pour la méthode de Gauss de degré 11 : GAC1-11-6 : 110 : * - les générateurs symétriques : 111 : REAL*8 X1D11,P1D11,X2D11,P2D11,X3D11,P3D11 112 : PARAMETER (X1D11=0.23861918608319690863050172168071D0) 113 : PARAMETER (P1D11=0.46791393457269104738987034398956D0) 114 : PARAMETER (X2D11=0.66120938646626451366139959501990D0) 115 : PARAMETER (P2D11=0.36076157304813860756983351383773D0) 116 : PARAMETER (X3D11=0.93246951420315202781230155449399D0) 117 : PARAMETER (P3D11=0.17132449237917034504029614217274D0) 118 : * 119 : INTEGER NOPG 120 : * 121 : * Executable statements 122 : * 123 : IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingase' 124 : * 125 : * Méthode de nom : NCC1-1-2 126 : * Sur un segment : Cubature d'ordre 1 à 2 points 127 : * espace de référence de dimension 1 128 : * 129 : * In INIPG : SEGINI PGCOUR 129 : CALL INIPG('NCC1-1-2','NEWTON-COTES','SEGMENT', 130 : $ 1,2,1, 131 : $ PGCOUR, 132 : $ IMPR,IRET) 133 : IF (IRET.NE.0) GOTO 9999 134 : NOPG=0 135 : XCOR(1)=-1.D0 135 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 136 : XCOR(1)=1.D0 136 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 137 : IF (IRET.NE.0) GOTO 9999 138 : SEGDES PGCOUR 139 : MYPGS.LISPG(**)=PGCOUR 140 : * 141 : * Méthode de nom : NCC1-3-3 142 : * Sur un segment : Cubature d'ordre 3 à 3 points 143 : * espace de référence de dimension 1 144 : * 145 : * In INIPG : SEGINI PGCOUR 145 : CALL INIPG('NCC1-3-3','NEWTON-COTES','SEGMENT', 146 : $ 3,3,1, 147 : $ PGCOUR, 148 : $ IMPR,IRET) 149 : IF (IRET.NE.0) GOTO 9999 150 : NOPG=0 151 : XCOR(1)=0.D0 152 : POIDS=4.D0/3.D0 152 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET) 153 : XCOR(1)=1.D0 154 : POIDS=1.D0/3.D0 154 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET) 155 : IF (IRET.NE.0) GOTO 9999 156 : SEGDES PGCOUR 157 : MYPGS.LISPG(**)=PGCOUR 158 : * 159 : * Méthode de nom : GAC1-1-1 160 : * Sur un segment : méthode de Gauss d'ordre 1 à 1 point 161 : * espace de référence de dimension 1 162 : * 163 : * In INIPG : SEGINI PGCOUR 163 : CALL INIPG('GAC1-1-1','GAUSS','SEGMENT', 164 : $ 1,1,1, 165 : $ PGCOUR, 166 : $ IMPR,IRET) 167 : IF (IRET.NE.0) GOTO 9999 168 : NOPG=0 169 : XCOR(1)=X0D1 169 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D1,IMPR,IRET) 170 : IF (IRET.NE.0) GOTO 9999 171 : SEGDES PGCOUR 172 : MYPGS.LISPG(**)=PGCOUR 173 : * 174 : * Méthode de nom : GAC1-3-2 175 : * Sur un segment : méthode de Gauss d'ordre 3 à 2 points 176 : * espace de référence de dimension 1 177 : * 178 : * In INIPG : SEGINI PGCOUR 178 : CALL INIPG('GAC1-3-2','GAUSS','SEGMENT', 179 : $ 3,2,1, 180 : $ PGCOUR, 181 : $ IMPR,IRET) 182 : IF (IRET.NE.0) GOTO 9999 183 : NOPG=0 184 : XCOR(1)=X1D3 184 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) 185 : IF (IRET.NE.0) GOTO 9999 186 : SEGDES PGCOUR 187 : MYPGS.LISPG(**)=PGCOUR 188 : * 189 : * Méthode de nom : GAC1-5-3 190 : * Sur un segment : méthode de Gauss d'ordre 5 à 3 points 191 : * espace de référence de dimension 1 192 : * 193 : * In INIPG : SEGINI PGCOUR 193 : CALL INIPG('GAC1-5-3','GAUSS','SEGMENT', 194 : $ 5,3,1, 195 : $ PGCOUR, 196 : $ IMPR,IRET) 197 : IF (IRET.NE.0) GOTO 9999 198 : NOPG=0 199 : XCOR(1)=X0D5 199 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D5,IMPR,IRET) 200 : IF (IRET.NE.0) GOTO 9999 201 : XCOR(1)=X1D5 201 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET) 202 : IF (IRET.NE.0) GOTO 9999 203 : SEGDES PGCOUR 204 : MYPGS.LISPG(**)=PGCOUR 205 : * 206 : * Méthode de nom : GAC1-7-4 207 : * Sur un segment : méthode de Gauss d'ordre 7 à 4 points 208 : * espace de référence de dimension 1 209 : * 210 : * In INIPG : SEGINI PGCOUR 210 : CALL INIPG('GAC1-7-4','GAUSS','SEGMENT', 211 : $ 7,4,1, 212 : $ PGCOUR, 213 : $ IMPR,IRET) 214 : IF (IRET.NE.0) GOTO 9999 215 : NOPG=0 216 : XCOR(1)=X1D7 216 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET) 217 : IF (IRET.NE.0) GOTO 9999 218 : XCOR(1)=X2D7 218 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET) 219 : IF (IRET.NE.0) GOTO 9999 220 : SEGDES PGCOUR 221 : MYPGS.LISPG(**)=PGCOUR 222 : * 223 : * Méthode de nom : GAC1-9-5 224 : * Sur un segment : méthode de Gauss d'ordre 9 à 5 points 225 : * espace de référence de dimension 1 226 : * 227 : * In INIPG : SEGINI PGCOUR 227 : CALL INIPG('GAC1-9-5','GAUSS','SEGMENT', 228 : $ 9,5,1, 229 : $ PGCOUR, 230 : $ IMPR,IRET) 231 : IF (IRET.NE.0) GOTO 9999 232 : NOPG=0 233 : XCOR(1)=X0D9 233 : CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D9,IMPR,IRET) 234 : IF (IRET.NE.0) GOTO 9999 235 : XCOR(1)=X1D9 235 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D9,IMPR,IRET) 236 : IF (IRET.NE.0) GOTO 9999 237 : XCOR(1)=X2D9 237 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D9,IMPR,IRET) 238 : IF (IRET.NE.0) GOTO 9999 239 : SEGDES PGCOUR 240 : MYPGS.LISPG(**)=PGCOUR 241 : * 242 : * Méthode de nom : GAC1-11-6 243 : * Sur un segment : méthode de Gauss d'ordre 11 à 6 points 244 : * espace de référence de dimension 1 245 : * 246 : * In INIPG : SEGINI PGCOUR 246 : CALL INIPG('GAC1-11-6','GAUSS','SEGMENT', 247 : $ 11,6,1, 248 : $ PGCOUR, 249 : $ IMPR,IRET) 250 : IF (IRET.NE.0) GOTO 9999 251 : NOPG=0 252 : XCOR(1)=X1D11 252 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D11,IMPR,IRET) 253 : IF (IRET.NE.0) GOTO 9999 254 : XCOR(1)=X2D11 254 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D11,IMPR,IRET) 255 : IF (IRET.NE.0) GOTO 9999 256 : XCOR(1)=X3D11 256 : CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P3D11,IMPR,IRET) 257 : IF (IRET.NE.0) GOTO 9999 258 : SEGDES PGCOUR 259 : MYPGS.LISPG(**)=PGCOUR 260 : * 261 : * Normal termination 262 : * 263 : IRET=0 264 : RETURN 265 : * 266 : * Format handling 267 : * 268 : * 269 : * Error handling 270 : * 271 : 9999 CONTINUE 272 : IRET=1 273 : WRITE(IOIMP,*) 'An error was detected in subroutine ingase' 274 : RETURN 275 : * 276 : * End of subroutine INGASE 277 : * 278 : END 279 : 280 : 281 :
© Cast3M 2003 - All rights reserved.
Disclaimer