Télécharger fptau.eso

Retour à la liste

Numérotation des lignes :

  1. C FPTAU SOURCE PV 08/09/11 21:15:54 6150
  2. SUBROUTINE FPTAU(IDIM,NP,NBEL,NUM,IPADT,XCOOR,
  3. & U,NPTI,NC,IPADS,NPTS,ANUT,AMUEF,
  4. & YP,AN,RO,IKR,AMU,IKM,UET,YPS,ITYPEL,TKE,NBINC)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8 (A-H,O-Z)
  7. C.......................................................................
  8. C Ut(NP,NBEL) module vitesse tangente = |U - Un|
  9. C Un(NP,NBEL) module vitesse normale = Un
  10. C YP(NP) distance la paroi
  11. C UET(*) uet aux noeuds frontieres
  12. C AN(IDIM,NP) normale
  13. C ANN module de la normale
  14. C NUM(NP) connectivite element appuyé
  15. C NMS(NS) connectivite relative (frontiere)
  16. C NMY(NY) connectivite complementaire
  17. C RO(*) densite
  18. C AMUEF(NP,NBEL) Viscosite effective aux noeuds/elt
  19. C
  20. C
  21. C1/ ce qu'il faudrait faire c'est attribuer une normale pour chaque point
  22. C de l'element appuye -> c'est fait
  23. C.......................................................................
  24.  
  25. -INC CCREEL
  26. DIMENSION LKKE(1000)
  27. real*8 ann(3),Xa(9),Ya(9),Za(9)
  28. DIMENSION XCOOR(*),XYZ(3,27)
  29. DIMENSION NUM(NP,NBEL),NMS(9),NMY(27)
  30. DIMENSION U(NPTI,NC),IPADT(*),VS(27,3),IPADS(*)
  31. DIMENSION ANUT(NPTI),AMUEF(NP,NBEL)
  32. DIMENSION YP(NP,NBEL),AN(IDIM,NP,NBEL)
  33. DIMENSION UT(27),UN(27)
  34. DIMENSION RO(*),AMU(*)
  35. DIMENSION UET(NPTS),YPS(NPTS),TKE(NPTS,*)
  36. DIMENSION NUMFA(27)
  37. DATA EPS/1.D-4/,NBIT/100/
  38. c-----------------------------------------------------------------------
  39. c------- Caracterisation des elements ----------------------------------
  40. c-----------------------------------------------------------------------
  41. DIMENSION NMF9(3,4),N1F9(3,4),N2F9(3,4)
  42. DATA NMF9/1,2,3,3,4,5,5,6,7,7,8,1/
  43. DATA N1F9/8,9,4,2,9,6,4,9,8,6,9,2/,N2F9/7,6,5,1,8,7,3,2,1,5,4,3/
  44. c-----------------------------------------------------------------------
  45. DIMENSION NMF7(3,3),N1F7(3,3),N2F7(1,3)
  46. DATA NMF7/1,2,3,3,4,5,5,6,1/
  47. DATA N1F7/6,7,4,2,7,6,4,7,2/,N2F7/5,1,3/
  48. c-----------------------------------------------------------------------
  49. DIMENSION NMF8(4,6),N1F8(4,6),N2F8(4,6)
  50. DATA NMF8/1,2,3,4,2,3,7,6,5,6,7,8,1,4,8,5,1,2,6,5,4,3,7,8/
  51. c? DATA N1F8/6,7,4,2,7,6,4,7,2/,N2F8/5,1,3/
  52. c-----------------------------------------------------------------------
  53. C CONSTANTES PHYSIQUES
  54. DATA CNU,AKER /0.09D0,0.41D0/
  55. c-----------------------------------------------------------------------
  56. c-----------------------------------------------------------------------
  57.  
  58. SQCNU=SQRT(CNU)
  59.  
  60. c write(6,*)' DEBUT FPTAU : NBINC=',NBINC
  61. c write(6,*)'NES=',nes,' NPG=',npg,' idim=',idim,' np=',np
  62. c write(6,*)'NPTS=',NPTS
  63.  
  64. DEUPI=1.D0
  65. cc IF(IAXI.NE.0)DEUPI=2.D0*XPI
  66. IERC=0
  67. ERRMAX=0.D0
  68.  
  69. KT=0
  70. DO 2 K=1,NBEL
  71. c write(6,*)' DEBUT BCL sur K ',K
  72. NS=0
  73. NY=0
  74. CALL INITI(NMS,NP,0)
  75. c CALL INITI(MMS,NP,0)
  76. CALL INITI(NMY,NP,0)
  77. c write(6,*)' ELEMENT K=',K
  78.  
  79. DO 20 I=1,NP
  80. J=NUM(I,K)
  81. JI=IPADT(NUM(I,K))
  82. JM=(1-IKM)*(JI-1)+1
  83. c JS=IPADS(NUM(I,K))
  84. c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=',
  85. c & IPADS(J)
  86. IF(IPADS(J).NE.0)THEN
  87. NS=NS+1
  88. NMS(NS)=I
  89. c MMS(NS)=I
  90. c write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=',
  91. c & IPADS(J),' NS=',NS,' ',NMS(NS)
  92. ELSE
  93. NY=NY+1
  94. NMY(NY)=I
  95. ENDIF
  96. DO 10 N=1,IDIM
  97. XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N)
  98. VS(I,N)=U(JI,N)
  99. 10 CONTINUE
  100. c UETK(I,K)=VET(JS)
  101. c UETK(I,K)=0.D0
  102. c ROS(I) = RO(JR)
  103. c AMUS(I)= AMU(JM)
  104. 20 CONTINUE
  105. c write(6,*)'K=',K,' MMS=',(mms(ijk),ijk=1,nS)
  106.  
  107. c write(6,*)'K=',K,' XYZ=',(XYZ(1,I),I=1,np)
  108. c write(6,*)'K=',K,' XYZ=',(XYZ(2,I),I=1,np)
  109. c write(6,*)'K=',K,' NP=',NP,' NS=',NS,' NY=',NY
  110. c write(6,*)'K=',K,' NUM ',(num(i,k),i=1,np)
  111. c write(6,*)'K=',K,' NMSO',(nms(i),i=1,nS)
  112. c write(6,*)'K=',K,' NMS ',(num(nms(i),k),i=1,nS)
  113. c write(6,*)'K=',K,' NMY ',(num(NMY(I),k),I=1,NY)
  114.  
  115. IF(NY.EQ.1)GO TO 2
  116. IF(NS.EQ.1)GO TO 2
  117. IF(NS.LE.2.AND.IDIM.EQ.3)GO TO 2
  118. KT=KT+1
  119. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  120. c._._ On Calcule la normale au SEG2 (NMS) et YP
  121. c._._ (La normale est la meme pour les deux pts)
  122. IF(IDIM.EQ.2.AND.NS.EQ.2)THEN
  123.  
  124. Xa(1)=XYZ(1,NMS(1))
  125. Ya(1)=XYZ(2,NMS(1))
  126. Xa(2)=XYZ(1,NMS(2))
  127. Ya(2)=XYZ(2,NMS(2))
  128. c a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2))
  129. c a=1.
  130. aa=Ya(1)-Ya(2)
  131. bb=Xa(1)-Xa(2)
  132. AN(1,NMS(1),K)=aa
  133. AN(2,NMS(1),K)=(-1.)*bb
  134. ann(1)=((aa*aa)+(bb*bb))**0.5D0
  135. AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1)
  136. AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1)
  137. DO 241 I=1,NY
  138. Xq=XYZ(1,NMY(I))
  139. Yq=XYZ(2,NMY(I))
  140. YP(NMY(I),K)=abs((aa*Xq)-(bb*Yq)+(bb*Ya(1))-(aa*Xa(1)))/ann(1)
  141. 241 CONTINUE
  142. c la normale est la meme pour tous les pts (SEG2)
  143. DO 251 I=2,NS
  144. AN(1,NMS(I),K)=AN(1,NMS(1),K)
  145. AN(2,NMS(I),K)=AN(2,NMS(1),K)
  146. 251 CONTINUE
  147. DO 252 I=1,NY
  148. AN(1,NMY(I),K)=AN(1,NMS(1),K)
  149. AN(2,NMY(I),K)=AN(2,NMS(1),K)
  150. 252 CONTINUE
  151. c._._ FIN On Calcule la normale au SEG2 (NMS) et YP
  152. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  153.  
  154. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  155. c._._ On Calcule la normale au SEG3 (NMS) et YP
  156. c._._ (La normale n'est pas la meme pour tous les pts)
  157. ELSEIF(IDIM.EQ.2.AND.NS.EQ.3)THEN
  158.  
  159. C._._ Mais quelle est donc la face?
  160. call initi(numfa,np,0)
  161. do 2300 I=1,NS
  162. numfa(nms(I))=1
  163. 2300 continue
  164. nfa=numfa(2)+numfa(4)*2+numfa(6)*3+numfa(8)*4
  165. c write(6,*)'NFA=',NFA
  166.  
  167. Xa(1)=XYZ(1,NMF9(1,NFA))
  168. Ya(1)=XYZ(2,NMF9(1,NFA))
  169. Xa(2)=XYZ(1,NMF9(2,NFA))
  170. Ya(2)=XYZ(2,NMF9(2,NFA))
  171. Xa(3)=XYZ(1,NMF9(3,NFA))
  172. Ya(3)=XYZ(2,NMF9(3,NFA))
  173.  
  174. a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2))
  175. AN(1,NUM(NMF9(1,NFA),K),K)=a
  176. AN(2,NUM(NMF9(1,NFA),K),K)=-1.D0
  177. ann(1)=(a*a+1.D0)**0.5D0
  178. AN(1,NUM(NMF9(1,NFA),K),K)=AN(1,NUM(NMF9(1,NFA),K),K)/ann(1)
  179. AN(2,NUM(NMF9(1,NFA),K),K)=AN(2,NUM(NMF9(1,NFA),K),K)/ann(1)
  180.  
  181. b=(Ya(1)-Ya(3))/(Xa(1)-Xa(3))
  182. AN(1,NUM(NMF9(2,NFA),K),K)=b
  183. AN(2,NUM(NMF9(2,NFA),K),K)=-1.D0
  184. ann(2)=(b*b+1.D0)**0.5D0
  185. AN(1,NUM(NMF9(2,NFA),K),K)=AN(1,NUM(NMF9(2,NFA),K),K)/ann(2)
  186. AN(2,NUM(NMF9(2,NFA),K),K)=AN(2,NUM(NMF9(2,NFA),K),K)/ann(2)
  187.  
  188. c=(Ya(3)-Ya(2))/(Xa(3)-Xa(2))
  189. AN(1,NUM(NMF9(3,NFA),K),K)=c
  190. AN(2,NUM(NMF9(3,NFA),K),K)=-1.D0
  191. ann(3)=(c*c+1.D0)**0.5D0
  192. AN(1,NUM(NMF9(3,NFA),K),K)=AN(1,NUM(NMF9(3,NFA),K),K)/ann(3)
  193. AN(2,NUM(NMF9(3,NFA),K),K)=AN(2,NUM(NMF9(3,NFA),K),K)/ann(3)
  194.  
  195. c write(6,*)'ann=',(ann(kk),kk=1,ns)
  196. DO 231 I=1,NS
  197. Xq=XYZ(1,N1F9(I,NFA))
  198. Yq=XYZ(2,N1F9(I,NFA))
  199. YP(N1F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I)
  200. Xq=XYZ(1,N2F9(I,NFA))
  201. Yq=XYZ(2,N2F9(I,NFA))
  202. YP(N2F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I)
  203. 231 CONTINUE
  204.  
  205. c write(6,*)'NMY=',(NMY(I),I=1,NY)
  206. c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
  207. c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
  208. c la normale n'est pas la meme pour tous les pts (SEG3)
  209. DO 232 I=1,NS
  210. AN(1,N1F9(I,NFA),K)=AN(1,NMS(I),K)
  211. AN(2,N1F9(I,NFA),K)=AN(2,NMS(I),K)
  212. AN(1,N2F9(I,NFA),K)=AN(1,NMS(I),K)
  213. AN(2,N2F9(I,NFA),K)=AN(2,NMS(I),K)
  214. 232 CONTINUE
  215. c._._ On Calcule la normale au SEG3 (NMS) et YP
  216. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  217.  
  218. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  219. c._._ On Calcule les normales aux CUB8,PRI6,TET4 et PYR5
  220. c._._ (La normale est la meme pour tous les pts)
  221. ELSEIF( ITYPEL.EQ.14
  222. & .OR.ITYPEL.EQ.16
  223. & .OR.ITYPEL.EQ.23
  224. & .OR.ITYPEL.EQ.25 )THEN
  225.  
  226. Xa(1)=XYZ(1,NMS(1))
  227. Ya(1)=XYZ(2,NMS(1))
  228. Za(1)=XYZ(3,NMS(1))
  229.  
  230. Xa(2)=XYZ(1,NMS(2))
  231. Ya(2)=XYZ(2,NMS(2))
  232. Za(2)=XYZ(3,NMS(2))
  233.  
  234. Xa(3)=XYZ(1,NMS(3))
  235. Ya(3)=XYZ(2,NMS(3))
  236. Za(3)=XYZ(3,NMS(3))
  237.  
  238. AN(1,NMS(1),K)=
  239. & (Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1))
  240. AN(2,NMS(1),K)=
  241. & (Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1))
  242. AN(3,NMS(1),K)=
  243. & (Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1))
  244. ann(1)=(AN(1,NMS(1),K)*AN(1,NMS(1),K)
  245. & +AN(2,NMS(1),K)*AN(2,NMS(1),K)
  246. & +AN(3,NMS(1),K)*AN(3,NMS(1),K))**0.5D0
  247. AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1)
  248. AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1)
  249. AN(3,NMS(1),K)=AN(3,NMS(1),K)/ann(1)
  250.  
  251. DO 431 I=2,NS
  252. AN(1,NMS(I),K)=AN(1,NMS(1),K)
  253. AN(2,NMS(I),K)=AN(2,NMS(1),K)
  254. AN(3,NMS(I),K)=AN(3,NMS(1),K)
  255. 431 CONTINUE
  256.  
  257. c write(6,*)'ann=',ann(1)
  258. DO 432 I=1,NY
  259. Xq=XYZ(1,NMY(I))
  260. Yq=XYZ(2,NMY(I))
  261. Zq=XYZ(3,NMY(I))
  262. YP(NMY(I),K)=ABS(
  263. & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K)
  264. & +(Zq-Za(1))*AN(3,NMS(1),K))
  265. c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)),
  266. c &' (Zq-Za(1))=',(Zq-Za(1))
  267. 432 CONTINUE
  268.  
  269. c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K)
  270.  
  271. c write(6,*)'NMY=',(NMY(I),I=1,NY)
  272. c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
  273. c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
  274. c la normale n'est pas la meme pour tous les pts (SEG3)
  275. DO 433 I=1,NY
  276. AN(1,NMY(I),K)=AN(1,NMS(1),K)
  277. AN(2,NMY(I),K)=AN(2,NMS(1),K)
  278. AN(3,NMY(I),K)=AN(3,NMS(1),K)
  279. 433 CONTINUE
  280. c._._ FIN Calcul des normales aux CUB8,PRI6,TET4 et PYR5
  281. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  282.  
  283. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  284. c._._ On Calcule les normales aux CU27,PR21,TE15 et PY19
  285. c._._ (La normale est la meme pour tous les pts)
  286. ELSEIF( ITYPEL.EQ.33
  287. & .OR.ITYPEL.EQ.34
  288. & .OR.ITYPEL.EQ.35
  289. & .OR.ITYPEL.EQ.36 )THEN
  290.  
  291. Xa(1)=XYZ(1,NMS(1))
  292. Ya(1)=XYZ(2,NMS(1))
  293. Za(1)=XYZ(3,NMS(1))
  294.  
  295. Xa(2)=XYZ(1,NMS(2))
  296. Ya(2)=XYZ(2,NMS(2))
  297. Za(2)=XYZ(3,NMS(2))
  298.  
  299. Xa(3)=XYZ(1,NMS(3))
  300. Ya(3)=XYZ(2,NMS(3))
  301. Za(3)=XYZ(3,NMS(3))
  302.  
  303. Xa(4)=XYZ(1,NMS(NS))
  304. Ya(4)=XYZ(2,NMS(NS))
  305. Za(4)=XYZ(3,NMS(NS))
  306.  
  307. ANx1=(Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1))
  308. ANy1=(Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1))
  309. ANz1=(Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1))
  310. ann(1)=(ANx1*ANx1+ANy1*ANy1+ANz1*ANz1)**0.5D0
  311.  
  312. ANx2=(Ya(2)-Ya(1))*(Za(4)-Za(1))-(Za(2)-Za(1))*(Ya(4)-Ya(1))
  313. ANy2=(Za(2)-Za(1))*(Xa(4)-Xa(1))-(Xa(2)-Xa(1))*(Za(4)-Za(1))
  314. ANz2=(Xa(2)-Xa(1))*(Za(4)-Za(1))-(Ya(2)-Ya(1))*(Xa(4)-Xa(1))
  315. ann(2)=(ANx2*ANx2+ANy2*ANy2+ANz2*ANz2)**0.5D0
  316.  
  317. IF(ann(1).GT.ann(2))THEN
  318. AN(1,NMS(1),K)=ANx1/ann(1)
  319. AN(2,NMS(1),K)=ANy1/ann(1)
  320. AN(3,NMS(1),K)=ANz1/ann(1)
  321. ELSE
  322. AN(1,NMS(1),K)=ANx2/ann(2)
  323. AN(2,NMS(1),K)=ANy2/ann(2)
  324. AN(3,NMS(1),K)=ANz2/ann(2)
  325. ann(1)=ann(2)
  326. ENDIF
  327.  
  328. DO 531 I=2,NS
  329. AN(1,NMS(I),K)=AN(1,NMS(1),K)
  330. AN(2,NMS(I),K)=AN(2,NMS(1),K)
  331. AN(3,NMS(I),K)=AN(3,NMS(1),K)
  332. 531 CONTINUE
  333.  
  334. c write(6,*)'ann=',ann(1)
  335. DO 532 I=1,NY
  336. Xq=XYZ(1,NMY(I))
  337. Yq=XYZ(2,NMY(I))
  338. Zq=XYZ(3,NMY(I))
  339. YP(NMY(I),K)=ABS(
  340. & (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K)
  341. & +(Zq-Za(1))*AN(3,NMS(1),K))
  342. c write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)),
  343. c &' (Zq-Za(1))=',(Zq-Za(1))
  344. 532 CONTINUE
  345.  
  346. c write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K)
  347.  
  348. c write(6,*)'NMS=',(NMS(I),I=1,NS)
  349. c write(6,*)'NMY=',(NMY(I),I=1,NY)
  350. c write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
  351. c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
  352. c la normale n'est pas la meme pour tous les pts (SEG3)
  353. DO 533 I=1,NY
  354. AN(1,NMY(I),K)=AN(1,NMS(1),K)
  355. AN(2,NMY(I),K)=AN(2,NMS(1),K)
  356. AN(3,NMY(I),K)=AN(3,NMS(1),K)
  357. 533 CONTINUE
  358. c._._ FIN Calcul des normales aux CU27,PR21,TE15 et PY19
  359. c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
  360. ELSE
  361.  
  362. write(6,*)'Type d elements non prevus dans FPU'
  363. write(6,*)'K=',K,' IDIM=',IDIM,' NS=',NS
  364.  
  365. RETURN
  366. ENDIF
  367.  
  368. c write(6,*)' FIN CALCUL NORMALES'
  369. c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY)
  370. DO 244 I=1,NP
  371. U0=0.D0
  372. U1=0.D0
  373. call initd(un,np,0.D0)
  374. DO 243 N=1,IDIM
  375. U0=U0+VS(I,N)*AN(N,I,K)
  376. U1=U1+(VS(I,N) - UN(I)*AN(N,I,K))**2.D0
  377. 243 CONTINUE
  378. UN(I)=U0
  379. UT(I)=U1**0.5D0 + 1.D-10
  380. 244 CONTINUE
  381. c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY)
  382.  
  383. c write(6,*)'AN=',AN(1,NMS(1),K),AN(2,NMS(1),K)
  384. c write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
  385. c write(6,*)'UN=',(UN(NMY(I)),I=1,NY)
  386. c write(6,*)'UT=',(UT(NMY(I)),I=1,NY)
  387. c write(6,*)'VS= '
  388. c write(6,1002)(VS(I,1),I=1,NP)
  389. c write(6,1002)(VS(I,2),I=1,NP)
  390.  
  391. cc CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
  392. cc & NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  393. c._._ On Calcule la matrice masse diagonale maillage appuye
  394. cc CALL INITD(DG,NPTA,0.D0)
  395. cc DO 221 I=1,NP
  396. cc AI=0.D0
  397. cc DO 220 LG=1,NPG
  398. cc PDR=PGSQ(LG)*DEUPI*RPG(LG)
  399. cc AI=AI+FN(I,LG)*PDR
  400. c220 CONTINUE
  401. cc JD=IPADA(NUM(I,K))
  402. cc DG(JD)=DG(JD)+AI
  403. c221 CONTINUE
  404. c._._ On Calcule la matrice masse diagonale maillage appuye
  405.  
  406. UETM=0.D0
  407. YPM =0.D0
  408. DO 225 II=1,NY
  409. I=NMY(II)
  410. IS=NMS(II)
  411. YP(I,K)=YP(I,K)+1.D-20
  412. c write(6,*)'YP ',YP(I,K)
  413.  
  414. JI=IPADT(NUM(I,K))
  415. c JS=IPADS(NUM(IS,K))
  416. JR=(1-IKR)*(JI-1)+1
  417. JM=(1-IKM)*(JI-1)+1
  418. c YPLUS=RO(JR)*YP(I,K)*UETK(I,K)/AMU(JM)
  419. YPLUS=RO(JR)*YP(I,K)*UT(I)/AMU(JM)
  420. YPLUS=YPLUS+1.D-10
  421.  
  422. c UET0=UETK(I,K)
  423. UET0=1.D-20
  424. c write(6,*)'M=0 YPLUS=',YPLUS,UT(I),AMU(JM),UET0
  425. w=0.5D0
  426.  
  427. DO 30 M=1,NBIT
  428. IF(YPLUS.LE.5.D0)THEN
  429. UPLUS=YPLUS
  430. UETN=UT(I)/UPLUS
  431. YPLUS=w*YPLUS + (1.D0-w)*RO(JR)*YP(I,K)*UETN/AMU(JM)
  432. ELSE
  433. c Reichardt
  434. UPLUS=1.D0/AKER*Log(1.D0+AKER*yplus)+
  435. & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0))
  436. UETN=UT(I)/UPLUS
  437. YPLUS=RO(JR)*YP(I,K)*UETN/AMU(JM)
  438. ENDIF
  439. resid=ABS(UETN-UET0)/UETN
  440. UET0=UETN
  441. IF(resid.LE.EPS)GO TO 31
  442. 30 CONTINUE
  443. write(6,*)'Operateur FPU, NON CONVERGENCE en 40 iterations: eps=',
  444. & eps,' residu ->',resid
  445. 31 CONTINUE
  446.  
  447. c write(6,*)' Convergence en ',m,' iterations'
  448. c UETK(I,K)=UETN
  449. c VET(JS)=UETN
  450. UETM = (UETM*FLOAT(II-1)+UETN)/FLOAT(II)
  451. YPM = (YPM *FLOAT(II-1)+YP(I,K))/FLOAT(II)
  452. AMUEF(I,K)=AMU(JM)
  453. IF(YPLUS.GT.5.D0)THEN
  454. AMUEF(I,K)=(RO(JR)*UETN*UETN*YP(I,K)/UT(I))+AMU(JM)
  455. ANUT(IPADT(NUM(I,K)))=AMUEF(I,K)
  456. ENDIF
  457.  
  458. 225 CONTINUE
  459. c write(6,*)' APRES 225 NY=',NY
  460.  
  461. DO 226 I=1,NS
  462. JI=IPADT(NUM(I,K))
  463. JM=(1-IKM)*(JI-1)+1
  464. AMUEF(I,K)=AMU(JM)
  465. ANUT(IPADT(NUM(NMS(I),K)))=AMUEF(I,K)
  466. UET(IPADS(NUM(NMS(I),K)))=UETM
  467. YPS(IPADS(NUM(NMS(I),K)))=YPM
  468. 226 CONTINUE
  469.  
  470. IF(NBINC.EQ.3)THEN
  471. DO 227 I=1,NS
  472. JI=IPADT(NUM(I,K))
  473. JM=(1-IKM)*(JI-1)+1
  474. TKE(IPADS(NUM(NMS(I),K)),1)=UETM*UETM/SQCNU
  475. TKE(IPADS(NUM(NMS(I),K)),2)=UETM*UETM*UETM/(AKER*YPM)
  476. 227 CONTINUE
  477. ENDIF
  478.  
  479. c write(6,*)'NS=',NS
  480. c write(6,1001)(LKKE(I),I=1,NS)
  481. c write(6,1001)(NUM(I,K),I=1,NP)
  482. c write(6,*)'UET K=',K,(UETK(I,K),I=1,NP)
  483. c write(6,*)'AMUEF K=',K,(AMUEF(I,K),I=1,NP)
  484.  
  485. c._._ On deduit aussi Tau aux sommets
  486. c._._ On deduit Nut aux points de Gauss
  487.  
  488. C.................................................................
  489.  
  490. C...... FIN Boucle sur les elements Indice K
  491. c write(6,*)'. FIN Boucle sur les elements Indice K ',K
  492. 2 CONTINUE
  493.  
  494. c write(6,*)'UET='
  495. c write(6,1001)(UET(i),i=1,npts)
  496.  
  497.  
  498. c write(6,*)' Nb d element traite KT=',KT
  499. c write(6,*)'RETOUR FPTAU'
  500. RETURN
  501. 1002 FORMAT(10(1X,1PE11.4))
  502. 1001 FORMAT(20(1X,I5))
  503. END
  504.  
  505.  
  506.  
  507.  

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