C INGJ20 SOURCE GOUNAND 21/06/02 21:16:46 11022 SUBROUTINE INGJ20(MYPGS,IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : INGJ20 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=2, \beta=0) à une dimension sur C l'intervalle [0,1] (ordre 1 à 11). C C On intègre donc \int_0^1 (1-x)^2 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 : tétraèdres, C pyramides. 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 : GJ20-1-1 : * REAL*8 X1D1,P1D1 PARAMETER (X1D1=0.25D0) PARAMETER (P1D1=1.D0/3.D0) * * Générateurs pour la méthode de Gauss-Jacobi de degré 3 : GJ20-3-2 : * REAL*8 X1D3,P1D3,X2D3,P2D3 PARAMETER (X1D3=0.12251482265544137786674043037115D0) PARAMETER (P1D3=0.23254745125350790274997694884235D0) PARAMETER (X2D3=0.54415184401122528879992623629551D0) PARAMETER (P2D3=0.10078588207982543058335638449098D0) * * Générateurs pour la méthode de Gauss-Jacobi de degré 5 : GJ20-5-3 : * REAL*8 X1D5,P1D5,X2D5,P2D5,X3D5,P3D5 PARAMETER (X1D5=0.72994024073149732155837979012005D-1) PARAMETER (P1D5=0.15713636106488661331834482221327D0) PARAMETER (X2D5=0.34700376603835188472176354340395D0) PARAMETER (P2D5=0.14624626925986602200351202036424D0) PARAMETER (X3D5=0.70500220988849838312239847758405D0) PARAMETER (P3D5=0.29950703008580698011476490755826D-1) * * Générateurs pour la méthode de Gauss-Jacobi de degré 7 : GJ20-7-4 : * REAL*8 X1D7,P1D7,X2D7,P2D7,X3D7,P3D7,X4D7,P4D7 PARAMETER (X1D7=0.48500549446997329297067257098986D-1) PARAMETER (P1D7=0.11088841561127798368323131746895D0) PARAMETER (X2D7=0.23860073755186230505898141272470D0) PARAMETER (P2D7=0.14345878979921420904832801427594D0) PARAMETER (X3D7=0.51704729510436750234057336938307D0) PARAMETER (P3D7=0.68633887172923075317376345041811D-1) PARAMETER (X4D7=0.79585141789677286330337796079324D0) PARAMETER (P4D7=0.10352240749918065284397656546638D-1) * * Générateurs pour la méthode de Gauss-Jacobi de degré 9 : GJ20-9-5 : * REAL*8 X1D9,P1D9,X2D9,P2D9,X3D9,P3D9,X4D9,P4D9,X5D9,P5D9 PARAMETER (X1D9=0.34578939918215091524457428631527D-1) PARAMETER (P1D9=0.81764784285770917904880732922335D-1) PARAMETER (X2D9=0.17348032077169572310459241798618D0) PARAMETER (P2D9=0.12619896189991148802883293516466D0) PARAMETER (X3D9=0.38988638706551932824089541038500D0) PARAMETER (P3D9=0.89200161221590000186254493070382D-1) PARAMETER (X4D9=0.63433347263088677234716388892062D0) PARAMETER (P4D9=0.32055600722961919254748930556631D-1) PARAMETER (X5D9=0.85105421294701641811622418741001D0) PARAMETER (P5D9=0.41138252030990079586162416192984D-2) * * Générateurs pour la méthode de Gauss-Jacobi de degré 11 : GJ20-11-6 : * REAL*8 X1D11,P1D11,X2D11,P2D11,X3D11,P3D11 REAL*8 X4D11,P4D11,X5D11,P5D11,X6D11,P6D11 PARAMETER (X1D11=0.25904555093667192754643606997233D-1) PARAMETER (P1D11=0.62538702726580937878526556468326D-1) PARAMETER (X2D11=0.13156394165798513398691085074097D0) PARAMETER (P2D11=0.10737649973678063260575568234795D0) PARAMETER (X3D11=0.30243691802289123274990557791855D0) PARAMETER (P3D11=0.94577186748541203568292051720050D-1) PARAMETER (X4D11=0.50903641316475208401103990516772D0) PARAMETER (P4D11=0.51289571129616210220129325076918D-1) PARAMETER (X5D11=0.71568112731171391876766262459361D0) PARAMETER (P5D11=0.15720297184945051327851262130019D-1) PARAMETER (X6D11=0.88680561617756186630126600601049D0) PARAMETER (P6D11=0.18310758068692977327784555900563D-2) * INTEGER NOPG * * Executable statements * IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingj20' * * Méthode de nom : GJ20-1-1 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 1 à 1 point * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-1-1','GAUSS-JACOBI20','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 : GJ20-3-2 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 3 à 2 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-3-2','GAUSS-JACOBI20','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 : GJ20-5-3 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 5 à 3 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-5-3','GAUSS-JACOBI20','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 : GJ20-7-4 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 7 à 4 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-7-4','GAUSS-JACOBI20','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 : GJ20-9-5 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 9 à 5 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-9-5','GAUSS-JACOBI20','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 : GJ20-11-6 * Sur un segment : méthode de Gauss-Jacobi \alpha=2 \beta=0 * d'ordre 11 à 6 points * espace de référence de dimension 1 * * In INIPG : SEGINI PGCOUR CALL INIPG('GJ20-11-6','GAUSS-JACOBI20','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 ingj20' RETURN * * End of subroutine INGJ20 * END