Télécharger pente3.eso

Retour à la liste

Numérotation des lignes :

  1. C PENTE3 SOURCE CHAT 05/01/13 02:12:24 5004
  2. SUBROUTINE PENTE3(NFAC,MELEFL,MLECEN,NCOMP,MPOCHP,
  3. & MPOGRA,MPOMIN,MPOMAX,MPOALP)
  4. C
  5. C************************************************************************
  6. C
  7. C PROJET : CASTEM 2000
  8. C
  9. C NOM : PENTE3
  10. C
  11. C DESCRIPTION : Cette subroutine est appellée par la subroutine
  12. C PENTE1 (calcul du gradient d'un CHPOINT 2D de type
  13. C CENTRE)
  14. C Elle contient la partie du calcul du limiteur.
  15. C
  16. C LANGUAGE : ESOPE 2000 avec extensions CISI
  17. C
  18. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
  19. C AUTEUR (Modif.) : R. MOREL, DRN/DMT/SEMT/TTMF
  20. C
  21. C************************************************************************
  22. C
  23. C
  24. C APPELES (E/S) : none
  25. C
  26. C APPELES (Calcul) : none
  27. C
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES : NFAC : nombre de faces
  32. C
  33. C MELEFL : pointeur du MELEME 'FACEL'
  34. C
  35. C MLECEN : pointeur de MLENTI qui contient la table
  36. C numerotation global/local de CENTREs
  37. C
  38. C NCOMP : nombre de composantes de CHPOINT dont on veut
  39. C calculer les pentes
  40. C
  41. C MPOCHP : pointeur de MPOVAL de CHPOINT dont on veut le
  42. C gradient
  43. C
  44. C MPOGRA : pointeur de MPOVAL du gradient
  45. C
  46. C MPOMIN : pointeur de MPOVAL du minimum sur le stencil
  47. C
  48. C MPOMAX : pointeur de MPOVAL du maximum sur le stencil
  49. C
  50. C SORTIES : MPOALP : pointeur de MELVAL du limiteur
  51. C
  52. C************************************************************************
  53. C
  54. C HISTORIQUE (Anomalies et modifications éventuelles)
  55. C
  56. C HISTORIQUE : Cree le 6-4-1998
  57. C
  58. C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D
  59. C
  60. C************************************************************************
  61. C
  62. IMPLICIT INTEGER(I-N)
  63. -INC CCOPTIO
  64. C
  65. C**** Variables de CCOPTIO
  66. C
  67. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  68. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  69. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  70. C & ,IECHO, IIMPI, IOSPI
  71. C & ,IDIM
  72. CC & ,MCOORD
  73. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  74. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  75. C & ,NORINC,NORVAL,NORIND,NORVAD
  76. C & ,NUCROU, IPSAUV
  77. C
  78. INTEGER NFAC,NCOMP
  79. & ,NLCF,NGCF,NGCEG,NLCEG,NGCED,NLCED
  80. & ,INDCEL,I1,IGR,IDIMP1
  81. C
  82. C**** Pour l'instant 2D
  83. C
  84. REAL*8 XG,YG,XD,YD,XC,YC,DXG,DYG,DXD,DYD,DG,DD,DT
  85. & ,PHI,PHIMAX,PHIMIN,DPHIST,GRADX,GRADY,DPHI0,DPHI1,ALPHA
  86. & ,VALEUR
  87. C
  88. C**** Extension 3D
  89. C
  90. REAL*8 ZG,ZD,ZC,DZG,DZD,GRADZ
  91. C
  92. C
  93. -INC SMCOORD
  94. -INC SMCHPOI
  95. -INC SMELEME
  96. -INC SMLENTI
  97. POINTEUR MPOMIN.MPOVAL, MPOMAX.MPOVAL, MPOALP.MPOVAL,
  98. & MPOCHP.MPOVAL, MPOGRA.MPOVAL
  99. POINTEUR MELEFL.MELEME, MLECEN.MLENTI
  100. C
  101. C**** N.B. Tous les pointeurs ici sont déjà activés!
  102. C
  103. C**** Maillage FACEL
  104. C Boucle sur les faces
  105. C
  106. C
  107. C**** NGCEG = Numero global centre d'elt "gauche"
  108. C NLCEG = Numero local centre d'elt "gauche"
  109. C NGCF = " global centre de face
  110. C NLCF = numero local centre de face
  111. C NGCEG = Numero global centre d'elt "gauche"
  112. C NLCEG = Numero local centre d'elt "gauche"
  113. C
  114. IDIMP1 = IDIM+1
  115. C
  116. C Cas de la dimension 2
  117. C
  118. IF (IDIM .EQ. 2) THEN
  119. DO NLCF = 1, NFAC
  120. NGCEG = MELEFL.NUM(1,NLCF)
  121. NLCEG = MLECEN.LECT(NGCEG)
  122. NGCF = MELEFL.NUM(2,NLCF)
  123. NGCED = MELEFL.NUM(3,NLCF)
  124. NLCED = MLECEN.LECT(NGCED)
  125. IF(NGCEG .EQ. NGCED)THEN
  126. C
  127. C********** Limiteur dans le cas mur
  128. C
  129. C
  130. C********** Coordonees et parametres geometriques
  131. C
  132. C
  133. C XCOOR : table de coordonnes de points;
  134. C pour le point de numero global NG
  135. C X_NG = XCOOR((NG-1)*(IDIM+1)+1)
  136. C Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
  137. C Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
  138. C (dans notre cas IDIM=2, i.e. 2D)
  139. C
  140. INDCEL = (NGCEG-1)*IDIMP1
  141. XG = XCOOR(INDCEL+1)
  142. YG = XCOOR(INDCEL+2)
  143. INDCEL = (NGCF-1)*IDIMP1
  144. XC = XCOOR(INDCEL+1)
  145. YC = XCOOR(INDCEL+2)
  146. DXG = XC - XG
  147. DYG = YC - YG
  148. C
  149. C********** Boucle sur le composantes
  150. C
  151. DO I1 = 1, NCOMP
  152. IGR = 2*I1 - 1
  153. C
  154. C************* Limiteur
  155. C
  156. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  157. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  158. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  159. DPHIST = PHIMAX - PHIMIN
  160. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  161. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  162. DPHI0 = GRADX * DXG + GRADY * DYG
  163. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  164. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  165. VALEUR = 1.0D0
  166. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  167. DPHI1 = PHIMAX - PHI
  168. VALEUR = DPHI1/DPHI0*0.5D0
  169. ELSE
  170. DPHI1 = PHIMIN - PHI
  171. VALEUR = DPHI1/DPHI0*0.5D0
  172. ENDIF
  173. ALPHA = MIN(ALPHA, VALEUR)
  174. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  175. ENDDO
  176. ELSE
  177. C
  178. C******* NGCEG != NGCED
  179. C
  180. INDCEL = (NGCEG-1)*IDIMP1
  181. XG = XCOOR(INDCEL+1)
  182. YG = XCOOR(INDCEL+2)
  183. INDCEL = (NGCED-1)*IDIMP1
  184. XD = XCOOR(INDCEL+1)
  185. YD = XCOOR(INDCEL+2)
  186. INDCEL = (NGCF-1)*IDIMP1
  187. XC = XCOOR(INDCEL+1)
  188. YC = XCOOR(INDCEL+2)
  189. C
  190. DXG = XC - XG
  191. DYG = YC - YG
  192. DXD = XC - XD
  193. DYD = YC - YD
  194. DG = DXG * DXG + DYG * DYG
  195. DG = SQRT(DG)
  196. DD = DXD * DXD + DYD * DYD
  197. DD = SQRT(DD)
  198. DT = DG + DD
  199. C
  200. C********** Boucle sur le composantes
  201. C
  202. DO I1 = 1, NCOMP
  203. IGR = 2*I1 - 1
  204. C
  205. C************* Limiteur a gauche
  206. C
  207. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  208. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  209. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  210. DPHIST = PHIMAX - PHIMIN
  211. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  212. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  213. DPHI0 = GRADX * DXG + GRADY * DYG
  214. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  215. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  216. VALEUR = 1.0D0
  217. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  218. DPHI1 = PHIMAX - PHI
  219. VALEUR = DPHI1/DPHI0*DG/DT
  220. ELSE
  221. DPHI1 = PHIMIN - PHI
  222. VALEUR = DPHI1/DPHI0*DG/DT
  223. ENDIF
  224.  
  225. ALPHA = MIN(ALPHA, VALEUR)
  226. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  227. C
  228. C************* Limiteur a droite
  229. C
  230. PHI = MPOCHP.VPOCHA(NLCED,I1)
  231. PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
  232. PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
  233. DPHIST = PHIMAX - PHIMIN
  234. GRADX = MPOGRA.VPOCHA(NLCED,IGR)
  235. GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
  236. DPHI0 = GRADX * DXD + GRADY * DYD
  237. ALPHA = MPOALP.VPOCHA(NLCED,I1)
  238. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  239. VALEUR = 1.0D0
  240. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  241. DPHI1 = PHIMAX - PHI
  242. VALEUR = DPHI1/DPHI0*DD/DT
  243. ELSE
  244. DPHI1 = PHIMIN - PHI
  245. VALEUR = DPHI1/DPHI0*DD/DT
  246. ENDIF
  247. ALPHA = MIN(ALPHA, VALEUR)
  248. MPOALP.VPOCHA(NLCED,I1) = ALPHA
  249. ENDDO
  250. ENDIF
  251. ENDDO
  252. C
  253. C Cas de la dimension 3
  254. C
  255. ELSE
  256. DO NLCF = 1, NFAC
  257. NGCEG = MELEFL.NUM(1,NLCF)
  258. NLCEG = MLECEN.LECT(NGCEG)
  259. NGCF = MELEFL.NUM(2,NLCF)
  260. NGCED = MELEFL.NUM(3,NLCF)
  261. NLCED = MLECEN.LECT(NGCED)
  262. IF(NGCEG .EQ. NGCED)THEN
  263. C
  264. C********** Limiteur dans le cas mur
  265. C
  266. C
  267. C********** Coordonees et parametres geometriques
  268. C
  269. C
  270. C XCOOR : table de coordonnes de points;
  271. C pour le point de numero global NG
  272. C X_NG = XCOOR((NG-1)*(IDIM+1)+1)
  273. C Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
  274. C Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
  275. C
  276. INDCEL = (NGCEG-1)*IDIMP1
  277. XG = XCOOR(INDCEL+1)
  278. YG = XCOOR(INDCEL+2)
  279. ZG = XCOOR(INDCEL+3)
  280. INDCEL = (NGCF-1)*IDIMP1
  281. XC = XCOOR(INDCEL+1)
  282. YC = XCOOR(INDCEL+2)
  283. ZC = XCOOR(INDCEL+3)
  284. DXG = XC - XG
  285. DYG = YC - YG
  286. DZG = ZC - ZG
  287. C
  288. C********** Boucle sur le composantes
  289. C
  290. DO I1 = 1, NCOMP
  291. IGR = IDIM*(I1-1)+1
  292. C
  293. C************* Limiteur
  294. C
  295. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  296. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  297. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  298. DPHIST = PHIMAX - PHIMIN
  299. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  300. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  301. GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
  302. DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
  303. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  304. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  305. VALEUR = 1.0D0
  306. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  307. DPHI1 = PHIMAX - PHI
  308. VALEUR = DPHI1/DPHI0*0.5D0
  309. ELSE
  310. DPHI1 = PHIMIN - PHI
  311. VALEUR = DPHI1/DPHI0*0.5D0
  312. ENDIF
  313. ALPHA = MIN(ALPHA, VALEUR)
  314. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  315. ENDDO
  316. ELSE
  317. C
  318. C******* NGCEG != NGCED
  319. C
  320. INDCEL = (NGCEG-1)*IDIMP1
  321. XG = XCOOR(INDCEL+1)
  322. YG = XCOOR(INDCEL+2)
  323. ZG = XCOOR(INDCEL+3)
  324. INDCEL = (NGCED-1)*IDIMP1
  325. XD = XCOOR(INDCEL+1)
  326. YD = XCOOR(INDCEL+2)
  327. ZD = XCOOR(INDCEL+3)
  328. INDCEL = (NGCF-1)*IDIMP1
  329. XC = XCOOR(INDCEL+1)
  330. YC = XCOOR(INDCEL+2)
  331. ZC = XCOOR(INDCEL+3)
  332. C
  333. DXG = XC - XG
  334. DYG = YC - YG
  335. DZG = ZC - ZG
  336. DXD = XC - XD
  337. DYD = YC - YD
  338. DZD = ZC - ZD
  339. DG = DXG * DXG + DYG * DYG + DZG * DZG
  340. DG = SQRT(DG)
  341. DD = DXD * DXD + DYD * DYD + DZD * DZD
  342. DD = SQRT(DD)
  343. DT = DG + DD
  344. C
  345. C********** Boucle sur le composantes
  346. C
  347. DO I1 = 1, NCOMP
  348. IGR = IDIM*(I1-1)+1
  349. C
  350. C************* Limiteur a gauche
  351. C
  352. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  353. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  354. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  355. DPHIST = PHIMAX - PHIMIN
  356. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  357. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  358. GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
  359. DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
  360. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  361. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  362. VALEUR = 1.0D0
  363. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  364. DPHI1 = PHIMAX - PHI
  365. VALEUR = DPHI1/DPHI0*DG/DT
  366. ELSE
  367. DPHI1 = PHIMIN - PHI
  368. VALEUR = DPHI1/DPHI0*DG/DT
  369. ENDIF
  370.  
  371. ALPHA = MIN(ALPHA, VALEUR)
  372. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  373. C
  374. C************* Limiteur a droite
  375. C
  376. PHI = MPOCHP.VPOCHA(NLCED,I1)
  377. PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
  378. PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
  379. DPHIST = PHIMAX - PHIMIN
  380. GRADX = MPOGRA.VPOCHA(NLCED,IGR)
  381. GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
  382. GRADZ = MPOGRA.VPOCHA(NLCED,IGR+2)
  383. DPHI0 = GRADX * DXD + GRADY * DYD + GRADZ * DZD
  384. ALPHA = MPOALP.VPOCHA(NLCED,I1)
  385. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  386. VALEUR = 1.0D0
  387. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  388. DPHI1 = PHIMAX - PHI
  389. VALEUR = DPHI1/DPHI0*DD/DT
  390. ELSE
  391. DPHI1 = PHIMIN - PHI
  392. VALEUR = DPHI1/DPHI0*DD/DT
  393. ENDIF
  394. ALPHA = MIN(ALPHA, VALEUR)
  395. MPOALP.VPOCHA(NLCED,I1) = ALPHA
  396. ENDDO
  397. ENDIF
  398. ENDDO
  399. ENDIF
  400. C
  401. RETURN
  402. END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  

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