Télécharger ingaqu.eso

Retour à la liste

Numérotation des lignes :

  1. C INGAQU SOURCE PV 07/11/23 21:17:18 5978
  2. SUBROUTINE INGAQU(MYPGS,IMPR,IRET)
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM : INGAQU
  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 carré (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, GCSINO, GCPASY, GCROIN, GCCESY, GCRESY
  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. *
  53. integer PGPRO1
  54. INTEGER IMPR,IRET
  55. *
  56. INTEGER DIMSRF
  57. PARAMETER(DIMSRF=2)
  58. REAL*8 XCOR(DIMSRF)
  59. *
  60. * Générateurs pour la cubature de degré 1 à 1 point : GAC2-1-1 :
  61. * - [ Origin ]
  62. REAL*8 X0D1,Y0D1,P0D1
  63. PARAMETER (X0D1= 0.D0)
  64. PARAMETER (Y0D1= 0.D0)
  65. PARAMETER (P0D1= 4.D0)
  66. *
  67. * Générateurs pour la cubature de degré 2 à 3 points : GAC2-2-3 :
  68. * - [ Partial symmetry ] (\sqrt{2./3.} ; 0.) ; (4./3.)
  69. REAL*8 X1D2,Y1D2,P1D2
  70. PARAMETER (X1D2= 0.81649658092772603273242802490196D0)
  71. PARAMETER (Y1D2= 0.D0)
  72. PARAMETER (P1D2= 4.D0/3.D0)
  73. * - [ Partial symmetry ] (-\sqrt{1./6.} ; \sqrt{1./2.}) ; (4./3.)
  74. REAL*8 X2D2,Y2D2,P2D2
  75. PARAMETER (X2D2=-0.40824829046386301636621401245098D0)
  76. PARAMETER (Y2D2= 0.70710678118654752440084436210485D0)
  77. PARAMETER (P2D2= 4.D0/3.D0)
  78. *
  79. * Générateurs pour la cubature de degré 3 à 4 points : GAC2-3-4A :
  80. * - [ Fully symmetric ] (\sqrt{2./3.} ; 0.) ; (1.)
  81. REAL*8 X1D3,Y1D3,P1D3
  82. PARAMETER (X1D3= 0.816496580927726032732428024901963D0)
  83. PARAMETER (Y1D3= 0.D0)
  84. PARAMETER (P1D3= 1.D0)
  85. *
  86. * Générateurs pour la cubature de degré 4 à 6 points : GAC2-4-6C :
  87. * - [ Origin ]
  88. REAL*8 X0D4,Y0D4,P0D4
  89. PARAMETER (X0D4= 0.D0)
  90. PARAMETER (Y0D4= 0.D0)
  91. PARAMETER (P0D4= 1.14285714285714285714285714285714D0)
  92. * - [ Partial symmetry ]
  93. REAL*8 X1D4,Y1D4,P1D4
  94. PARAMETER (X1D4= 0.966091783079295904914577610477984D0)
  95. PARAMETER (Y1D4= 0.D0)
  96. PARAMETER (P1D4= 0.439560439560439560439560439560439D0)
  97. * - [ Partial symmetry ]
  98. REAL*8 X2D4,Y2D4,P2D4
  99. PARAMETER (X2D4= 0.455603727836192842249584745989760D0)
  100. PARAMETER (Y2D4= 0.851914653304600491301835640366669D0)
  101. PARAMETER (P2D4= 0.566072207007532105264050459451812D0)
  102. * - [ Partial symmetry ]
  103. REAL*8 X3D4,Y3D4,P3D4
  104. PARAMETER (X3D4=-0.731629951573134529368035491840613D0)
  105. PARAMETER (Y3D4= 0.630912788976754026243413202993658D0)
  106. PARAMETER (P3D4= 0.642719001783676685944740749339395D0)
  107. *
  108. * Générateurs pour la cubature de degré 5 à 7 points : GAC2-5-7A :
  109. * - [ Origin ]
  110. REAL*8 X0D5,Y0D5,P0D5
  111. PARAMETER (X0D5= 0.D0)
  112. PARAMETER (Y0D5= 0.D0)
  113. PARAMETER (P0D5= 1.14285714285714285714285714285714D0)
  114. * - [ Central symmetry ]
  115. REAL*8 X1D5,Y1D5,P1D5
  116. PARAMETER (X1D5= 0.966091783079295904914577610477984D0)
  117. PARAMETER (Y1D5= 0.D0)
  118. PARAMETER (P1D5= 0.317460317460317460317460317460317D0)
  119. * - [ Rectangular symmetry ]
  120. REAL*8 X2D5,Y2D5,P2D5
  121. PARAMETER (X2D5= 0.577350269189625764509148780501957D0)
  122. PARAMETER (Y2D5= 0.774596669241483377035853079956479D0)
  123. PARAMETER (P2D5= 0.555555555555555555555555555555555D0)
  124. *
  125. * Générateurs pour la cubature de degré 6 à 10 points : GAC2-6-10C :
  126. * - [ Partial symmetry ]
  127. REAL*8 X1D6,Y1D6,P1D6
  128. PARAMETER (X1D6= 0.836405633697625604537901440708439D0)
  129. PARAMETER (Y1D6= 0.D0)
  130. PARAMETER (P1D6= 0.455343245714173438553541721312115D0)
  131. * - [ Partial symmetry ]
  132. REAL*8 X2D6,Y2D6,P2D6
  133. PARAMETER (X2D6=-0.357460165391307184886174384445965D0)
  134. PARAMETER (Y2D6= 0.D0)
  135. PARAMETER (P2D6= 0.827395973202965492356696938299033D0)
  136. * - [ Partial symmetry ]
  137. REAL*8 X3D6,Y3D6,P3D6
  138. PARAMETER (X3D6= 0.872101531193130590032875340374412D0)
  139. PARAMETER (Y3D6= 0.888764014654764533853170471609380D0)
  140. PARAMETER (P3D6= 0.144000884599645384928603086560623D0)
  141. * - [ Partial symmetry ]
  142. REAL*8 X4D6,Y4D6,P4D6
  143. PARAMETER (X4D6= 0.305985162155426661489382735682569D0)
  144. PARAMETER (Y4D6= 0.604857639464685027698532295335587D0)
  145. PARAMETER (P4D6= 0.668259104262665149929832544411851D0)
  146. * - [ Partial symmetry ]
  147. REAL*8 X5D6,Y5D6,P5D6
  148. PARAMETER (X5D6=-0.410270899466657838725680144555973D0)
  149. PARAMETER (Y5D6= 0.955447506641063747757528360774788D0)
  150. PARAMETER (P5D6= 0.225474004890679357411055647855252D0)
  151. * - [ Partial symmetry ]
  152. REAL*8 X6D6,Y6D6,P6D6
  153. PARAMETER (X6D6=-0.872869311156879339251259208496329D0)
  154. PARAMETER (Y6D6= 0.565459993438754045894375389139462D0)
  155. PARAMETER (P6D6= 0.320896396788440642275389391366698D0)
  156. *
  157. * Générateurs pour la cubature de degré 7 à 12 points : GAC2-7-12A :
  158. * - [ Fully symmetric ]
  159. REAL*8 X1D7,Y1D7,P1D7
  160. PARAMETER (X1D7= 0.925820099772551461566566776583999D0)
  161. PARAMETER (Y1D7= 0.D0)
  162. PARAMETER (P1D7= 0.241975308641975308641975308641975D0)
  163. * - [ Fully symmetric ]
  164. REAL*8 X2D7,Y2D7,P2D7
  165. PARAMETER (X2D7= 0.380554433208315656379106359086394D0)
  166. PARAMETER (Y2D7= 0.380554433208315656379106359086394D0)
  167. PARAMETER (P2D7= 0.520592916667394457139919432046731D0)
  168. * - [ Fully symmetric ]
  169. REAL*8 X3D7,Y3D7,P3D7
  170. PARAMETER (X3D7= 0.805979782918598743707856181350744D0)
  171. PARAMETER (Y3D7= 0.805979782918598743707856181350744D0)
  172. PARAMETER (P3D7= 0.237431774690630234218105259311293D0)
  173. *
  174. INTEGER NOPG
  175. *
  176. * Executable statements
  177. *
  178. IF (IMPR.GT.6) WRITE(IOIMP,*) 'Entrée dans ingaqu'
  179. *
  180. * Méthode de nom : NCC2-1-4
  181. * Sur un carré : cubature d'ordre 1 à 4 points
  182. * espace de référence de dimension 2
  183. *
  184. * In INIPG : SEGINI PGCOUR
  185. CALL INIPG('NCC2-1-4','NEWTON-COTES','CARRE',
  186. $ 1,4,2,
  187. $ PGCOUR,
  188. $ IMPR,IRET)
  189. IF (IRET.NE.0) GOTO 9999
  190. NOPG=0
  191. XCOR(1)=1.D0
  192. XCOR(2)=1.D0
  193. CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET)
  194. IF (IRET.NE.0) GOTO 9999
  195. SEGDES PGCOUR
  196. MYPGS.LISPG(**)=PGCOUR
  197. *
  198. * Méthode de nom : NCC2-3-9
  199. * Sur un carré : cubature d'ordre 3 à 9 points
  200. * espace de référence de dimension 2
  201. *
  202. * In INIPG : SEGINI PGCOUR
  203. CALL INIPG('NCC2-3-9','NEWTON-COTES','CARRE',
  204. $ 3,9,2,
  205. $ PGCOUR,
  206. $ IMPR,IRET)
  207. IF (IRET.NE.0) GOTO 9999
  208. CALL FIPG('NCC1-3-3',MYPGS,PGPRO1,IMPR,IRET)
  209. IF (IRET.NE.0) GOTO 9999
  210. CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET)
  211. IF (IRET.NE.0) GOTO 9999
  212. SEGDES PGCOUR
  213. MYPGS.LISPG(**)=PGCOUR
  214. *
  215. * Méthode de nom : GAC2-1-1
  216. * Sur un carré : cubature d'ordre 1 à 1 point
  217. * espace de référence de dimension 2
  218. *
  219. * In INIPG : SEGINI PGCOUR
  220. CALL INIPG('GAC2-1-1','GAUSS','CARRE',
  221. $ 1,1,2,
  222. $ PGCOUR,
  223. $ IMPR,IRET)
  224. IF (IRET.NE.0) GOTO 9999
  225. NOPG=0
  226. XCOR(1)=X0D1
  227. XCOR(2)=Y0D1
  228. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D1,IMPR,IRET)
  229. IF (IRET.NE.0) GOTO 9999
  230. SEGDES PGCOUR
  231. MYPGS.LISPG(**)=PGCOUR
  232. *
  233. * Méthode de nom : GAC2-2-3
  234. * Sur un carré : cubature d'ordre 2 à 3 points
  235. * espace de référence de dimension 2
  236. *
  237. * In INIPG : SEGINI PGCOUR
  238. CALL INIPG('GAC2-2-3','GAUSS','CARRE',
  239. $ 2,3,2,
  240. $ PGCOUR,
  241. $ IMPR,IRET)
  242. IF (IRET.NE.0) GOTO 9999
  243. NOPG=0
  244. XCOR(1)=X1D2
  245. XCOR(2)=Y1D2
  246. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D2,IMPR,IRET)
  247. IF (IRET.NE.0) GOTO 9999
  248. XCOR(1)=X2D2
  249. XCOR(2)=Y2D2
  250. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P2D2,IMPR,IRET)
  251. IF (IRET.NE.0) GOTO 9999
  252. SEGDES PGCOUR
  253. MYPGS.LISPG(**)=PGCOUR
  254. *
  255. * Méthode de nom : GAC2-3-4A
  256. * Sur un carré : cubature d'ordre 3 à 4 points
  257. * espace de référence de dimension 2
  258. *
  259. * In INIPG : SEGINI PGCOUR
  260. CALL INIPG('GAC2-3-4A','GAUSS','CARRE',
  261. $ 3,4,2,
  262. $ PGCOUR,
  263. $ IMPR,IRET)
  264. IF (IRET.NE.0) GOTO 9999
  265. NOPG=0
  266. XCOR(1)=X1D3
  267. XCOR(2)=Y1D3
  268. CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D3,IMPR,IRET)
  269. IF (IRET.NE.0) GOTO 9999
  270. SEGDES PGCOUR
  271. MYPGS.LISPG(**)=PGCOUR
  272. *
  273. * Méthode de nom : GPC2-3-4
  274. * Sur un carré : méthode gauss-produit d'ordre 3 à 4 points
  275. * espace de référence de dimension 2
  276. *
  277. * In INIPG : SEGINI PGCOUR
  278. CALL INIPG('GPC2-3-4','GAUSS-PRODUIT','CARRE',
  279. $ 3,4,2,
  280. $ PGCOUR,
  281. $ IMPR,IRET)
  282. IF (IRET.NE.0) GOTO 9999
  283. CALL FIPG('GAC1-3-2',MYPGS,PGPRO1,IMPR,IRET)
  284. IF (IRET.NE.0) GOTO 9999
  285. CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET)
  286. IF (IRET.NE.0) GOTO 9999
  287. SEGDES PGCOUR
  288. MYPGS.LISPG(**)=PGCOUR
  289. *
  290. * Méthode de nom : GAC2-4-6C
  291. * Sur un carré : cubature d'ordre 4 à 6 points
  292. * espace de référence de dimension 2
  293. *
  294. * In INIPG : SEGINI PGCOUR
  295. CALL INIPG('GAC2-4-6C','GAUSS','CARRE',
  296. $ 4,6,2,
  297. $ PGCOUR,
  298. $ IMPR,IRET)
  299. IF (IRET.NE.0) GOTO 9999
  300. NOPG=0
  301. XCOR(1)=X0D4
  302. XCOR(2)=Y0D4
  303. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D4,IMPR,IRET)
  304. IF (IRET.NE.0) GOTO 9999
  305. XCOR(1)=X1D4
  306. XCOR(2)=Y1D4
  307. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D4,IMPR,IRET)
  308. IF (IRET.NE.0) GOTO 9999
  309. XCOR(1)=X2D4
  310. XCOR(2)=Y2D4
  311. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P2D4,IMPR,IRET)
  312. IF (IRET.NE.0) GOTO 9999
  313. XCOR(1)=X3D4
  314. XCOR(2)=Y3D4
  315. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P3D4,IMPR,IRET)
  316. IF (IRET.NE.0) GOTO 9999
  317. SEGDES PGCOUR
  318. MYPGS.LISPG(**)=PGCOUR
  319. *
  320. * Méthode de nom : GAC2-5-7A
  321. * Sur un carré : cubature d'ordre 5 à 7 points
  322. * espace de référence de dimension 2
  323. *
  324. * In INIPG : SEGINI PGCOUR
  325. CALL INIPG('GAC2-5-7A','GAUSS','CARRE',
  326. $ 5,7,2,
  327. $ PGCOUR,
  328. $ IMPR,IRET)
  329. IF (IRET.NE.0) GOTO 9999
  330. NOPG=0
  331. XCOR(1)=X0D5
  332. XCOR(2)=Y0D5
  333. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P0D5,IMPR,IRET)
  334. IF (IRET.NE.0) GOTO 9999
  335. XCOR(1)=X1D5
  336. XCOR(2)=Y1D5
  337. CALL GCCESY(PGCOUR,NOPG,DIMSRF,XCOR,P1D5,IMPR,IRET)
  338. IF (IRET.NE.0) GOTO 9999
  339. XCOR(1)=X2D5
  340. XCOR(2)=Y2D5
  341. CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D5,IMPR,IRET)
  342. IF (IRET.NE.0) GOTO 9999
  343. SEGDES PGCOUR
  344. MYPGS.LISPG(**)=PGCOUR
  345. *
  346. * Méthode de nom : GPC2-5-9
  347. * Sur un carré : méthode gauss-produit d'ordre 5 à 9 points
  348. * espace de référence de dimension 2
  349. *
  350. * In INIPG : SEGINI PGCOUR
  351. CALL INIPG('GPC2-5-9','GAUSS-PRODUIT','CARRE',
  352. $ 5,9,2,
  353. $ PGCOUR,
  354. $ IMPR,IRET)
  355. IF (IRET.NE.0) GOTO 9999
  356. CALL FIPG('GAC1-5-3',MYPGS,PGPRO1,IMPR,IRET)
  357. IF (IRET.NE.0) GOTO 9999
  358. CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET)
  359. IF (IRET.NE.0) GOTO 9999
  360. SEGDES PGCOUR
  361. MYPGS.LISPG(**)=PGCOUR
  362. *
  363. * Méthode de nom : GAC2-6-10C
  364. * Sur un carré : cubature d'ordre 6 à 10 points
  365. * espace de référence de dimension 2
  366. *
  367. * In INIPG : SEGINI PGCOUR
  368. CALL INIPG('GAC2-6-10C','GAUSS','CARRE',
  369. $ 6,10,2,
  370. $ PGCOUR,
  371. $ IMPR,IRET)
  372. IF (IRET.NE.0) GOTO 9999
  373. NOPG=0
  374. XCOR(1)=X1D6
  375. XCOR(2)=Y1D6
  376. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P1D6,IMPR,IRET)
  377. IF (IRET.NE.0) GOTO 9999
  378. XCOR(1)=X2D6
  379. XCOR(2)=Y2D6
  380. CALL GCSINO(PGCOUR,NOPG,DIMSRF,XCOR,P2D6,IMPR,IRET)
  381. IF (IRET.NE.0) GOTO 9999
  382. XCOR(1)=X3D6
  383. XCOR(2)=Y3D6
  384. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P3D6,IMPR,IRET)
  385. IF (IRET.NE.0) GOTO 9999
  386. XCOR(1)=X4D6
  387. XCOR(2)=Y4D6
  388. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P4D6,IMPR,IRET)
  389. IF (IRET.NE.0) GOTO 9999
  390. XCOR(1)=X5D6
  391. XCOR(2)=Y5D6
  392. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P5D6,IMPR,IRET)
  393. IF (IRET.NE.0) GOTO 9999
  394. XCOR(1)=X6D6
  395. XCOR(2)=Y6D6
  396. CALL GCPASY(PGCOUR,NOPG,DIMSRF,XCOR,P6D6,IMPR,IRET)
  397. IF (IRET.NE.0) GOTO 9999
  398. SEGDES PGCOUR
  399. MYPGS.LISPG(**)=PGCOUR
  400. *
  401. * Méthode de nom : GAC2-7-12A
  402. * Sur un carré : cubature d'ordre 7 à 12 points
  403. * espace de référence de dimension 2
  404. *
  405. * In INIPG : SEGINI PGCOUR
  406. CALL INIPG('GAC2-7-12A','GAUSS','CARRE',
  407. $ 7,12,2,
  408. $ PGCOUR,
  409. $ IMPR,IRET)
  410. IF (IRET.NE.0) GOTO 9999
  411. NOPG=0
  412. XCOR(1)=X1D7
  413. XCOR(2)=Y1D7
  414. CALL GCROIN(PGCOUR,NOPG,DIMSRF,XCOR,P1D7,IMPR,IRET)
  415. IF (IRET.NE.0) GOTO 9999
  416. XCOR(1)=X2D7
  417. XCOR(2)=Y2D7
  418. CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P2D7,IMPR,IRET)
  419. IF (IRET.NE.0) GOTO 9999
  420. XCOR(1)=X3D7
  421. XCOR(2)=Y3D7
  422. CALL GCRESY(PGCOUR,NOPG,DIMSRF,XCOR,P3D7,IMPR,IRET)
  423. IF (IRET.NE.0) GOTO 9999
  424. SEGDES PGCOUR
  425. MYPGS.LISPG(**)=PGCOUR
  426. *
  427. * Méthode de nom : GPC2-7-16
  428. * Sur un carré : méthode gauss-produit d'ordre 7 à 16 points
  429. * espace de référence de dimension 2
  430. *
  431. * In INIPG : SEGINI PGCOUR
  432. CALL INIPG('GPC2-7-16','GAUSS-PRODUIT','CARRE',
  433. $ 7,16,2,
  434. $ PGCOUR,
  435. $ IMPR,IRET)
  436. IF (IRET.NE.0) GOTO 9999
  437. CALL FIPG('GAC1-7-4',MYPGS,PGPRO1,IMPR,IRET)
  438. IF (IRET.NE.0) GOTO 9999
  439. CALL PROPG(PGPRO1,PGPRO1,PGCOUR,IMPR,IRET)
  440. IF (IRET.NE.0) GOTO 9999
  441. SEGDES PGCOUR
  442. MYPGS.LISPG(**)=PGCOUR
  443. *
  444. * Normal termination
  445. *
  446. IRET=0
  447. RETURN
  448. *
  449. * Format handling
  450. *
  451. *
  452. * Error handling
  453. *
  454. 9999 CONTINUE
  455. IRET=1
  456. WRITE(IOIMP,*) 'An error was detected in subroutine ingaqu'
  457. RETURN
  458. *
  459. * End of subroutine INGAQU
  460. *
  461. END
  462.  
  463.  
  464.  
  465.  
  466.  

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