Télécharger chaga1.eso

Retour à la liste

Numérotation des lignes :

chaga1
  1. C CHAGA1 SOURCE CB215821 20/11/04 21:15:28 10766
  2.  
  3. C=======================================================================
  4. C= C H A G A 1 =
  5. C= ----------- =
  6. C= =
  7. C= Fonction : =
  8. C= ---------- =
  9. C= =
  10. C= Creation du MCHAML de SOURCE de chaleur GAUSSIENNE decrit par les =
  11. C= caracteristiques du modele passees en entree. =
  12. C= =
  13. C= =
  14. C= Entrees : =
  15. C= --------- =
  16. C= =
  17. C= IPMCHA : pointeur sur segment MCHAML de la sous-zone traitee =
  18. C= =
  19. C= IPMAIL : pointeur sur segment MELEME du maillage ou creer la =
  20. C= source de chaleur =
  21. C= =
  22. C= IPINTE : pointeur sur segment MINTE de l'element finis =
  23. C= =
  24. C= IK : indicateur du "type" de distribution gaussienne : =
  25. C= IK = 1 : source gaussienne ISOTROPE =
  26. C= IK = 2 : source gaussienne ISOTROPE-TRANSVERSE =
  27. C= =
  28. C= =
  29. C= Sortie : =
  30. C= -------- =
  31. C= =
  32. C= IPQGAU : pointeur sur segment MELVAL des valeurs de la source de =
  33. C= chaleur GAUSSIENNE =
  34. C= =
  35. C= XQ0 : valeur de la quantite de chaleur imposee, utilisee en =
  36. C= sortie dans chamas pour "renormaliser" le champ a QTOT =
  37. C= =
  38. C=======================================================================
  39.  
  40. SUBROUTINE CHAGA1(IPMCHA,IPMAIL,IPINTE,IK,IPQGAU,XQ0)
  41.  
  42. IMPLICIT INTEGER(I-N)
  43. IMPLICIT REAL*8(A-H,O-Z)
  44.  
  45. -INC PPARAM
  46. -INC CCOPTIO
  47. -INC CCREEL
  48. -INC SMCOORD
  49. -INC SMELEME
  50. -INC SMCHAML
  51. -INC SMINTE
  52.  
  53. REAL*8 XCO1(3),XCD1(3),XCX1(3)
  54. CHARACTER*8 TYPC1
  55. CHARACTER*4 MOCAR(5)
  56. CHARACTER*(LOCOMP) NOMC1
  57.  
  58. DATA MOCAR /'QTOT','ORIG','RGAU','DIRE','ZGAU'/
  59.  
  60.  
  61. IPQGAU = 0
  62.  
  63. C==== Recuperation des caracteristiques de la source de chaleur
  64. C
  65. C Initialisation des valeurs des parametres de la source
  66. C XQ0 : valeur de QTOT
  67. C NPO1, NPD1 : numeros points origine & direction
  68. C RX0, ZX0 : valeurs de RGAU et ZGAU
  69. XQ0 = 0.D0
  70. NPO1 = 0
  71. NPD1 = 0
  72. RX0 = 0.D0
  73. ZX0 = 0.D0
  74. C
  75. C NCAR1 : nombre de caracteristiques a lire
  76. NCAR1 = 3
  77. IF (IK.EQ.2) NCAR1 = 5
  78. C
  79. MCHAM1=IPMCHA
  80. C SEGACT,MCHAM1
  81. NCO1 = MCHAM1.IELVAL(/1)
  82. DO 100 I=1,NCAR1
  83. NOMC1 = MOCAR(I)
  84. CALL PLACE(MCHAM1.NOMCHE,NCO1,IPLAC,NOMC1)
  85. c write(6,*) 'NONC1 =',NOMC1
  86. c write(6,*) 'MCHAM1.NOMCHE =',(MCHAM1.NOMCHE(ii),ii=1,NCO1)
  87. IF (IPLAC.EQ.0) THEN
  88. GOTO 997
  89. ELSE
  90. TYPC1 = MCHAM1.TYPCHE(IPLAC)(1:8)
  91. MELVA1 = MCHAM1.IELVAL(IPLAC)
  92. IF (MELVA1.VELCHE(/1).GT.1
  93. & .OR.MELVA1.VELCHE(/2).GT.1) GOTO 999
  94. C
  95. C Lecture QTOT
  96. IF (I.EQ.1) THEN
  97. IF (TYPC1.NE.'REAL*8 ') GOTO 998
  98. SEGACT,MELVA1
  99. XQ0 = MELVA1.VELCHE(1,1)
  100. C SEGDES,MELVA1
  101. C
  102. C Lecture ORIG
  103. ELSEIF (I.EQ.2) THEN
  104. IF (TYPC1.NE.'POINTEUR') GOTO 998
  105. MELVA1 = MCHAM1.IELVAL(IPLAC)
  106. SEGACT,MELVA1
  107. NPO1 = MELVA1.IELCHE(1,1)
  108. C SEGDES,MELVA1
  109. C
  110. C Lecture RGAU
  111. ELSEIF (I.EQ.3) THEN
  112. IF (TYPC1.NE.'REAL*8 ') GOTO 998
  113. MELVA1 = MCHAM1.IELVAL(IPLAC)
  114. SEGACT,MELVA1
  115. RX0 = MELVA1.VELCHE(1,1)
  116. C SEGDES,MELVA1
  117. C
  118. C Lecture DIRE
  119. ELSEIF (I.EQ.4) THEN
  120. IF (TYPC1.NE.'POINTEUR') GOTO 998
  121. MELVA1 = MCHAM1.IELVAL(IPLAC)
  122. SEGACT,MELVA1
  123. NPD1 = MELVA1.IELCHE(1,1)
  124. C SEGDES,MELVA1
  125. C
  126. C Lecture ZGAU
  127. ELSEIF (I.EQ.5) THEN
  128. IF (TYPC1.NE.'REAL*8 ') GOTO 998
  129. MELVA1 = MCHAM1.IELVAL(IPLAC)
  130. SEGACT,MELVA1
  131. ZX0 = MELVA1.VELCHE(1,1)
  132. C SEGDES,MELVA1
  133. ENDIF
  134. ENDIF
  135. 100 CONTINUE
  136. C SEGDES,MCHAM1
  137. C
  138. C Remplissage du vecteur de coordonnes de ORIG
  139. IDIMP1 = IDIM+1
  140. IF (NPO1.EQ.0.OR.NPO1.GT.nbpts) THEN
  141. NOMC1='ORIG'
  142. GOTO 998
  143. ELSE
  144. DO 120 ID1=1,IDIM
  145. XCO1(ID1)=XCOOR((NPO1-1)*IDIMP1+ID1)
  146. 120 CONTINUE
  147. ENDIF
  148. C
  149. C Remplissage du vecteur de coordonnes de DIRE
  150. C XN1 : norme de DIRE au carre
  151. XN1 = 1.D0
  152. IF (IK.EQ.1) THEN
  153. ZX0 = RX0
  154. XCD1(1)=1.D0
  155. IF (IDIM.GT.1) THEN
  156. DO 130 ID1=2,IDIM
  157. XCD1(ID1)=0.D0
  158. 130 CONTINUE
  159. ENDIF
  160. ELSEIF (IK.EQ.2) THEN
  161. XN1 = 0.D0
  162. DO 140 ID1=1,IDIM
  163. XCD1(ID1)=XCOOR((NPD1-1)*IDIMP1+ID1)
  164. XN1 = XCD1(ID1) ** 2 + XN1
  165. 140 CONTINUE
  166. ELSE
  167. WRITE(IOIMP,*) ' *** Dans CHAGA1, cas (IK) non prevu'
  168. CALL ERREUR(5)
  169. RETURN
  170. ENDIF
  171. C write(6,*) 'XQ0 =',XQ0
  172. C write(6,*) 'RX0 =',RX0
  173. C write(6,*) 'ZX0 =',ZX0
  174. C write(6,*) 'NPO1 =',NPO1
  175. C write(6,*) 'NPD1 =',NPD1
  176. C write(6,*) 'XCO1 =',(XCO1(i),i=1,IDIM)
  177. C write(6,*) 'XCD1 =',(XCD1(i),i=1,IDIM)
  178.  
  179.  
  180. C==== Construction du MELVAL de la distribution Gaussienne de chaleur
  181.  
  182. C Activation du maillage
  183. IPT1 = IPMAIL
  184. SEGACT,IPT1
  185. NBPT1 = IPT1.NUM(/1)
  186. NBEL1 = IPT1.NUM(/2)
  187.  
  188. C Petite verif. sous-zone elementaire
  189. NBS1 = IPT1.LISOUS(/1)
  190. IF (NBS1.NE.0) THEN
  191. WRITE(6,*) 'Dans CHAGA1 : le maillage a plusieurs sous-zones ?'
  192. CALL ERREUR(21)
  193. RETURN
  194. ENDIF
  195. C
  196. C Activation du segment MINTE
  197. C write(6,*) 'IPINTE=',IPINTE
  198. MINTE1 = IPINTE
  199. SEGACT,MINTE1
  200.  
  201. C Creation du MELVAL
  202. N1PTEL = MINTE1.POIGAU(/1)
  203. N1EL = NBEL1
  204. N2PTEL = 0
  205. N2EL = 0
  206. SEGINI,MELVA2
  207.  
  208. C Facteur de normalisation de la puissance :
  209. IF (IDIM.EQ.2.AND.IFOUR.NE.0) THEN
  210. XQ1 = 4.D0 * XQ0 / (RX0 * ZX0 * XPI)
  211. ELSEIF (IDIM.EQ.2.AND.IFOUR.EQ.0) THEN
  212. C write (6,*) ' *** Cas axis'
  213. RC2 = SQRT (2.D0)
  214. RCPI = SQRT (XPI)
  215. XIK1 = 0.25D0 * RX0 * RX0
  216. XIK1 = XIK1 * EXP(-2.D0 * XCO1(1) * XCO1(1) / RX0 / RX0)
  217. XIK2 = 0.5D0 * RX0 * XCO1(1) * RCPI / RC2
  218. XIK3 = RX0 * XCO1(1) * RCPI / (RC2 * RC2 * RC2)
  219. XIK3 = XIK3 * ERF(XCO1(1) * RC2 / RX0)
  220. XQ1 = XPI * ZX0 * RCPI / RC2
  221. XQ1 = XQ1 * (XIK1 + XIK2 + XIK3)
  222. XQ1 = XQ0 / XQ1
  223. ELSEIF (IDIM.EQ.3) THEN
  224. RX1 = RX0 * RX0 * ZX0
  225. XQ1 = SQRT ((2.D0 ** 5) / (XPI ** 3))
  226. XQ1 = XQ1 * XQ0 / RX1
  227. ELSE
  228. GOTO 996
  229. ENDIF
  230.  
  231. C Boucle sur les elements du maillage
  232. DO 200 IEL1=1,NBEL1
  233. C
  234. C Boucle sur les pts de Gauss
  235. C et calcul de la fonction Gaussienne
  236. C RR1 : distance a l'origine ORIG au carre
  237. C ZX1 : distance a la direction DIRE au carre
  238. C RX1 : rayon dans le plan orthog. a DIRE au carre
  239. DO 201 IG1=1,N1PTEL
  240. RR1 = 0.D0
  241. RX1 = 0.D0
  242. ZX1 = 0.D0
  243. DO 203 ID1=1,IDIM
  244. C XPG1 : coordonnees ID1 du Pt de Gauss
  245. XPG1 = 0.D0
  246. DO 202 JPT1=1,NBPT1
  247. NPTJM1 = IPT1.NUM(JPT1,IEL1)-1
  248. XPTJ1 = XCOOR(NPTJM1*IDIMP1+ID1)
  249. VNJ1 = MINTE1.SHPTOT(1,JPT1,IG1)
  250. XPG1 = XPG1 + XPTJ1*VNJ1
  251. 202 CONTINUE
  252. XCX1(ID1) = XPG1 - XCO1(ID1)
  253. RR1 = XCX1(ID1) ** 2 + RR1
  254. ZX1 = XCX1(ID1) * XCD1(ID1) + ZX1
  255. 203 CONTINUE
  256. ZX1 = ZX1 * ZX1 / XN1
  257. RX1 = RR1 - ZX1
  258. c write(6,*) 'RR1,ZX1,RX1=',RR1,ZX1,RX1
  259. c write(6,*) 'XCX1(i)=',(XCX1(i),i=1,idim)
  260. XS1 = RX1 / RX0 / RX0
  261. XS1 = ZX1 / ZX0 / ZX0 + XS1
  262. XS1 = EXP(-2.D0 * XS1)
  263. MELVA2.VELCHE(IG1,IEL1) = XS1 * XQ1
  264. 201 CONTINUE
  265. 200 CONTINUE
  266. C SEGDES,IPT1,MINTE1,MELVA2
  267. C
  268. C==== Affectation du resultat et sortie
  269. IPQGAU = MELVA2
  270. C
  271. RETURN
  272.  
  273. C==== Gestion erreurs et fin
  274.  
  275. C Option indisponible
  276. 996 CONTINUE
  277. INTERR(1) = IDIM
  278. CALL ERREUR(709)
  279. RETURN
  280.  
  281. C Il manque la donnee d'une composante
  282. 997 CONTINUE
  283. CALL ERREUR(80)
  284. RETURN
  285.  
  286. C Le type de la composante NOMC1 n'est pas celui attendu
  287. 998 CONTINUE
  288. MOTERR(1:8)=NOMC1
  289. CALL ERREUR(679)
  290. RETURN
  291.  
  292. C La composante est attendue constante par sous-zones
  293. 999 CONTINUE
  294. MOTERR(1:4)=NOMC1
  295. MOTERR(5:20)='CARACTERISTIQUES'
  296. CALL ERREUR(1065)
  297. RETURN
  298.  
  299. END
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  

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