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 PPARAM
  64. -INC CCOPTIO
  65. C
  66. C**** Variables de CCOPTIO
  67. C
  68. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  69. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  70. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  71. C & ,IECHO, IIMPI, IOSPI
  72. C & ,IDIM
  73. CC & ,MCOORD
  74. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  75. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  76. C & ,NORINC,NORVAL,NORIND,NORVAD
  77. C & ,NUCROU, IPSAUV
  78. C
  79. INTEGER NFAC,NCOMP
  80. & ,NLCF,NGCF,NGCEG,NLCEG,NGCED,NLCED
  81. & ,INDCEL,I1,IGR,IDIMP1
  82. C
  83. C**** Pour l'instant 2D
  84. C
  85. REAL*8 XG,YG,XD,YD,XC,YC,DXG,DYG,DXD,DYD,DG,DD,DT
  86. & ,PHI,PHIMAX,PHIMIN,DPHIST,GRADX,GRADY,DPHI0,DPHI1,ALPHA
  87. & ,VALEUR
  88. C
  89. C**** Extension 3D
  90. C
  91. REAL*8 ZG,ZD,ZC,DZG,DZD,GRADZ
  92. C
  93. C
  94. -INC SMCOORD
  95. -INC SMCHPOI
  96. -INC SMELEME
  97. -INC SMLENTI
  98. POINTEUR MPOMIN.MPOVAL, MPOMAX.MPOVAL, MPOALP.MPOVAL,
  99. & MPOCHP.MPOVAL, MPOGRA.MPOVAL
  100. POINTEUR MELEFL.MELEME, MLECEN.MLENTI
  101. C
  102. C**** N.B. Tous les pointeurs ici sont déjà activés!
  103. C
  104. C**** Maillage FACEL
  105. C Boucle sur les faces
  106. C
  107. C
  108. C**** NGCEG = Numero global centre d'elt "gauche"
  109. C NLCEG = Numero local centre d'elt "gauche"
  110. C NGCF = " global centre de face
  111. C NLCF = numero local centre de face
  112. C NGCEG = Numero global centre d'elt "gauche"
  113. C NLCEG = Numero local centre d'elt "gauche"
  114. C
  115. IDIMP1 = IDIM+1
  116. C
  117. C Cas de la dimension 2
  118. C
  119. IF (IDIM .EQ. 2) THEN
  120. DO NLCF = 1, NFAC
  121. NGCEG = MELEFL.NUM(1,NLCF)
  122. NLCEG = MLECEN.LECT(NGCEG)
  123. NGCF = MELEFL.NUM(2,NLCF)
  124. NGCED = MELEFL.NUM(3,NLCF)
  125. NLCED = MLECEN.LECT(NGCED)
  126. IF(NGCEG .EQ. NGCED)THEN
  127. C
  128. C********** Limiteur dans le cas mur
  129. C
  130. C
  131. C********** Coordonees et parametres geometriques
  132. C
  133. C
  134. C XCOOR : table de coordonnes de points;
  135. C pour le point de numero global NG
  136. C X_NG = XCOOR((NG-1)*(IDIM+1)+1)
  137. C Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
  138. C Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
  139. C (dans notre cas IDIM=2, i.e. 2D)
  140. C
  141. INDCEL = (NGCEG-1)*IDIMP1
  142. XG = XCOOR(INDCEL+1)
  143. YG = XCOOR(INDCEL+2)
  144. INDCEL = (NGCF-1)*IDIMP1
  145. XC = XCOOR(INDCEL+1)
  146. YC = XCOOR(INDCEL+2)
  147. DXG = XC - XG
  148. DYG = YC - YG
  149. C
  150. C********** Boucle sur le composantes
  151. C
  152. DO I1 = 1, NCOMP
  153. IGR = 2*I1 - 1
  154. C
  155. C************* Limiteur
  156. C
  157. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  158. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  159. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  160. DPHIST = PHIMAX - PHIMIN
  161. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  162. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  163. DPHI0 = GRADX * DXG + GRADY * DYG
  164. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  165. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  166. VALEUR = 1.0D0
  167. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  168. DPHI1 = PHIMAX - PHI
  169. VALEUR = DPHI1/DPHI0*0.5D0
  170. ELSE
  171. DPHI1 = PHIMIN - PHI
  172. VALEUR = DPHI1/DPHI0*0.5D0
  173. ENDIF
  174. ALPHA = MIN(ALPHA, VALEUR)
  175. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  176. ENDDO
  177. ELSE
  178. C
  179. C******* NGCEG != NGCED
  180. C
  181. INDCEL = (NGCEG-1)*IDIMP1
  182. XG = XCOOR(INDCEL+1)
  183. YG = XCOOR(INDCEL+2)
  184. INDCEL = (NGCED-1)*IDIMP1
  185. XD = XCOOR(INDCEL+1)
  186. YD = XCOOR(INDCEL+2)
  187. INDCEL = (NGCF-1)*IDIMP1
  188. XC = XCOOR(INDCEL+1)
  189. YC = XCOOR(INDCEL+2)
  190. C
  191. DXG = XC - XG
  192. DYG = YC - YG
  193. DXD = XC - XD
  194. DYD = YC - YD
  195. DG = DXG * DXG + DYG * DYG
  196. DG = SQRT(DG)
  197. DD = DXD * DXD + DYD * DYD
  198. DD = SQRT(DD)
  199. DT = DG + DD
  200. C
  201. C********** Boucle sur le composantes
  202. C
  203. DO I1 = 1, NCOMP
  204. IGR = 2*I1 - 1
  205. C
  206. C************* Limiteur a gauche
  207. C
  208. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  209. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  210. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  211. DPHIST = PHIMAX - PHIMIN
  212. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  213. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  214. DPHI0 = GRADX * DXG + GRADY * DYG
  215. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  216. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  217. VALEUR = 1.0D0
  218. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  219. DPHI1 = PHIMAX - PHI
  220. VALEUR = DPHI1/DPHI0*DG/DT
  221. ELSE
  222. DPHI1 = PHIMIN - PHI
  223. VALEUR = DPHI1/DPHI0*DG/DT
  224. ENDIF
  225.  
  226. ALPHA = MIN(ALPHA, VALEUR)
  227. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  228. C
  229. C************* Limiteur a droite
  230. C
  231. PHI = MPOCHP.VPOCHA(NLCED,I1)
  232. PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
  233. PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
  234. DPHIST = PHIMAX - PHIMIN
  235. GRADX = MPOGRA.VPOCHA(NLCED,IGR)
  236. GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
  237. DPHI0 = GRADX * DXD + GRADY * DYD
  238. ALPHA = MPOALP.VPOCHA(NLCED,I1)
  239. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  240. VALEUR = 1.0D0
  241. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  242. DPHI1 = PHIMAX - PHI
  243. VALEUR = DPHI1/DPHI0*DD/DT
  244. ELSE
  245. DPHI1 = PHIMIN - PHI
  246. VALEUR = DPHI1/DPHI0*DD/DT
  247. ENDIF
  248. ALPHA = MIN(ALPHA, VALEUR)
  249. MPOALP.VPOCHA(NLCED,I1) = ALPHA
  250. ENDDO
  251. ENDIF
  252. ENDDO
  253. C
  254. C Cas de la dimension 3
  255. C
  256. ELSE
  257. DO NLCF = 1, NFAC
  258. NGCEG = MELEFL.NUM(1,NLCF)
  259. NLCEG = MLECEN.LECT(NGCEG)
  260. NGCF = MELEFL.NUM(2,NLCF)
  261. NGCED = MELEFL.NUM(3,NLCF)
  262. NLCED = MLECEN.LECT(NGCED)
  263. IF(NGCEG .EQ. NGCED)THEN
  264. C
  265. C********** Limiteur dans le cas mur
  266. C
  267. C
  268. C********** Coordonees et parametres geometriques
  269. C
  270. C
  271. C XCOOR : table de coordonnes de points;
  272. C pour le point de numero global NG
  273. C X_NG = XCOOR((NG-1)*(IDIM+1)+1)
  274. C Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
  275. C Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
  276. C
  277. INDCEL = (NGCEG-1)*IDIMP1
  278. XG = XCOOR(INDCEL+1)
  279. YG = XCOOR(INDCEL+2)
  280. ZG = XCOOR(INDCEL+3)
  281. INDCEL = (NGCF-1)*IDIMP1
  282. XC = XCOOR(INDCEL+1)
  283. YC = XCOOR(INDCEL+2)
  284. ZC = XCOOR(INDCEL+3)
  285. DXG = XC - XG
  286. DYG = YC - YG
  287. DZG = ZC - ZG
  288. C
  289. C********** Boucle sur le composantes
  290. C
  291. DO I1 = 1, NCOMP
  292. IGR = IDIM*(I1-1)+1
  293. C
  294. C************* Limiteur
  295. C
  296. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  297. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  298. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  299. DPHIST = PHIMAX - PHIMIN
  300. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  301. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  302. GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
  303. DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
  304. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  305. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  306. VALEUR = 1.0D0
  307. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  308. DPHI1 = PHIMAX - PHI
  309. VALEUR = DPHI1/DPHI0*0.5D0
  310. ELSE
  311. DPHI1 = PHIMIN - PHI
  312. VALEUR = DPHI1/DPHI0*0.5D0
  313. ENDIF
  314. ALPHA = MIN(ALPHA, VALEUR)
  315. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  316. ENDDO
  317. ELSE
  318. C
  319. C******* NGCEG != NGCED
  320. C
  321. INDCEL = (NGCEG-1)*IDIMP1
  322. XG = XCOOR(INDCEL+1)
  323. YG = XCOOR(INDCEL+2)
  324. ZG = XCOOR(INDCEL+3)
  325. INDCEL = (NGCED-1)*IDIMP1
  326. XD = XCOOR(INDCEL+1)
  327. YD = XCOOR(INDCEL+2)
  328. ZD = XCOOR(INDCEL+3)
  329. INDCEL = (NGCF-1)*IDIMP1
  330. XC = XCOOR(INDCEL+1)
  331. YC = XCOOR(INDCEL+2)
  332. ZC = XCOOR(INDCEL+3)
  333. C
  334. DXG = XC - XG
  335. DYG = YC - YG
  336. DZG = ZC - ZG
  337. DXD = XC - XD
  338. DYD = YC - YD
  339. DZD = ZC - ZD
  340. DG = DXG * DXG + DYG * DYG + DZG * DZG
  341. DG = SQRT(DG)
  342. DD = DXD * DXD + DYD * DYD + DZD * DZD
  343. DD = SQRT(DD)
  344. DT = DG + DD
  345. C
  346. C********** Boucle sur le composantes
  347. C
  348. DO I1 = 1, NCOMP
  349. IGR = IDIM*(I1-1)+1
  350. C
  351. C************* Limiteur a gauche
  352. C
  353. PHI = MPOCHP.VPOCHA(NLCEG,I1)
  354. PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
  355. PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
  356. DPHIST = PHIMAX - PHIMIN
  357. GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
  358. GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
  359. GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
  360. DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
  361. ALPHA = MPOALP.VPOCHA(NLCEG,I1)
  362. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  363. VALEUR = 1.0D0
  364. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  365. DPHI1 = PHIMAX - PHI
  366. VALEUR = DPHI1/DPHI0*DG/DT
  367. ELSE
  368. DPHI1 = PHIMIN - PHI
  369. VALEUR = DPHI1/DPHI0*DG/DT
  370. ENDIF
  371.  
  372. ALPHA = MIN(ALPHA, VALEUR)
  373. MPOALP.VPOCHA(NLCEG,I1) = ALPHA
  374. C
  375. C************* Limiteur a droite
  376. C
  377. PHI = MPOCHP.VPOCHA(NLCED,I1)
  378. PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
  379. PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
  380. DPHIST = PHIMAX - PHIMIN
  381. GRADX = MPOGRA.VPOCHA(NLCED,IGR)
  382. GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
  383. GRADZ = MPOGRA.VPOCHA(NLCED,IGR+2)
  384. DPHI0 = GRADX * DXD + GRADY * DYD + GRADZ * DZD
  385. ALPHA = MPOALP.VPOCHA(NLCED,I1)
  386. IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
  387. VALEUR = 1.0D0
  388. ELSEIF(DPHI0 .GT. 0.0D0)THEN
  389. DPHI1 = PHIMAX - PHI
  390. VALEUR = DPHI1/DPHI0*DD/DT
  391. ELSE
  392. DPHI1 = PHIMIN - PHI
  393. VALEUR = DPHI1/DPHI0*DD/DT
  394. ENDIF
  395. ALPHA = MIN(ALPHA, VALEUR)
  396. MPOALP.VPOCHA(NLCED,I1) = ALPHA
  397. ENDDO
  398. ENDIF
  399. ENDDO
  400. ENDIF
  401. C
  402. RETURN
  403. END
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  

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