Télécharger ksprus.eso

Retour à la liste

Numérotation des lignes :

  1. C KSPRUS SOURCE PV 16/11/17 22:00:29 9180
  2. SUBROUTINE KSPRUS
  3. &(MELEME,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,
  4. & ANUKK,IDCEN,IKOMP,NPTU,IPADU,IZTU1,TETAN,NP,NEL,CMD)
  5.  
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. C*****************************************************************************
  9. C
  10. C Ce SP calcule les matrices elementaires de divergence alias C
  11. C
  12. C IKAS=1 KMCT calcul de Ct (Div U)
  13. C IKAS=2 KMAC calcul de C uniquement (Grad p)
  14. C IKAS=3 KCCT calcul de C assemblage pour C et Ct
  15. C
  16. C*****************************************************************************
  17. CHARACTER*8 NOM0
  18.  
  19. -INC CCOPTIO
  20. -INC CCGEOME
  21. -INC SMCOORD
  22. -INC SIZFFB
  23. POINTEUR IZF1.IZFFM
  24. -INC SMLENTI
  25. POINTEUR IPADU.MLENTI
  26. -INC SMELEME
  27. -INC SMCHPOI
  28. POINTEUR IZTGG1.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL
  29. -INC CCREEL
  30. PARAMETER (LRV=64)
  31. SEGMENT PETROV
  32. REAL*8 WT(LRV,NP,NPG,IDIM),WS(LRV,NP,NPG,IDIM),HK(LRV,IDIM,NP,NPG)
  33. REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG)
  34. REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV)
  35. REAL*8 AJK(LRV,IDIM,IDIM,NPG)
  36. ENDSEGMENT
  37. DIMENSION KIPM(3),XYZ1(24)
  38. C
  39. C OPERATEUR PRESSION
  40. C
  41. DEUPI=1.D0
  42. IF(IAXI.NE.0)DEUPI=2.D0*XPI
  43. IF(IKAS.EQ.3)DEUPI=-DEUPI
  44.  
  45. NBEL = NEL
  46. C write(6,*)' DEUPI=',deupi
  47. IF(IDIM.EQ.2)IPM3=IPM1
  48. KIPM(1)=IPM1
  49. KIPM(2)=IPM2
  50. KIPM(3)=IPM3
  51. SEGACT MELEME
  52. SEGACT IPM1*MOD,IPM2*MOD,IPM3*MOD
  53. SEGACT IPADU,IZTU1,TETAN,IZTGG1
  54.  
  55. IF(KPRE.NE.2)THEN
  56. C
  57. C ON TRAITE LE CAS DIV U -> T sommet
  58. C
  59. IF(INEFMD.EQ.3)THEN
  60. IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'PRP0'
  61. IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'PRP1'
  62. IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'PFP1'
  63. ELSEIF(INEFMD.EQ.2)THEN
  64. IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'MCP0'
  65. IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'MCP1'
  66. IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'MCF1'
  67. ELSEIF(INEFMD.EQ.4)THEN
  68. NOM0=NOMS(ITYPEL)//' '
  69. ENDIF
  70. c write(*,*)' Kpruss : NOM0=',nom0,' ikas=',ikas,' nel=',nel
  71. c write(*,*)' IK1',IK1,' COEF',IZTGG1.VPOCHA(1,1),'ANUK',ANUKK
  72. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  73.  
  74. SEGACT IZFFM*MOD
  75. IZHR=KZHR(1)
  76. SEGACT IZHR*MOD
  77. NES=GR(/1)
  78. NPG=GR(/3)
  79. IZF1=KTP(1)
  80. SEGACT IZF1*MOD
  81. * ON modifie pour avoir les bonnes fonctions de forme. GBM 18/10/99
  82. MP1=FN(/1)
  83. *
  84. * On fait des packets
  85. *
  86. NNN=MOD(NBEL,LRV)
  87. IF(NNN.EQ.0) THEN
  88. NPACK=NBEL/LRV
  89. ELSE
  90. NPACK=1+(NBEL-NNN)/LRV
  91. ENDIF
  92. KPACKD=1
  93. KPACKF=NPACK
  94. *
  95. * boucle sur les packets, on active la structure petrov
  96. *
  97. SEGINI PETROV
  98. DO 7001 KPACK=KPACKD,KPACKF
  99. C
  100. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  101. C
  102. C 1. Calcul des limites du paquet courant.
  103. KDEB=1+(KPACK-1)*LRV
  104. KFIN=MIN(NBEL,KDEB+LRV-1)
  105. C
  106. DO 7002 K=KDEB,KFIN
  107. KP=K-KDEB+1
  108. NK = K + K0
  109. NKR=(1-IK1)*(NK-1)+1
  110. COEFK(KP)=IZTGG1.VPOCHA(NKR,1) + 1.D-30
  111. * WRITE(*,*) IZTGG1.VPOCHA(NKR,1), COEFK(KP)
  112. ANUK(KP)=ANUKK/COEFK(KP)
  113. 7002 CONTINUE
  114. C
  115. C Une seule composante de projection car div se projete
  116. C sur un scalaire.
  117. C
  118. c write (*,*) 'IPADU', IPADU.LECT(/1)
  119. CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  120. & NES,NP,NPG,IAXI,XCOOR,
  121. & MELEME.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  122. & TETAN.VPOCHA,IPADU.LECT,IZTU1.VPOCHA,IPADU.LECT,NPTU,ANUK,
  123. & WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL,
  124. & DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  125.  
  126. DO 30 KE=KDEB,KFIN
  127. NK=KE-KDEB+1
  128.  
  129. * DO 12 I=1,NP
  130. * J=NUM(I,KE)
  131. * DO 12 N=1,IDIM
  132. * XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N)
  133. * 12 CONTINUE
  134.  
  135. * CALL CALJBR
  136. * &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  137. C write(6,*)' Retour caljbr ',mp1,np,npg
  138.  
  139.  
  140.  
  141. DO 324 K=1,IDIM
  142. IPM4=KIPM(K)
  143. DO 323 M=1,MP1
  144. DO 323 I=1,NP
  145. U=0.D0
  146. DO 333 L=1,NPG
  147. * ON modifie la formule GBM 18/10/99
  148. U=U+WT(NK,M,L,1)*HK(NK,K,I,L)*PGSK(NK,L)*DEUPI*RPGK(NK,L)
  149. * U=U+FN(M,L)*HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L)
  150. 333 CONTINUE
  151.  
  152. IF(IAXI.NE.0.AND.K.EQ.1)THEN
  153. DO 334 L=1,NPG
  154. * modifié par GBM 18/10/99
  155. U=U+WT(NK,M,L,1)*FN(I,L)*PGSK(NK,L)*DEUPI
  156. * U=U+FN(M,L)*FN(I,L)*PGSQ(L)*DEUPI
  157. 334 CONTINUE
  158. ENDIF
  159. C
  160. C K1=1+(1-IK1)*(NK-1)
  161. C attention valable que pour coef au centre
  162. U=U*COEFK(NK)
  163. * WRITE(*,*) NK,M,I,COEFK(NK)
  164. C
  165. if(ikas.ne.2)then
  166. IPM4.AM(KE,I,M)=IPM4.AM(KE,I,M)+U
  167. else
  168. IPM4.AM(KE,M,I)=IPM4.AM(KE,M,I)+U
  169. endif
  170.  
  171. * write(*,*) KE,PGSQ
  172. 323 CONTINUE
  173. 324 CONTINUE
  174. 30 CONTINUE
  175. 7001 CONTINUE
  176.  
  177. SEGSUP IZHR,IZFFM,PETROV
  178.  
  179.  
  180. C CAS MACRO CENTRE
  181.  
  182. ELSEIF(KPRE.EQ.2)THEN
  183. * initialisation de npg au pif par TC
  184. NPG=1
  185. SEGINI PETROV
  186. NOM0=NOMS(ITYPEL)//' '
  187. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  188.  
  189. SEGACT IZFFM*MOD
  190. IZHR=KZHR(1)
  191. SEGACT IZHR*MOD
  192. NES=GR(/1)
  193. NPG=GR(/3)
  194. IZF1=KTP(1)
  195. SEGACT IZF1*MOD
  196. MPG=IZF1.FN(/2)
  197. NP=GR(/2)
  198.  
  199. DO 40 KE=1,NEL
  200.  
  201. NK=K0+KE
  202.  
  203. IX=0
  204. DO 42 I=1,NP
  205. J=MELEME.NUM(I,KE)
  206. DO 42 N=1,IDIM
  207. IX=IX+1
  208. XYZ1(IX)=XCOOR((J-1)*(IDIM+1) +N)
  209. 42 CONTINUE
  210.  
  211. CALL CALJBR(FN,GR,PG,XYZ1,HR,PGSQ,RPG,NES,
  212. & IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  213.  
  214. DO 424 K=1,IDIM
  215. IPM4=KIPM(K)
  216.  
  217. DO 423 I=1,NP
  218. U=0.D0
  219. DO 433 L=1,NPG
  220. U=U+HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L)
  221. 433 CONTINUE
  222.  
  223. IF(IAXI.NE.0.AND.K.EQ.1)THEN
  224. DO 434 L=1,NPG
  225. U=U+FN(I,L)*PGSQ(L)*DEUPI
  226. 434 CONTINUE
  227. ENDIF
  228.  
  229. K1=1+(1-IK1)*(NK-1)
  230. U=U*IZTGG1.VPOCHA(K1,1)
  231.  
  232. if(ikas.ne.2)then
  233. IPM4.AM(KE,I,1)=IPM4.AM(KE,I,1)+U
  234. else
  235. IPM4.AM(KE,1,I)=IPM4.AM(KE,1,I)+U
  236. endif
  237.  
  238. 423 CONTINUE
  239. 424 CONTINUE
  240.  
  241. 40 CONTINUE
  242.  
  243. SEGSUP IZHR,IZFFM,PETROV
  244.  
  245. ENDIF
  246.  
  247. RETURN
  248. 1002 FORMAT(10(1X,1PE11.4))
  249. 1040 FORMAT(1X,'CALCUL MATRICE AM ',I4/10(1X,1PE11.4))
  250. END
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  

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