Télécharger shapx.eso

Retour à la liste

Numérotation des lignes :

  1. C SHAPX SOURCE BP208322 18/04/12 21:15:12 9802
  2. c
  3. SUBROUTINE SHAPX(MELE,LV7,NBRMAX,NBENRJ,TLREEL,SHP,IRET)
  4. C
  5. C=======================================================================
  6. C
  7. C Entree: SHP = FONCTIONS DE FORME STANDARD ET LEUR DERIVEES
  8. C LV7 = FONCTIONS DE NIVEAU ET LEUR DERIVEES
  9. C LV7(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  10. C Sortie: SHP = TOUTES LES FONCTIONS DE FORME DE L'EF (std + enrichissement)
  11. c ET LEUR DERIVEES DANS LE REPERE GLOBAL
  12. C
  13. C
  14. C=====PARTIE DECLARATIVE================================================
  15. C
  16. implicit real*8(a-h,o-z)
  17. implicit integer(i-n)
  18. -INC CCOPTIO
  19. -INC CCREEL
  20. -INC SMLREEL
  21. C
  22. DIMENSION SHP(6,*)
  23. c PARAMETER (NBRMAX=5)
  24. REAL*8 LV7(NBRMAX,2,6,*)
  25. INTEGER TLREEL(NBRMAX,*)
  26. c DIMENSION LV7(NBRMAX,2,6,*)
  27. c DIMENSION TLREEL(NBRMAX,*)
  28. REAL*8 PSI,PHI,H,R,THETA
  29. C
  30. C
  31. c WRITE(*,*) '#### ENTREE DANS SHAPX ####'
  32. C
  33. c
  34. C=====AIGUILLAGE selon element==========================================
  35. c
  36. if(MELE.eq.263) GOTO 263
  37. if(MELE.eq.264) GOTO 264
  38. c Erreur pour les elements non programmés
  39. call erreur (21)
  40. RETURN
  41. C
  42. CTY if(MELE.eq.263) GOTO 263
  43. CTY Erreur pour les elements non programmés
  44. CTY call erreur (21)
  45. CTY RETURN
  46. c
  47. c
  48. C=====XQ4R==============================================================
  49. c
  50. 263 continue
  51. NBNN = 4
  52. c KMAX =LV7(/1)
  53. c
  54. c
  55. C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++
  56. KENR = 1
  57.  
  58. C>>>> boucle sur les noeuds
  59. DO 200 I=1,NBNN
  60.  
  61. MLREEL=TLREEL(KENR,I)
  62. if(MLREEL.eq.0) goto 200
  63.  
  64. C on calcule les fonctions utiles
  65. PHI = LV7(KENR,2,1,I)
  66. H = SIGN(1.D0,PHI)
  67. C
  68. c on modifie les fonctions de forme et leurs derivees
  69. II = (KENR*NBNN) + I
  70. DO IND=1,3
  71. SHP(IND,II) = H*SHP(IND,I)
  72. ENDDO
  73. C
  74. 200 CONTINUE
  75. C<<<<<<<<<<<<<
  76. C
  77. C
  78. C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++
  79. C
  80. c DO 300 KENR=2,NBRMAX
  81. DO 300 KENR=2,NBENRJ
  82. IENR = 4*(KENR-2) + 2
  83.  
  84. C>>>>>>> boucle sur les noeuds
  85. DO 305 I=1,NBNN
  86. C
  87. MLREEL=TLREEL(KENR,I)
  88. c write(*,*) 'KENR,I=',KENR,I,' MLREEL=',MLREEL
  89. if(MLREEL.eq.0) goto 305
  90.  
  91. C on recupere les level set et leur derivees
  92. PSI = LV7(KENR,1,1,I)
  93. PSIX = LV7(KENR,1,2,I)
  94. PSIY = LV7(KENR,1,3,I)
  95. PHI = LV7(KENR,2,1,I)
  96. PHIX = LV7(KENR,2,2,I)
  97. PHIY = LV7(KENR,2,3,I)
  98.  
  99. C on calcule les fonctions utiles
  100. H = SIGN(1.D0,PHI)
  101. R = ( (PSI**2) + (PHI**2) )**0.5
  102. THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI))))
  103. SIN1T = SIN(THETA)
  104. COS1T = COS(THETA)
  105. SIN05T = SIN(THETA/2.)
  106. COS05T = COS(THETA/2.)
  107. C
  108. C et leurs derivées
  109. C RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5)
  110. C RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5)
  111. C THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2)
  112. C THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2)
  113.  
  114. c on calcule les fonctions de forme d ENRICHISSEMENT + derivees
  115.  
  116. C+++ Fb = (r**0.5)*sin(theta/2)
  117. II = (IENR*NBNN) + I
  118.  
  119. C Fonction de forme Ni*Fb
  120. FB = (R**0.5) * SIN05T
  121. SHP(1,II) = SHP(1,I) * FB
  122. C Derivees (Ni*Fb),X = (Ni,X*Fb) + (Ni*Fb,X)
  123. SHP(2,II) = ( SHP(2,I) * FB ) +
  124. $ ( SHP(1,I)*0.5/(R**0.5) *
  125. $ ( ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  126. $ + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  127. C Derivees (Ni*Fb),y = (Ni,y*Fb) + (Ni*Fb,y)
  128. SHP(3,II) = ( SHP(3,I) * FB ) +
  129. $ ( SHP(1,I)*0.5/(R**0.5) *
  130. $ ( ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  131. $ + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  132.  
  133.  
  134. C+++ Fc = (r**0.5)*cos(theta/2)
  135. II = ((IENR+1)*NBNN) + I
  136.  
  137. C Fonction de forme Ni*Fc
  138. FC = (R**0.5) * COS05T
  139. SHP(1,II) = SHP(1,I) * FC
  140. C Derivees (Ni*Fc),x = (Ni,x*Fc) + (Ni*Fc,x)
  141. SHP(2,II) = ( SHP(2,I) * FC ) +
  142. $ ( SHP(1,I)*0.5/(R**0.5) *
  143. $ ( ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  144. $ - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  145. C Derivees (Ni*Fc),y = (Ni,y*Fc) + (Ni*Fc,y)
  146. SHP(3,II) = ( SHP(3,I) * FC) +
  147. $ ( SHP(1,I)*0.5/(R**0.5) *
  148. $ ( ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  149. $ - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  150.  
  151.  
  152. C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta)
  153. II = ((IENR+2)*NBNN) + I
  154.  
  155. C Fonction de forme Ni*Fd
  156. FD = (R**0.5) * SIN05T * SIN1T
  157. SHP(1,II)= SHP(1,I) * FD
  158. C Derivees (Ni*Fc),X = (Ni,X*Fc) + (Ni*Fc,X)
  159. SHP(2,II) = ( SHP(2,I) * FD ) +
  160. $ ( SHP(1,I)*0.5/(R**0.5) *
  161. $ ( ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  162. $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
  163. $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  164. C Derivees (Ni*Fc),Y = (Ni,Y*Fc) + (Ni*Fc,Y)
  165. SHP(3,II) = ( SHP(3,I) * FD ) +
  166. $ ( SHP(1,I)*0.5/(R**0.5) *
  167. $ ( ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  168. $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
  169. $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  170.  
  171.  
  172. C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta)
  173. II = ((IENR+3)*NBNN) + I
  174.  
  175. C Fonction de forme Ni*Fe
  176. FE = (R**0.5) * COS05T * SIN1T
  177. SHP(1,II)= SHP(1,I) * FE
  178. C Derivees (Ni*Fc),x = (Ni,x*Fe) + (Ni*Fe,x)
  179. SHP(2,II) = ( SHP(2,I) * FE ) +
  180. $ ( SHP(1,I)*0.5/(R**0.5) *
  181. $ ( ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  182. $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
  183. $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  184. C Derivees (Ni*Fe),y = (Ni,y*Fe) + (Ni*Fe,y)
  185. SHP(3,II) = ( SHP(3,I) * FE ) +
  186. $ ( SHP(1,I)*0.5/(R**0.5) *
  187. $ ( ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  188. $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
  189. $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  190.  
  191. C
  192. C
  193. 305 CONTINUE
  194. C<<<<<<<<<<<<<
  195.  
  196. 300 CONTINUE
  197.  
  198. GOTO 999
  199.  
  200. C_____________________________________________
  201. C FONCTIONS DE FORME TRIDIMENSIONNELLES
  202. C
  203. C=====XC8R==============================================================
  204. c
  205. 264 continue
  206. C
  207. NBNN = 8
  208. c KMAX =LV7(/1)
  209. c
  210. c
  211. C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++
  212. KENR = 1
  213.  
  214. C>>>> boucle sur les noeuds
  215. DO 400 I=1,NBNN
  216.  
  217. MLREEL=TLREEL(KENR,I)
  218. if(MLREEL.eq.0) goto 400
  219.  
  220. C on calcule les fonctions utiles
  221. PHI = LV7(KENR,2,1,I)
  222. H = SIGN(1.D0,PHI)
  223. C
  224. c on modifie les fonctions de forme et leurs derivees
  225. II = (KENR*NBNN) + I
  226. DO IND=1,4
  227. SHP(IND,II) = H*SHP(IND,I)
  228. ENDDO
  229. C
  230. 400 CONTINUE
  231. C<<<<<<<<<<<<<
  232. C
  233. C
  234. C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++
  235. C
  236. c DO 500 KENR=2,NBRMAX
  237. DO 500 KENR=2,NBENRJ
  238. IENR = 4*(KENR-2) + 2
  239.  
  240. C>>>>>>> boucle sur les noeuds
  241. DO 505 I=1,NBNN
  242. C
  243. MLREEL=TLREEL(KENR,I)
  244. c write(*,*) 'KENR,I=',KENR,I,' MLREEL=',MLREEL
  245. if(MLREEL.eq.0) goto 505
  246.  
  247. C on recupere les level set et leur derivees
  248. PSI = LV7(KENR,1,1,I)
  249. PSIX = LV7(KENR,1,2,I)
  250. PSIY = LV7(KENR,1,3,I)
  251. PSIZ = LV7(KENR,1,4,I)
  252. PHI = LV7(KENR,2,1,I)
  253. PHIX = LV7(KENR,2,2,I)
  254. PHIY = LV7(KENR,2,3,I)
  255. PHIZ = LV7(KENR,2,4,I)
  256.  
  257. C on calcule les fonctions utiles
  258. H = SIGN(1.D0,PHI)
  259. R = ( (PSI**2) + (PHI**2) )**0.5
  260. THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI))))
  261. SIN1T = SIN(THETA)
  262. COS1T = COS(THETA)
  263. SIN05T = SIN(THETA/2.)
  264. COS05T = COS(THETA/2.)
  265. C
  266. C et leurs derivées
  267. C RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5)
  268. C RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5)
  269. C THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2)
  270. C THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2)
  271.  
  272. c on calcule les fonctions de forme d ENRICHISSEMENT + derivees
  273.  
  274. C+++ Fb = (r**0.5)*sin(theta/2)
  275. II = (IENR*NBNN) + I
  276.  
  277. C Fonction de forme Ni*Fb
  278. FB = (R**0.5) * SIN05T
  279. SHP(1,II) = SHP(1,I) * FB
  280. C Derivees (Ni*Fb),X = (Ni,X*Fb) + (Ni*Fb,X)
  281. SHP(2,II) = ( SHP(2,I) * FB ) +
  282. $ ( SHP(1,I)*0.5/(R**0.5) *
  283. $ ( ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  284. $ + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  285. C Derivees (Ni*Fb),y = (Ni,y*Fb) + (Ni*Fb,y)
  286. SHP(3,II) = ( SHP(3,I) * FB ) +
  287. $ ( SHP(1,I)*0.5/(R**0.5) *
  288. $ ( ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  289. $ + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  290. C Derivees (Ni*Fb),z = (Ni,z*Fb) + (Ni*Fb,z)
  291. SHP(4,II) = ( SHP(4,I) * FB ) +
  292. $ ( SHP(1,I)*0.5/(R**0.5) *
  293. $ ( ( SIN05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
  294. $ + ( COS05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
  295. C
  296. C+++ Fc = (r**0.5)*cos(theta/2)
  297. II = ((IENR+1)*NBNN) + I
  298.  
  299. C Fonction de forme Ni*Fc
  300. FC = (R**0.5) * COS05T
  301. SHP(1,II) = SHP(1,I) * FC
  302. C Derivees (Ni*Fc),x = (Ni,x*Fc) + (Ni*Fc,x)
  303. SHP(2,II) = ( SHP(2,I) * FC ) +
  304. $ ( SHP(1,I)*0.5/(R**0.5) *
  305. $ ( ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  306. $ - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  307. C Derivees (Ni*Fc),y = (Ni,y*Fc) + (Ni*Fc,y)
  308. SHP(3,II) = ( SHP(3,I) * FC) +
  309. $ ( SHP(1,I)*0.5/(R**0.5) *
  310. $ ( ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  311. $ - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  312. C Derivees (Ni*Fc),z = (Ni,z*Fc) + (Ni*Fc,z)
  313. SHP(4,II) = ( SHP(4,I) * FC) +
  314. $ ( SHP(1,I)*0.5/(R**0.5) *
  315. $ ( ( COS05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
  316. $ - ( SIN05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
  317.  
  318.  
  319. C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta)
  320. II = ((IENR+2)*NBNN) + I
  321.  
  322. C Fonction de forme Ni*Fd
  323. FD = (R**0.5) * SIN05T * SIN1T
  324. SHP(1,II)= SHP(1,I) * FD
  325. C Derivees (Ni*Fd),X = (Ni,X*Fd) + (Ni*Fd,X)
  326. SHP(2,II) = ( SHP(2,I) * FD ) +
  327. $ ( SHP(1,I)*0.5/(R**0.5) *
  328. $ ( ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  329. $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
  330. $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  331. C Derivees (Ni*Fd),Y = (Ni,Y*Fd) + (Ni*Fd,Y)
  332. SHP(3,II) = ( SHP(3,I) * FD ) +
  333. $ ( SHP(1,I)*0.5/(R**0.5) *
  334. $ ( ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  335. $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
  336. $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  337. C Derivees (Ni*Fd),z = (Ni,z*Fd) + (Ni*Fd,z)
  338. SHP(4,II) = ( SHP(4,I) * FD ) +
  339. $ ( SHP(1,I)*0.5/(R**0.5) *
  340. $ ( ( SIN05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
  341. $ + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
  342. $ ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
  343.  
  344.  
  345. C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta)
  346. II = ((IENR+3)*NBNN) + I
  347.  
  348. C Fonction de forme Ni*Fe
  349. FE = (R**0.5) * COS05T * SIN1T
  350. SHP(1,II)= SHP(1,I) * FE
  351. C Derivees (Ni*Fe),x = (Ni,x*Fe) + (Ni*Fe,x)
  352. SHP(2,II) = ( SHP(2,I) * FE ) +
  353. $ ( SHP(1,I)*0.5/(R**0.5) *
  354. $ ( ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
  355. $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
  356. $ ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
  357. C Derivees (Ni*Fe),y = (Ni,y*Fe) + (Ni*Fe,y)
  358. SHP(3,II) = ( SHP(3,I) * FE ) +
  359. $ ( SHP(1,I)*0.5/(R**0.5) *
  360. $ ( ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
  361. $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
  362. $ ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
  363. C Derivees (Ni*Fe),z = (Ni,z*Fe) + (Ni*Fe,z)
  364. SHP(4,II) = ( SHP(4,I) * FE ) +
  365. $ ( SHP(1,I)*0.5/(R**0.5) *
  366. $ ( ( COS05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
  367. $ - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
  368. $ ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
  369.  
  370. C
  371. C
  372. 505 CONTINUE
  373. C<<<<<<<<<<<<<
  374.  
  375. 500 CONTINUE
  376.  
  377.  
  378. C_____________________________________________
  379. C FIN DU PROGRAMME
  380. 999 CONTINUE
  381. C
  382.  
  383. IRET = 1
  384. c WRITE(*,*) 'C EST FINI POUR SHAPX'
  385.  
  386. RETURN
  387. END
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  

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