Télécharger chaga1.eso

Retour à la liste

Numérotation des lignes :

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

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