Télécharger smtp.eso

Retour à la liste

Numérotation des lignes :

smtp
  1. C SMTP SOURCE CB215821 24/04/12 21:17:15 11897
  2. SUBROUTINE SMTP
  3. C-----------------------------------------------------------------------
  4. C 1) Calcul du second membre du systeme en trace de charge dans le cas
  5. C de la résolution des équations de Darcy par EFMH avec le modèle DARCY.
  6. C 2) Calcul du second membre du systeme en trace de charge et/ou d'une
  7. C contribution matricielle provenant de la convection dans le cas de la
  8. C résolution d'une équation de diffusion-convection par EFMH avec le
  9. C modèle DARCY.
  10. C-----------------------------------------------------------------------
  11. C
  12. C---------------------------
  13. C Phrase d'appel (GIBIANE) :
  14. C---------------------------
  15. C
  16. C (RIG2) (CHP3) = 'SMTP' MMODEL TABLE1 RIG1 (TABLE2) (CHP1) (CHP2) ;
  17. C
  18. C------------------------
  19. C Operandes et resultat :
  20. C------------------------
  21. C
  22. C MMODEL : Objet modele décrivant la formulation.
  23. C TABLE1 : TABLE DOMAINE contenant les maillages et les connectivités.
  24. C RIG1 : Matrices hybrides elementaires de DARCY crees par MHYB.
  25. C TABLE2 : TABLE DARCY_TRANSITOIRE contenant les infos transitoires.
  26. C CHP1 : CHPO centre de composante SOUR, source moyenne par élément.
  27. C CHP2 : CHPO face de composante FLUX, flux de vitesse convective.
  28. C CHP3 : CHPO face de composante FLUX, second membre du pb en TH.
  29. C RIG2 : RIGIDITE, contribution du terme convectif au LHS en TH.
  30. C
  31. C-----------------------------------------------------------------------
  32. C
  33. C Langage : ESOPE + FORTRAN77
  34. C
  35. C Auteurs : 09/93 F.DABBENE - Cas permanent
  36. C 09/94 X.NOUVELLON - Extension au transitoire
  37. C 12/95 F.DABBENE - Extension à la convection
  38. C
  39. C-----------------------------------------------------------------------
  40. IMPLICIT INTEGER(I-N)
  41. IMPLICIT REAL*8 (A-H,O-Z)
  42. *
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC SMCHAML
  47. -INC SMCHPOI
  48. -INC SMELEME
  49. -INC SMMODEL
  50. -INC SMRIGID
  51. -INC SMTABLE
  52. *
  53. SEGMENT IPMAHY
  54. INTEGER MAHYBR(NSOUS)
  55. ENDSEGMENT
  56. *
  57. LOGICAL LOGRE,LOGIN
  58. CHARACTER*4 NOMTOT(1)
  59. CHARACTER*6 CHAR6
  60. CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  61. *
  62. * Initialisations
  63. *
  64. ISOU1 = 0
  65. IHN1 = 0
  66. ITPN1 = 0
  67. IVALIN = 0
  68. XVALIN = 0.D0
  69. LOGIN = .TRUE.
  70. IOBIN = 0
  71. TAPIND = 'MOT '
  72. *
  73. * Lecture du MMODEL
  74. *
  75. CALL LIROBJ('MMODEL ',IPMODE,1,IRET)
  76. CALL ACTOBJ('MMODEL ',IPMODE,1)
  77. IF (IERR.NE.0) RETURN
  78. MMODEL = IPMODE
  79. *
  80. * Lecture de la TABLE domaine
  81. *
  82. IPTABL = 0
  83. CALL LEKMOD(IPMODE,IPTABL,IRET)
  84. IF (IERR.NE.0) RETURN
  85. CHARIN = 'MAILLAGE'
  86. TYPOBJ = 'MAILLAGE'
  87. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  88. IF (IERR.NE.0) RETURN
  89. IPGEOM = IOBRE
  90. CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  91. IF (IERR.NE.0) RETURN
  92. IELTFA = IOBRE
  93. CALL LEKTAB(IPTABL,'CENTRE',IOBRE)
  94. IF (IERR.NE.0) RETURN
  95. ICENTR = IOBRE
  96. CALL LEKTAB(IPTABL,'FACE',IOBRE)
  97. IF (IERR.NE.0) RETURN
  98. IPFACE = IOBRE
  99. *
  100. * Lecture de RIGIDITE
  101. *
  102. CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  103. IF (IERR.NE.0) RETURN
  104. MRIGID = IPRIGI
  105. *
  106. * Lecture éventuelle de la TABLE Darcy_transitoire
  107. *
  108. IPTAB1 = 0
  109. CALL LIRTAB('DARCY_TRANSITOIRE',IPTAB1,0,IRET)
  110. IF (IERR.NE.0) RETURN
  111. IF (IRET.EQ.1) THEN
  112. THEMIN = -1.D-12
  113. THEMAX = 1.D0 - THEMIN
  114. TYPOBJ = 'FLOTTANT'
  115. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  116. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  117. IF (IERR.NE.0) RETURN
  118. THETA = XVALRE
  119. IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  120. REAERR(1) = REAL(THETA)
  121. C REAERR(2) = REAL(0.D0)
  122. C REAERR(3) = REAL(1.D0)
  123. MOTERR(1:8) = ' THETA '
  124. CALL ERREUR(42)
  125. RETURN
  126. ENDIF
  127. IF (THETA.LT.0.D0) THETA=0.D0
  128. IF (THETA.GT.1.D0) THETA=1.D0
  129. TYPOBJ = ' '
  130. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA_CONVECTION',
  131. S LOGIN,IOBIN,
  132. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  133. IF (IERR.NE.0) RETURN
  134. IF (TYPOBJ.EQ.'FLOTTANT') THEN
  135. THETAC = XVALRE
  136. IF ((THETAC.LT.THEMIN).OR.(THETAC.GT.THEMAX)) THEN
  137. REAERR(1) = REAL(THETAC)
  138. REAERR(2) = REAL(0.D0)
  139. REAERR(3) = REAL(1.D0)
  140. MOTERR(1:8) = 'TETACONV'
  141. CALL ERREUR(42)
  142. RETURN
  143. ENDIF
  144. IF (THETAC.LT.0.D0) THETA=0.D0
  145. IF (THETAC.GT.1.D0) THETA=1.D0
  146. ELSE
  147. THETAC = THETA
  148. ENDIF
  149. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  150. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  151. IF (IERR.NE.0) RETURN
  152. DELTAT = XVALRE
  153. IF (DELTAT.LE.0.D0) THEN
  154. REAERR(1) = REAL(DELTAT)
  155. REAERR(2) = REAL(0.D0)
  156. MOTERR(1:8) = ' DELTAT '
  157. CALL ERREUR(41)
  158. RETURN
  159. ENDIF
  160. TYPOBJ = 'MCHAML '
  161. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  162. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  163. CALL ACTOBJ(TYPOBJ,IOBRE,1)
  164. IF (IERR.NE.0) RETURN
  165. IPCK = IOBRE
  166. TYPOBJ = 'CHPOINT '
  167. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN,
  168. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  169. CALL ACTOBJ(TYPOBJ,IOBRE,1)
  170. IF (IERR.NE.0) RETURN
  171. ITP = IOBRE
  172. CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN,
  173. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  174. IF (IERR.NE.0) RETURN
  175. IH = IOBRE
  176. ELSE
  177. THETA = 0.D0
  178. DELTAT = 0.D0
  179. IPCK = 0
  180. ITP = 0
  181. IH = 0
  182. ENDIF
  183. *
  184. * Lecture des champoints facultatifs SOUR et FLUX
  185. * prise en compte des termes sources
  186. * prise en compte de la convection
  187. *
  188. ISOUR = 0
  189. CALL LIROBJ('CHPOINT ',ISOUR,0,IRET)
  190. IF(IRET .EQ. 1) CALL ACTOBJ('CHPOINT ',ISOUR,1)
  191. IF (IERR.NE.0) RETURN
  192. *
  193. IFLUX = 0
  194. CALL LIROBJ('CHPOINT ',IFLUX,0,IRET)
  195. IF(IRET .EQ. 1) CALL ACTOBJ('CHPOINT ',IFLUX,1)
  196. IF (IERR.NE.0) RETURN
  197. *
  198. * Tri des deux champoints facultatifs
  199. *
  200. IFAC = ISOUR + IFLUX
  201. IF (IFAC.EQ.0.AND.IPTAB1.EQ.0) THEN
  202. MOTERR(1:8) = 'CHAMPOIN'
  203. CALL ERREUR(37)
  204. RETURN
  205. ENDIF
  206. IF (IFAC.NE.0) THEN
  207. NBCOMP = 0
  208. INDIC = 0
  209. CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
  210. IF (NOMTOT(1).NE.'SOUR'.AND.NOMTOT(1).NE.'FLUX') THEN
  211. MOTERR(1:4) = NOMTOT(1)
  212. CALL ERREUR(197)
  213. RETURN
  214. ELSE IF (NOMTOT(1).EQ.'FLUX') THEN
  215. IFAC = IFLUX
  216. IFLUX = ISOUR
  217. ISOUR = IFAC
  218. ENDIF
  219. ENDIF
  220. *
  221. * Dans le cas de la convection, récup. MCHAML orientation normale
  222. *
  223. IF (IFLUX.NE.0) THEN
  224. CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  225. IF (IERR.NE.0) RETURN
  226. ENDIF
  227. *
  228. * Test du CHAMPOINT de composante SOUR de spg CENTRE
  229. *
  230. IF (ISOUR.NE.0) THEN
  231. NBCOMP = 1
  232. INDIC = 1
  233. NOMTOT(1) = 'SOUR'
  234. CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
  235. IF (IERR.NE.0) RETURN
  236. ENDIF
  237. *
  238. * Test du CHAMPOINT de composante FLUX de spg FACE
  239. *
  240. IF (IFLUX.NE.0) THEN
  241. NBCOMP = 1
  242. INDIC = 1
  243. NOMTOT(1) = 'FLUX'
  244. CALL QUEPOI(IFLUX,IPFACE,INDIC,NBCOMP,NOMTOT)
  245. IF (IERR.NE.0) RETURN
  246. ENDIF
  247. *
  248. * Test du CHAMPOINT de composante TH de spg FACES
  249. *
  250. IF (ITP.NE.0) THEN
  251. NBCOMP = 1
  252. INDIC = 1
  253. NOMTOT(1) = 'TH'
  254. CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT)
  255. IF (IERR.NE.0) RETURN
  256. ENDIF
  257. *
  258. * Test du CHAMPOINT de composante H de spg CENTRE
  259. *
  260. IF (IH.NE.0) THEN
  261. NBCOMP = 1
  262. INDIC = 1
  263. NOMTOT(1) = 'H'
  264. CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT)
  265. IF (IERR.NE.0) RETURN
  266. ENDIF
  267. *
  268. * Test de la formulation
  269. *
  270. SEGACT MMODEL
  271. NSOUS = KMODEL(/1)
  272. SEGINI IPMAHY
  273. IDARCY = 0
  274. IPT1=IPGEOM
  275. IPT2=IPGEOM
  276. DO 10 ISOUS=1,NSOUS
  277. IF(NSOUS.GT.1)THEN
  278. SEGACT IPT2
  279. IPT1=IPT2.LISOUS(ISOUS)
  280. ENDIF
  281. IMODEL = KMODEL(ISOUS)
  282. SEGACT IMODEL
  283. LETYPE = FORMOD(1)
  284. IF (LETYPE.EQ.'DARCY') THEN
  285. IDARCY = IDARCY + 1
  286. MAHYBR(ISOUS) = IPT1
  287. ENDIF
  288. 10 CONTINUE
  289. IF (IDARCY.EQ.0) THEN
  290. MOTERR = LETYPE
  291. CALL ERREUR(193)
  292. GOTO 100
  293. ENDIF
  294. *
  295. * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  296. *
  297. MELEME = IELTFA
  298. SEGACT MELEME
  299. LZONES = LISOUS(/1)
  300. IF (LZONES.EQ.0) LZONES = 1
  301. IPT1 = IPGEOM
  302. SEGACT IPT1
  303. DO 30 ISOUS=1,NSOUS
  304. IMAMEL = MAHYBR(ISOUS)
  305. IF (IMAMEL.NE.0) THEN
  306. DO 20 ISZ=1,LZONES
  307. IF (LZONES.EQ.1) THEN
  308. IPT2 = IPT1
  309. IPT3 = MELEME
  310. ELSE
  311. IPT2 = IPT1.LISOUS(ISZ)
  312. IPT3 = LISOUS(ISZ)
  313. ENDIF
  314. IF (IPT2.EQ.IMAMEL) THEN
  315. MAHYBR(ISOUS) = IPT3
  316. GOTO 30
  317. ENDIF
  318. 20 CONTINUE
  319. IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  320. INTERR(1) = ISOUS
  321. CALL ERREUR(698)
  322. GOTO 100
  323. ENDIF
  324. ENDIF
  325. 30 CONTINUE
  326. *
  327. * Test du sous-type de la matrice de rigiditée récupérée
  328. *
  329. SEGACT MRIGID
  330. LETYPE = MTYMAT
  331. IF (LETYPE.NE.'DARCY') THEN
  332. MOTERR(1:8) = 'RIGIDITE'
  333. MOTERR(9:16) = 'DARCY '
  334. CALL ERREUR(79)
  335. SEGDES MRIGID
  336. GOTO 100
  337. ENDIF
  338. *
  339. * Controle des pointeurs de MELEME de la rigidité
  340. *
  341. DO 40 ISOUS=1,NSOUS
  342. IMAMEL = MAHYBR(ISOUS)
  343. IF (IMAMEL.NE.0) THEN
  344. IPTEST = IRIGEL(1,ISOUS)
  345. IF (MAHYBR(ISOUS).NE.IRIGEL(1,ISOUS)) THEN
  346. MOTERR(1:8) = 'DARCY '
  347. MOTERR(9:16) = 'ELTFA '
  348. INTERR(1) = ISOUS
  349. CALL ERREUR(698)
  350. SEGDES MRIGID
  351. GOTO 100
  352. ENDIF
  353. ENDIF
  354. 40 CONTINUE
  355. SEGDES MRIGID
  356. *
  357. * Vérification du support du MCHAML % au MMODEL
  358. * Controle du nombre de composantes reelles
  359. *
  360. IF (IPCK.NE.0) THEN
  361. ISUP = 2
  362. ICOND = 1
  363. CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  364. IF (IRET.GT.1) GOTO 100
  365. MCHELM = IPCK
  366. SEGACT MCHELM
  367. DO 50 ISOUS=1,NSOUS
  368. IF (MAHYBR(ISOUS).NE.0) THEN
  369. MCHAML = ICHAML(ISOUS)
  370. SEGACT MCHAML
  371. N2 = TYPCHE(/2)
  372. IF (N2.GT.1) THEN
  373. CALL ERREUR(320)
  374. GOTO 100
  375. ENDIF
  376. CHAR6 = TYPCHE(1)(1:6)
  377. IF (CHAR6.NE.'REAL*8') THEN
  378. MOTERR(1:8) = NOMCHE(1)
  379. CALL ERREUR(679)
  380. GOTO 100
  381. ENDIF
  382. ENDIF
  383. 50 CONTINUE
  384. ENDIF
  385. *
  386. * Construction de MCHPOI
  387. *
  388. SEGDES IPMAHY
  389. CALL SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK,
  390. S ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1)
  391. *
  392. IF (MCHPOI.NE.0) THEN
  393. CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  394. CALL ECROBJ('CHPOINT ',MCHPOI)
  395. ENDIF
  396. IF (IRI1.NE.0) CALL ECROBJ('RIGIDITE',IRI1)
  397. *
  398. * Ménage
  399. *
  400. 100 CONTINUE
  401. SEGSUP IPMAHY
  402. END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  

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