Télécharger ingatr.eso

Retour à la liste

Numérotation des lignes :

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

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