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

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