C INGJ10 SOURCE GOUNAND 21/06/02 21:16:45 11022 SUBROUTINE INGJ10(MYPGS,IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : INGJ10 C PROJET : Noyau linéaire NLIN C DESCRIPTION : Remplit le segment des méthodes d'intégration C avec des méthodes d'intégration numérique de C Gauss-Jacobi (\alpha=1, \beta=0) à une dimension sur C l'intervalle [0,1] (ordre 1 à 11). C C On intègre donc \int_0^1 (1-x) f(x) dx de manière C approchée. C C Ces méthodes sont utilisés pour générer des formules C produits pour les éléments de type cônes : triangles. C C REFERENCES : Numerical recipes (sous-programme gaujac modifié) C on a recalculé les poids et points de Gauss en REAL*16 C donc avec environ 32 (plutôt 31) chiffres significatifs C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : INIPG, GCSINO C APPELE PAR : INPGS C*********************************************************************** C ENTREES : - C ENTREES/SORTIES : MYPGS (actif en *MOD) C SORTIES : - C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 31/05/00, version initiale C HISTORIQUE : v1, 31/05/00, création C HISTORIQUE : C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC TNLIN *-INC SPOGAU POINTEUR MYPGS.POGAUS POINTEUR PGCOUR.POGAU * INTEGER IMPR,IRET * INTEGER DIMSRF PARAMETER(DIMSRF=1) REAL*8 XCOR(DIMSRF) * * Générateurs pour la méthode de Gauss-Jacobi de degré 1 : GJ10-1-1 : * REAL*8 X1D1,P1D1 PARAMETER (X1D1=1.D0/3.D0) PARAMETER (P1D1=0.5D0) * * Générateurs pour la méthode de Gauss-Jacobi de degré 3 : GJ10-3-2 : * REAL*8 X1D3,P1D3,X2D3,P2D3 PARAMETER (X1D3=0.15505102572168219018027159252941D0) PARAMETER (P1D3=0.31804138174397716939436900207515D0) PARAMETER (X2D3=0.64494897427831780981972840747060D0) PARAMETER (P2D3=0.18195861825602283060563099792484D0) * * Générateurs pour la méthode de Gauss-Jacobi de degré 5 : GJ10-5-3 : * REAL*8 X1D5,P1D5,X2D5,P2D5,X3D5,P3D5 PARAMETER (X1D5=0.88587959512703947395546143769455D-1) PARAMETER (P1D5=0.20093191373895963077219813326462D0) PARAMETER (X2D5=0.40946686444073471086492625206883D0) PARAMETER (P2D5=0.22924110635958624669392059455632D0) PARAMETER (X3D5=0.78765946176084705602524188987600D0) PARAMETER (P3D5=0.69826979901454122533881272179078D-1) * * Générateurs pour la méthode de Gauss-Jacobi de degré 7 : GJ10-7-4 : * REAL*8 X1D7,P1D7,X2D7,P2D7,X3D7,P3D7,X4D7,P4D7 PARAMETER (X1D7=0.57104196114517682193121192554116D-1) PARAMETER (P1D7=0.13550691343148811620826417407793D0) PARAMETER (X2D7=0.27684301363812382768004599768562D0) PARAMETER (P2D7=0.20346456801027136079140447593585D0) PARAMETER (X3D7=0.58359043236891682005669766866292D0) PARAMETER (P3D7=0.12984754760823244082645620288963D0) PARAMETER (X4D7=0.86024013565621944784791291887512D0) PARAMETER (P4D7=0.31180970950008082173875147096569D-1) * * Générateurs pour la méthode de Gauss-Jacobi de degré 9 : GJ10-9-5 : * REAL*8 X1D9,P1D9,X2D9,P2D9,X3D9,P3D9,X4D9,P4D9,X5D9,P5D9 PARAMETER (X1D9=0.39809857051468742340806690093331D-1) PARAMETER (P1D9=0.96781590226651679274360971636169D-1) PARAMETER (X2D9=0.19801341787360817253579213679530D0) PARAMETER (P2D9=0.16717463809436956549167562309770D0) PARAMETER (X3D9=0.43797481024738614400501252000523D0) PARAMETER (P3D9=0.14638698708466980869803786935596D0) PARAMETER (X4D9=0.69546427335363609451461482372117D0) PARAMETER (P4D9=0.73908870072616670350633219341705D-1) PARAMETER (X5D9=0.90146491420117357387650110211225D0) PARAMETER (P5D9=0.15747914521692276185292316568490D-1) * * Générateurs pour la méthode de Gauss-Jacobi de degré 11 : GJ10-11-6 : * REAL*8 X1D11,P1D11,X2D11,P2D11,X3D11,P3D11 REAL*8 X4D11,P4D11,X5D11,P5D11,X6D11,P6D11 PARAMETER (X1D11=0.29316427159784891972050276913165D-1) PARAMETER (P1D11=0.72310330725508683655454326124839D-1) PARAMETER (X2D11=0.14807859966848429184997685249598D0) PARAMETER (P2D11=0.13554249723151861684069039663804D0) PARAMETER (X3D11=0.33698469028115429909705297208078D0) PARAMETER (P3D11=0.14079255378819892811907683907092D0) PARAMETER (X4D11=0.55867151877155013208139334180552D0) PARAMETER (P4D11=0.98661150890655264120584510548346D-1) PARAMETER (X5D11=0.76923386203005450091688336011564D0) PARAMETER (P5D11=0.43955165550508975508176624305422D-1) PARAMETER (X6D11=0.92694567131974111485187396581969D0) PARAMETER (P6D11=0.87383018136095317560173033123964D-2) * INTEGER NOPG * * Executable statements * IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingj10' * * Méthode de nom : GJ10-1-1 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 1 à 1 point * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-1-1','GAUSS-JACOBI10','SEGMENT', $ 1,1,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D1 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D1,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Méthode de nom : GJ10-3-2 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 3 à 2 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-3-2','GAUSS-JACOBI10','SEGMENT', $ 3,2,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D3 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X2D3 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D3,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Méthode de nom : GJ10-5-3 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 5 à 3 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-5-3','GAUSS-JACOBI10','SEGMENT', $ 5,3,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D5 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X2D5 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X3D5 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P3D5,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Méthode de nom : GJ10-7-4 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 7 à 4 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-7-4','GAUSS-JACOBI10','SEGMENT', $ 7,4,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D7 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X2D7 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X3D7 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P3D7,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X4D7 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P4D7,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Méthode de nom : GJ10-9-5 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 9 à 5 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-9-5','GAUSS-JACOBI10','SEGMENT', $ 9,5,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D9 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D9,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X2D9 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D9,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X3D9 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P3D9,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X4D9 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P4D9,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X5D9 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P5D9,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Méthode de nom : GJ10-11-6 * Sur un segment : méthode de Gauss-Jacobi \alpha=1 \beta=0 * d'ordre 11 à 6 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ10-11-6','GAUSS-JACOBI10','SEGMENT', $ 11,6,1, $ PGCOUR, $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 NOPG=0 XCOR(1)=X1D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X2D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X3D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P3D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X4D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P4D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X5D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P5D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 XCOR(1)=X6D11 CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P6D11,IMPR,IRET) IF (IRET.NE.0) GOTO 9999 SEGDES PGCOUR MYPGS.LISPG(**)=PGCOUR * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ingj10' RETURN * * End of subroutine INGJ10 * END