Télécharger la plaquette Cast3M


fr en

Télécharger ingatr.eso

Retour à la liste

Numérotation des lignes : 

ingatr

  1. C INGATR    SOURCE    GOUNAND   21/06/02    21:16:44     11022          
  2.       SUBROUTINE INGATR(MYPGS,IMPR,IRET)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : INGATR
  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 triangle (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, FIPG, CPROPG
  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.  
  37. -INC PPARAM
  38. -INC CCOPTIO
  39. -INC TNLIN      
  40. *-INC SPOGAU
  41.       POINTEUR MYPGS.POGAUS
  42.       POINTEUR PGCOUR.POGAU
  43.       POINTEUR PGPRO1.POGAU
  44.       POINTEUR PGPRO2.POGAU
  45. *
  46.       INTEGER IMPR,IRET
  47. *   On travaille en coordonnées barycentriques.
  48. *   D'après le Dhatt et Touzot
  49. *   L1 = 1 - \ksi - \eta
  50. *   L2 = \ksi
  51. *   L3 = \eta
  52.       INTEGER DIMSRF
  53.       PARAMETER(DIMSRF=2)
  54.       REAL*8 XCOR(DIMSRF+1)
  55.       REAL*8 POIDS
  56. *
  57. *  Générateurs pour la cubature de degré 1 à 1 point  : GAT2-1-1 :
  58. *   - [ Fully symmetric ]
  59.       REAL*8 X1D1,Y1D1,P1D1
  60.       PARAMETER (X1D1=1.D0/3.D0)
  61.       PARAMETER (Y1D1=1.D0/3.D0)
  62.       PARAMETER (P1D1=0.5D0)
  63. *
  64. *  Générateurs pour la cubature de degré 2 à 3 points : GAT2-2-3A :
  65. *   - [ Fully symmetric ]
  66.       REAL*8 X1D2,Y1D2,P1D2
  67.       PARAMETER (X1D2=0.5D0)
  68.       PARAMETER (Y1D2=0.5D0)
  69.       PARAMETER (P1D2=1.D0/6.D0)
  70. *
  71. * Pour le degré 3, méthode produit conique à 4 points : GPT2-3-4
  72. *
  73. *  Générateurs pour la cubature de degré 4 à 6 points : GAT2-4-6A :
  74. *   - [ Fully symmetric ]
  75.       REAL*8 X1D4,Y1D4,P1D4
  76.       PARAMETER (X1D4=0.0915762135097707434595714634022015D0)
  77.       PARAMETER (Y1D4=0.0915762135097707434595714634022015D0)
  78.       PARAMETER (P1D4=0.0549758718276609338191631624501052D0)
  79. *   - [ Fully symmetric ]
  80.       REAL*8 X2D4,Y2D4,P2D4
  81.       PARAMETER (X2D4=0.445948490915964886318329253883051D0)
  82.       PARAMETER (Y2D4=0.445948490915964886318329253883051D0)
  83.       PARAMETER (P2D4=0.111690794839005732847503504216561D0)
  84. *
  85. *  Générateurs pour la cubature de degré 5 à 7 points : GAT2-5-7 :
  86. *   - [ Fully symmetric ]
  87.       REAL*8 X1D5,Y1D5,P1D5
  88.       PARAMETER (X1D5=1.D0/3.D0)
  89.       PARAMETER (Y1D5=1.D0/3.D0)
  90.       PARAMETER (P1D5=0.1125D0)
  91. *   - [ Fully symmetric ]
  92.       REAL*8 X2D5,Y2D5,P2D5
  93.       PARAMETER (X2D5=0.101286507323456338800987361915123D0)
  94.       PARAMETER (Y2D5=0.101286507323456338800987361915123D0)
  95.       PARAMETER (P2D5=0.0629695902724135762978419727500906D0)
  96. *   - [ Fully symmetric ]
  97.       REAL*8 X3D5,Y3D5,P3D5
  98.       PARAMETER (X3D5=0.470142064105115089770441209513447D0)
  99.       PARAMETER (Y3D5=0.470142064105115089770441209513447D0)
  100.       PARAMETER (P3D5=0.0661970763942530903688246939165759D0)
  101. *
  102. * Pas de degré 6, même nb. de points que degré 7 avec la
  103. * qualité PI (coeffs positifs, points à l'intérieur).
  104. *
  105. *  Générateurs pour la cubature de degré 7 à 12 points : GAT2-7-12 :
  106. *   - [ Ro3 invariant ]
  107.       REAL*8 X1D7,Y1D7,P1D7
  108.       PARAMETER (X1D7=0.0623822650944021181736830009963499D0)
  109.       PARAMETER (Y1D7=0.0675178670739160854425571310508685D0)
  110.       PARAMETER (P1D7=0.0265170281574362514287541804607391D0)
  111. *   - [ Ro3 invariant ]
  112.       REAL*8 X2D7,Y2D7,P2D7
  113.       PARAMETER (X2D7=0.0552254566569266117374791902756449D0)
  114.       PARAMETER (Y2D7=0.321502493851981822666307849199202D0)
  115.       PARAMETER (P2D7=0.0438814087144460550367699031392875D0)
  116. *   - [ Ro3 invariant ]
  117.       REAL*8 X3D7,Y3D7,P3D7
  118.       PARAMETER (X3D7=0.0343243029450971464696306424839376D0)
  119.       PARAMETER (Y3D7=0.660949196186735657611980310197799D0)
  120.       PARAMETER (P3D7=0.0287750427849815857384454969002185D0)
  121. *   - [ Ro3 invariant ]
  122.       REAL*8 X4D7,Y4D7,P4D7
  123.       PARAMETER (X4D7=0.515842334353591779257463386826430D0)
  124.       PARAMETER (Y4D7=0.277716166976391782569581871393723D0)
  125.       PARAMETER (P4D7=0.0674931870098027744626970861664214D0)
  126. *
  127.       INTEGER NOPG,IPG
  128. *
  129. * Executable statements
  130. *
  131.       IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingatr'
  132. *
  133. * Méthode de nom : NCT2-1-3
  134. * Sur un triangle : cubature d'ordre 1 à 3 points
  135. *                   espace de référence de dimension 2
  136. *
  137. * In INIPG : SEGINI PGCOUR
  138.       CALL INIPG('NCT2-1-3','NEWTON-COTES','TRIANGLE',
  139.      $     1,3,2,
  140.      $     PGCOUR,
  141.      $     IMPR,IRET)
  142.       IF (IRET.NE.0) GOTO 9999
  143.       NOPG=0
  144.       XCOR(1)=1.D0
  145.       XCOR(2)=0.D0
  146.       XCOR(3)=0.D0
  147.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET)
  148.       IF (IRET.NE.0) GOTO 9999
  149.       SEGDES PGCOUR
  150.       MYPGS.LISPG(**)=PGCOUR
  151. *
  152. * Méthode de nom : NCT2-3-7
  153. * Sur un triangle : cubature d'ordre 3 à 7 points
  154. *                   espace de référence de dimension 2
  155. *
  156. * In INIPG : SEGINI PGCOUR
  157.       CALL INIPG('NCT2-3-7','NEWTON-COTES','TRIANGLE',
  158.      $     3,7,2,
  159.      $     PGCOUR,
  160.      $     IMPR,IRET)
  161.       IF (IRET.NE.0) GOTO 9999
  162.       NOPG=0
  163.       XCOR(1)=1.D0
  164.       XCOR(2)=0.D0
  165.       XCOR(3)=0.D0
  166.       POIDS=0.025D0
  167.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET)
  168.       IF (IRET.NE.0) GOTO 9999
  169.       XCOR(1)=0.5D0
  170.       XCOR(2)=0.5D0
  171.       XCOR(3)=0.D0
  172.       POIDS=0.2D0/3.D0
  173.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET)
  174.       IF (IRET.NE.0) GOTO 9999
  175.       XCOR(1)=1.D0/3.D0
  176.       XCOR(2)=1.D0/3.D0
  177.       XCOR(3)=1-XCOR(1)-XCOR(2)
  178.       POIDS=0.225D0
  179.       CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,POIDS,IMPR,IRET)
  180.       IF (IRET.NE.0) GOTO 9999
  181.  
  182.       SEGDES PGCOUR
  183.       MYPGS.LISPG(**)=PGCOUR
  184. *
  185. * Méthode de nom : GAT2-1-1
  186. * Sur un triangle : cubature d'ordre 1 à 1 point
  187. *                   espace de référence de dimension 2
  188. *
  189. * In INIPG : SEGINI PGCOUR
  190.       CALL INIPG('GAT2-1-1','GAUSS','TRIANGLE',
  191.      $     1,1,2,
  192.      $     PGCOUR,
  193.      $     IMPR,IRET)
  194.       IF (IRET.NE.0) GOTO 9999
  195.       NOPG=0
  196.       XCOR(1)=X1D1
  197.       XCOR(2)=Y1D1
  198.       XCOR(3)=1-XCOR(1)-XCOR(2)
  199.       CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D1,IMPR,IRET)
  200.       IF (IRET.NE.0) GOTO 9999
  201.       SEGDES PGCOUR
  202.       MYPGS.LISPG(**)=PGCOUR
  203. *
  204. * Méthode de nom : GAT2-2-3A
  205. * Sur un triangle : cubature d'ordre 2 à 3 points
  206. *                   espace de référence de dimension 2
  207. *
  208. * In INIPG : SEGINI PGCOUR
  209.       CALL INIPG('GAT2-2-3A','GAUSS','TRIANGLE',
  210.      $     2,3,2,
  211.      $     PGCOUR,
  212.      $     IMPR,IRET)
  213.       IF (IRET.NE.0) GOTO 9999
  214.       NOPG=0
  215.       XCOR(1)=X1D2
  216.       XCOR(2)=Y1D2
  217.       XCOR(3)=1-XCOR(1)-XCOR(2)
  218.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET)
  219.       IF (IRET.NE.0) GOTO 9999
  220.       SEGDES PGCOUR
  221.       MYPGS.LISPG(**)=PGCOUR
  222. *
  223. * Méthode de nom : GPT2-3-4
  224. * Sur un triangle : méthode gauss-produit conique d'ordre 3 à 4 points
  225. *                   espace de référence de dimension 2
  226. *
  227. * In INIPG : SEGINI PGCOUR
  228.       CALL INIPG('GPT2-3-4','GAUSS-PRODUIT','TRIANGLE',
  229.      $     3,4,2,
  230.      $     PGCOUR,
  231.      $     IMPR,IRET)
  232.       IF (IRET.NE.0) GOTO 9999
  233.       CALL FIPG('GAC1-3-2',MYPGS,PGPRO1,IMPR,IRET)
  234.       IF (IRET.NE.0) GOTO 9999
  235. *
  236. * Il faut remettre les points de Gauss-Legendre sur [0,1]
  237. * à la place de [-1,1]
  238. *
  239.       SEGINI,PGPRO2=PGPRO1
  240.       DO 1 IPG=1,PGPRO2.XCOPG(/2)
  241.          PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0
  242.          PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0
  243.  1    CONTINUE
  244.       SEGDES PGPRO2
  245.       PGPRO1=PGPRO2
  246.       CALL FIPG('GJ10-3-2',MYPGS,PGPRO2,IMPR,IRET)
  247.       IF (IRET.NE.0) GOTO 9999
  248.       CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET)
  249.       IF (IRET.NE.0) GOTO 9999
  250.       SEGDES PGCOUR
  251.       MYPGS.LISPG(**)=PGCOUR
  252. *
  253. * Méthode de nom : GAT2-4-6A
  254. * Sur un triangle : cubature d'ordre 4 à 6 points
  255. *                   espace de référence de dimension 2
  256. *
  257. * In INIPG : SEGINI PGCOUR
  258.       CALL INIPG('GAT2-4-6A','GAUSS','TRIANGLE',
  259.      $     4,6,2,
  260.      $     PGCOUR,
  261.      $     IMPR,IRET)
  262.       IF (IRET.NE.0) GOTO 9999
  263.       NOPG=0
  264.       XCOR(1)=X1D4
  265.       XCOR(2)=Y1D4
  266.       XCOR(3)=1-XCOR(1)-XCOR(2)
  267.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D4,IMPR,IRET)
  268.       IF (IRET.NE.0) GOTO 9999
  269.       XCOR(1)=X2D4
  270.       XCOR(2)=Y2D4
  271.       XCOR(3)=1-XCOR(1)-XCOR(2)
  272.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D4,IMPR,IRET)
  273.       IF (IRET.NE.0) GOTO 9999
  274.       SEGDES PGCOUR
  275.       MYPGS.LISPG(**)=PGCOUR
  276. *
  277. * Méthode de nom : GAT2-5-7
  278. * Sur un triangle : cubature d'ordre 5 à 7 points
  279. *                   espace de référence de dimension 2
  280. *
  281. * In INIPG : SEGINI PGCOUR
  282.       CALL INIPG('GAT2-5-7','GAUSS','TRIANGLE',
  283.      $     5,7,2,
  284.      $     PGCOUR,
  285.      $     IMPR,IRET)
  286.       IF (IRET.NE.0) GOTO 9999
  287.       NOPG=0
  288.       XCOR(1)=X1D5
  289.       XCOR(2)=Y1D5
  290.       XCOR(3)=1-XCOR(1)-XCOR(2)
  291.       CALL GTSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET)
  292.       IF (IRET.NE.0) GOTO 9999
  293.       XCOR(1)=X2D5
  294.       XCOR(2)=Y2D5
  295.       XCOR(3)=1-XCOR(1)-XCOR(2)
  296.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET)
  297.       IF (IRET.NE.0) GOTO 9999
  298.       XCOR(1)=X3D5
  299.       XCOR(2)=Y3D5
  300.       XCOR(3)=1-XCOR(1)-XCOR(2)
  301.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P3D5,IMPR,IRET)
  302.       IF (IRET.NE.0) GOTO 9999
  303.       SEGDES PGCOUR
  304.       MYPGS.LISPG(**)=PGCOUR
  305. *
  306. * Méthode de nom : GPT2-5-9
  307. * Sur un triangle : méthode gauss-produit conique d'ordre 5 à 9 points
  308. *                   espace de référence de dimension 2
  309. *
  310. * In INIPG : SEGINI PGCOUR
  311.       CALL INIPG('GPT2-5-9','GAUSS-PRODUIT','TRIANGLE',
  312.      $     5,9,2,
  313.      $     PGCOUR,
  314.      $     IMPR,IRET)
  315.       IF (IRET.NE.0) GOTO 9999
  316.       CALL FIPG('GAC1-5-3',MYPGS,PGPRO1,IMPR,IRET)
  317.       IF (IRET.NE.0) GOTO 9999
  318. *
  319. * Il faut remettre les points de Gauss-Legendre sur [0,1]
  320. * à la place de [-1,1]
  321. *
  322.       SEGINI,PGPRO2=PGPRO1
  323.       DO 3 IPG=1,PGPRO2.XCOPG(/2)
  324.          PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0
  325.          PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0
  326.  3    CONTINUE
  327.       SEGDES PGPRO2
  328.       PGPRO1=PGPRO2
  329.       CALL FIPG('GJ10-5-3',MYPGS,PGPRO2,IMPR,IRET)
  330.       IF (IRET.NE.0) GOTO 9999
  331.       CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET)
  332.       IF (IRET.NE.0) GOTO 9999
  333.       SEGDES PGCOUR
  334.       MYPGS.LISPG(**)=PGCOUR
  335. *
  336. * Méthode de nom : GAT2-7-12
  337. * Sur un triangle : cubature d'ordre 7 à 12 points
  338. *                   espace de référence de dimension 2
  339. *
  340. * In INIPG : SEGINI PGCOUR
  341.       CALL INIPG('GAT2-7-12','GAUSS','TRIANGLE',
  342.      $     7,12,2,
  343.      $     PGCOUR,
  344.      $     IMPR,IRET)
  345.       IF (IRET.NE.0) GOTO 9999
  346.       NOPG=0
  347.       XCOR(1)=X1D7
  348.       XCOR(2)=Y1D7
  349.       XCOR(3)=1-XCOR(1)-XCOR(2)
  350.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET)
  351.       IF (IRET.NE.0) GOTO 9999
  352.       XCOR(1)=X2D7
  353.       XCOR(2)=Y2D7
  354.       XCOR(3)=1-XCOR(1)-XCOR(2)
  355.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET)
  356.       IF (IRET.NE.0) GOTO 9999
  357.       XCOR(1)=X3D7
  358.       XCOR(2)=Y3D7
  359.       XCOR(3)=1-XCOR(1)-XCOR(2)
  360.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P3D7,IMPR,IRET)
  361.       IF (IRET.NE.0) GOTO 9999
  362.       XCOR(1)=X4D7
  363.       XCOR(2)=Y4D7
  364.       XCOR(3)=1-XCOR(1)-XCOR(2)
  365.       CALL GTRO3I(PGCOUR,NOPG,DIMSRF,XCOR,P4D7,IMPR,IRET)
  366.       IF (IRET.NE.0) GOTO 9999
  367.       SEGDES PGCOUR
  368.       MYPGS.LISPG(**)=PGCOUR
  369. *
  370. * Méthode de nom : GPT2-7-16
  371. * Sur un triangle : méthode gauss-produit conique d'ordre 7 à 16 points
  372. *                   espace de référence de dimension 2
  373. *
  374. * In INIPG : SEGINI PGCOUR
  375.       CALL INIPG('GPT2-7-16','GAUSS-PRODUIT','TRIANGLE',
  376.      $     7,16,2,
  377.      $     PGCOUR,
  378.      $     IMPR,IRET)
  379.       IF (IRET.NE.0) GOTO 9999
  380.       CALL FIPG('GAC1-7-4',MYPGS,PGPRO1,IMPR,IRET)
  381.       IF (IRET.NE.0) GOTO 9999
  382. *
  383. * Il faut remettre les points de Gauss-Legendre sur [0,1]
  384. * à la place de [-1,1]
  385. *
  386.       SEGINI,PGPRO2=PGPRO1
  387.       DO 5 IPG=1,PGPRO2.XCOPG(/2)
  388.          PGPRO2.XCOPG(1,IPG)=(1.D0+PGPRO2.XCOPG(1,IPG))*0.5D0
  389.          PGPRO2.XPOPG(IPG)=PGPRO2.XPOPG(IPG)*0.5D0
  390.  5    CONTINUE
  391.       SEGDES PGPRO2
  392.       PGPRO1=PGPRO2
  393.       CALL FIPG('GJ10-7-4',MYPGS,PGPRO2,IMPR,IRET)
  394.       IF (IRET.NE.0) GOTO 9999
  395.       CALL CPROPG(PGPRO1,PGPRO2,PGCOUR,IMPR,IRET)
  396.       IF (IRET.NE.0) GOTO 9999
  397.       SEGDES PGCOUR
  398.       MYPGS.LISPG(**)=PGCOUR
  399. *
  400. * Normal termination
  401. *
  402.       IRET=0
  403.       RETURN
  404. *
  405. * Format handling
  406. *
  407. *
  408. * Error handling
  409. *
  410.  9999 CONTINUE
  411.       IRET=1
  412.       WRITE(IOIMP,*) 'An error was detected in subroutine ingatr'
  413.       RETURN
  414. *
  415. * End of subroutine INGATR
  416. *
  417.       END
  418.  
  419.  
  420.  
  421.  
  422.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales