Télécharger pente2.eso

Retour à la liste

Numérotation des lignes :

  1. C PENTE2 SOURCE KK2000 14/04/10 21:15:27 8032
  2. SUBROUTINE PENTE2(IOP2,INORM,MLECEN,MLEFAC,
  3. & MLENCL,IMCHAM,
  4. & NCOMP,MPOCHP,MPOVCL,MPOGRA,
  5. & MPOMIN,MPOMAX)
  6. C
  7. C************************************************************************
  8. C
  9. C PROJET : CASTEM 2000
  10. C
  11. C NOM : PENTE2
  12. C
  13. C DESCRIPTION : Cette subroutine est appellée par la subroutine
  14. C PENTE1 (calcul du gradient d'un CHPOINT de type
  15. C CENTRE)
  16. C Elle contient la partie du calcul du gradient.
  17. C
  18. C LANGUAGE : ESOPE 2000 avec extensions CISI
  19. C
  20. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
  21. C AUTEUR (modif.) : R. MOREL, DRN/DMT/SEMT/TTMF
  22. C
  23. C************************************************************************
  24. C
  25. C
  26. C APPELES (Outils) : LICHT, ERREUR
  27. C
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES : IOP2 : Entier indiquant le type de reconstruction (cf PENT)
  32. C
  33. C INORM : Pointeur sur le CHPOINT des normales aux faces
  34. C
  35. C MLECEN : table numerotation global/local pour les
  36. C points CENTRE (segment MLENTI)
  37. C
  38. C MLEFAC : table numerotation global/local pour les
  39. C points FACE
  40. C
  41. C MLENCL : table numerotation global/local pour les
  42. C SPG du CHPOINT de C.L. ou
  43. C table de 0 si ce CHPOINT n'existe pas
  44. C
  45. C IMCHAM : MCHAML de type 'GRADGEO'
  46. C
  47. C NCOMP : Nombre de composantes du CHPOINT dont
  48. C on veut le gradient
  49. C
  50. C MPOCHP : MPOVAL du CHPOINT dont on veut le gradient
  51. C
  52. C MPOVCL : MPOVAL du CHPOINT de C.L.
  53. C
  54. C SORTIES : MPOGRA : MPOVAL du gradient
  55. C
  56. C MPOMAX : MPOVAL du max sur le stencil
  57. C
  58. C MPOMIN : MPOVAL du min sur le stencil
  59. C
  60. C
  61. C************************************************************************
  62. C
  63. C HISTORIQUE (Anomalies et modifications éventuelles)
  64. C
  65. C HISTORIQUE : Cree le 4-6-1998
  66. C
  67. C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D
  68. C Modifie le 28-04-2000 pour reconstruction quadatique
  69. C exacte (A. Beccantini)
  70. C HISTORIQUE : Modifie le 05.07.2003: elimination reconstruction
  71. C quadratique exacte.
  72. C Algorithme de BArth-Jespersen remplacé par une
  73. C reconstruction lineaire exacte qui se base sur
  74. C LSM
  75. C************************************************************************
  76. C
  77. IMPLICIT INTEGER(I-N)
  78. IMPLICIT REAL*8(A-H,O-Z)
  79.  
  80. -INC CCOPTIO
  81. -INC SMCHAML
  82. -INC SMELEME
  83. -INC SMCHPOI
  84. -INC SMLENTI
  85. C
  86. POINTEUR INORM.MCHPOI
  87. POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLEFAC.MLENTI
  88. POINTEUR MPOMAX.MPOVAL, MPOMIN.MPOVAL, MPONOR.MPOVAL,
  89. & MPOVCL.MPOVAL, MPOCHP.MPOVAL, MPOGRA.MPOVAL
  90. POINTEUR MELVX.MELVAL, MELVY.MELVAL, MELVZ.MELVAL,
  91. & MELVXX.MELVAL, MELVYY.MELVAL, MELVZZ.MELVAL,
  92. & MELVXY.MELVAL, MELVXZ.MELVAL, MELVYZ.MELVAL
  93. C
  94. INTEGER IOP2,NCOMP,NLELEM,NVOI,I3,IMCHAM
  95. & ,NGCE,NLCE,NTOTV,NGCV,NLCV,ICOM,IGEOM,K,NBELEM,NBNN
  96. & ,NLCF,NLCL,NMAI
  97. REAL*8 VAL,VCOMP(3),VN
  98. REAL*8 NX,NY,NZ
  99. CHARACTER*8 TYPE
  100. C
  101. C**** Segments déjà activés
  102. C
  103. C MPOVCL
  104. C MLENCL
  105. C MLECEN
  106. C MLEFAC
  107. C MPOCHP
  108. C MPOGRA*MOD
  109. C MPOMIN*MOD
  110. C MPOMAX*MOD
  111. C
  112. C**** N.B. MPOMAX.VPOCHA, MPOMIN.VPOCHA déjà initialisé:
  113. C
  114. C MPOMAX.VPOCHA = MPOCHP.VPOCHA
  115. C MPOMIN.VPOCHA = MPOCHP.VPOCHA
  116. C
  117. CALL LICHT(INORM,MPONOR,TYPE,IGEOM)
  118. C
  119. C**** En LICHT
  120. C SEGACT*MOD MPONOR
  121. C MPONOR: pointeur sur les normales aux faces
  122. C
  123. MCHELM = IMCHAM
  124. SEGACT MCHELM
  125. NMAI=MCHELM.ICHAML(/1)
  126. C
  127. C*********************************************
  128. C Premier cas : on traite un champ scalaire en
  129. C reconstruction non-quadratique
  130. C*********************************************
  131. C
  132. C
  133. IF(IOP2.EQ.1) THEN
  134. C
  135. C**** N.B: IOP2 = 1 -> EULESCAL
  136. C
  137. DO 10 K=1,NMAI
  138. C
  139. C********** Boucle sur les differents maillages elementaires
  140. C IPT1 = maillage duale
  141. C
  142. IPT1=MCHELM.IMACHE(K)
  143. SEGACT IPT1
  144. NBNN=IPT1.NUM(/1)
  145. NTOTV=NBNN
  146. C
  147. C********** NTOTV = Nombre de voisins
  148. C
  149. NBELEM = IPT1.NUM(/2)
  150. MCHAM1 = MCHELM.ICHAML(K)
  151. SEGACT MCHAM1
  152. MELVX = MCHAM1.IELVAL(1)
  153. MELVY = MCHAM1.IELVAL(2)
  154. SEGACT MELVX
  155. SEGACT MELVY
  156. IF (IDIM .EQ. 3) THEN
  157. MELVZ = MCHAM1.IELVAL(3)
  158. SEGACT MELVZ
  159. ENDIF
  160. SEGDES MCHAM1
  161. DO NLELEM = 1, NBELEM, 1
  162. C
  163. C************* Boucle sur les ELTs
  164. C NLELEM = Numero local d'ELT dans le maillage 1
  165. C
  166. C NGCE = numero global centre d'elt
  167. C NLCE = numero local centre d'elt
  168. C
  169. C MLECEN = table numerotation global/local pour les
  170. C points CENTRE (segment MLENTI)
  171. C
  172. C MLEFAC = table numerotation global/local pour les
  173. C points FACE
  174. C
  175. C MLENCL = table numerotation global/local pour les
  176. C SPG du CHPOINT de C.L. ou
  177. C table de 0 si ce CHPOINT n'existe pas
  178. C
  179. NGCE = IPT1.NUM(NBNN,NLELEM)
  180. NLCE = MLECEN.LECT(NGCE)
  181. DO NVOI = 1, NTOTV
  182. C
  183. C**************** Boucle sur les voisins de NGCE
  184. C
  185. C NGCV = Numero global voisin
  186. C NLCV = " local de NGCV dans le maillage CENTRE
  187. C NLCL = " local de NGCV dans le SPG du CHPOINT
  188. C des C.L.
  189. C
  190. NGCV = IPT1.NUM(NVOI,NLELEM)
  191. NLCV = MLECEN.LECT(NGCV)
  192. IF(NLCV.EQ.0)THEN
  193. NLCL = MLENCL.LECT(NGCV)
  194. C
  195. C**************** NGCV n'est pas un point centre, mais un point face,
  196. C et on est sur les bords
  197. C
  198. DO I3 = 1, NCOMP
  199. C
  200. C********************** I3 : numero de composante du CHPOI
  201. C ICOM : numero de composante du gradient
  202. C
  203. ICOM = IDIM*(I3 -1)+1
  204. IF (NLCL .EQ. 0) THEN
  205. C
  206. C************************* NGCV n'est pas un point du SPG de C.L.
  207. C
  208. VAL = MPOCHP.VPOCHA(NLCE,I3)
  209. ELSE
  210. C
  211. C************************* NGCV est un point du SPG de C.L.
  212. C
  213. VAL = MPOVCL.VPOCHA(NLCL,I3)
  214. ENDIF
  215. MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE
  216. $ ,ICOM)+VAL * MELVX.VELCHE(NVOI,NLELEM)
  217. MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  218. $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI,NLELEM)
  219. IF (IDIM .EQ. 3) THEN
  220. MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  221. $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI
  222. $ ,NLELEM)
  223. ENDIF
  224. MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  225. $ ),VAL)
  226. MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  227. $ ),VAL)
  228. ENDDO
  229. ELSE
  230. C
  231. C********** NGCV est un point centre
  232. C
  233. DO I3 = 1, NCOMP
  234. ICOM = IDIM*(I3 -1)+1
  235. VAL = MPOCHP.VPOCHA(NLCV,I3)
  236. MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM
  237. $ )+VAL * MELVX.VELCHE(NVOI,NLELEM)
  238. MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  239. $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM)
  240. IF (IDIM .EQ. 3) THEN
  241. MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  242. $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM
  243. $ )
  244. ENDIF
  245. MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  246. $ ),VAL)
  247. MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  248. $ ),VAL)
  249. ENDDO
  250. ENDIF
  251. ENDDO
  252. C
  253. C******* Fin boucle sur les voisins de NGCE
  254. C
  255. ENDDO
  256. SEGDES MELVX
  257. SEGDES MELVY
  258. IF (IDIM .EQ. 3) SEGDES MELVZ
  259. SEGDES IPT1
  260. 10 CONTINUE
  261. C
  262. C**********************************************************
  263. C Deuxieme cas : on traite un champ vectoriel (IOP2=2)
  264. C**********************************************************
  265. C
  266. C
  267. ELSEIF(IOP2 .EQ. 2)THEN
  268. DO 20 K=1,NMAI
  269. IPT1=MCHELM.IMACHE(K)
  270. SEGACT IPT1
  271. NBNN=IPT1.NUM(/1)
  272. NTOTV=NBNN
  273. NBELEM=IPT1.NUM(/2)
  274. MCHAM1=MCHELM.ICHAML(K)
  275. SEGACT MCHAM1
  276. MELVX=MCHAM1.IELVAL(1)
  277. MELVY=MCHAM1.IELVAL(2)
  278. SEGACT MELVX
  279. SEGACT MELVY
  280. IF (IDIM .EQ. 3) THEN
  281. MELVZ=MCHAM1.IELVAL(3)
  282. SEGACT MELVZ
  283. ENDIF
  284. SEGDES MCHAM1
  285. DO NLELEM = 1, NBELEM
  286. NGCE = IPT1.NUM(NBNN,NLELEM)
  287. NLCE = MLECEN.LECT(NGCE)
  288. DO NVOI = 1, NTOTV
  289. NGCV = IPT1.NUM(NVOI,NLELEM)
  290. NLCV = MLECEN.LECT(NGCV)
  291. IF(NLCV .EQ. 0)THEN
  292. NLCL = MLENCL.LECT(NGCV)
  293. C
  294. C************* NGCV n'est pas un point centre, mais un point face.
  295. C
  296. IF(NLCL .NE. 0)THEN
  297. C
  298. C************************* NGCV est un point du SPG de C.L.
  299. C
  300. DO I3 = 1, NCOMP
  301. VCOMP(I3) = MPOVCL.VPOCHA(NLCL,I3)
  302. ENDDO
  303. ELSE
  304. NLCF = MLEFAC.LECT(NGCV)
  305. IF (IDIM .EQ.2) THEN
  306. VCOMP(1) = MPOCHP.VPOCHA(NLCE,1)
  307. VCOMP(2) = MPOCHP.VPOCHA(NLCE,2)
  308. NX = MPONOR.VPOCHA(NLCF,1)
  309. NY = MPONOR.VPOCHA(NLCF,2)
  310. VN = VCOMP(1)*NX+VCOMP(2)*NY
  311. VCOMP(1) = VCOMP(1)-2*NX*VN
  312. VCOMP(2) = VCOMP(2)-2*NY*VN
  313. ELSE
  314. VCOMP(1) = MPOCHP.VPOCHA(NLCE,1)
  315. VCOMP(2) = MPOCHP.VPOCHA(NLCE,2)
  316. VCOMP(3) = MPOCHP.VPOCHA(NLCE,3)
  317. NX = MPONOR.VPOCHA(NLCF,1)
  318. NY = MPONOR.VPOCHA(NLCF,2)
  319. NZ = MPONOR.VPOCHA(NLCF,3)
  320. VN = VCOMP(1)*NX+VCOMP(2)*NY+VCOMP(3)*NZ
  321. VCOMP(1) = VCOMP(1)-2*NX*VN
  322. VCOMP(2) = VCOMP(2)-2*NY*VN
  323. VCOMP(3) = VCOMP(3)-2*NZ*VN
  324. ENDIF
  325. ENDIF
  326. DO I3 = 1, NCOMP
  327. ICOM = IDIM*(I3 -1)+1
  328. VAL = VCOMP(I3)
  329. MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE
  330. $ ,ICOM)+VAL * MELVX.VELCHE(NVOI
  331. $ ,NLELEM)
  332. MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  333. $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI
  334. $ ,NLELEM)
  335. IF (IDIM .EQ. 3) THEN
  336. MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  337. $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI
  338. $ ,NLELEM)
  339. ENDIF
  340. MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  341. $ ),VAL)
  342. MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  343. $ ),VAL)
  344. ENDDO
  345. ELSE
  346. C
  347. C********** NGCV est un point centre
  348. C
  349. DO I3 = 1, NCOMP
  350. ICOM = IDIM*(I3 -1)+1
  351. VAL = MPOCHP.VPOCHA(NLCV,I3)
  352. MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM
  353. $ )+VAL * MELVX.VELCHE(NVOI,NLELEM)
  354. MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  355. $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM)
  356. IF (IDIM .EQ. 3) THEN
  357. MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  358. $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM
  359. $ )
  360. ENDIF
  361. MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  362. $ ),VAL)
  363. MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  364. $ ),VAL)
  365. ENDDO
  366. ENDIF
  367. ENDDO
  368. C
  369. C******* Fin boucle sur les voisins de NGCE
  370. C
  371. ENDDO
  372. SEGDES MELVX
  373. SEGDES MELVY
  374. IF (IDIM .EQ. 3) SEGDES MELVZ
  375. SEGDES IPT1
  376. 20 CONTINUE
  377. ENDIF
  378. SEGDES MPONOR
  379. SEGDES MCHELM
  380. RETURN
  381. END
  382.  
  383.  
  384.  
  385.  
  386.  

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