Télécharger hybp1.eso

Retour à la liste

Numérotation des lignes :

  1. C HYBP1 SOURCE CHAT 09/10/09 21:18:48 6519
  2. SUBROUTINE HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE,
  3. S ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,IRET)
  4. C-----------------------------------------------------------------------
  5. C Calcul de la charge au centre de chaque maille
  6. C-----------------------------------------------------------------------
  7. C
  8. C---------------------------
  9. C Parametres Entree/Sortie :
  10. C---------------------------
  11. C
  12. C E/ IPMAHY : Segment contenant le pointeur vers le meleme des
  13. C connectivites elements/faces pour les zones du MMODEL
  14. C ou on a defini DARCY.
  15. C E/ ICENTR : Pointeur vers l'objet maillage CENTRE
  16. C E/ IPFACE : Pointeur vers l'objet maillage FACE -1
  17. C E/ IPDARC : Matrice de sous type DARCY (contient RE ).
  18. C E/ ITPN : CHPO face des Traces de charge TH au temps n
  19. C E/ THETA : Parametre de discretisation temporelle (Theta-methode)
  20. C E/ DELTAT : Pas de temps
  21. C E/ IPCK : MCHAML donnant pour chaque element Ck|K|
  22. C E/ ITP : CHPO face des traces de charge TH au temps n-1
  23. C E/ IH : CHPO centre des charges H au temps n-1
  24. C E/ ISOUR : CHPO centre SOUR contenant un eventuel terme source
  25. C E/ IPFORC : Matrice de sous type MASSE (contient RE ).
  26. C E/ IFORC : CHPO centre force contenant la force volumique
  27. C /S IRET : CHPO centre H resultat contenant la charge au temps n
  28. C
  29. C----------------------
  30. C Variables en COMMON :
  31. C----------------------
  32. C
  33. C E/ IDIM : cf CCOPTIO
  34. C E/ IFOMOD : cf CCOPTIO
  35. C E/ NIFOUR : cf CCOPTIO
  36. C
  37. C----------------------
  38. C Tableaux de travail :
  39. C----------------------
  40. C
  41. C ICPR(I)=J : Le noeud I a le numero local J
  42. C Correspondance numerotation globale/locale
  43. C NNGOT : Nombre de noeuds total du domaine
  44. C
  45. C-----------------------------------------------------------------------
  46. C
  47. C Langage : ESOPE + FORTRAN77
  48. C
  49. C Auteurs : 09/93 F.DABBENE - Cas permanent
  50. C 09/94 X.NOUVELLON - Extension au cas transitoire
  51. C 02/96 L.V.BENET - Prise en compte de forces de volume
  52. C
  53. C-----------------------------------------------------------------------
  54. IMPLICIT INTEGER(I-N)
  55. IMPLICIT REAL*8 (A-H,O-Z)
  56. *
  57. -INC CCOPTIO
  58. -INC SMCHAML
  59. -INC SMCHPOI
  60. POINTEUR MCHPO5.MCHPOI,MCHPO6.MCHPOI,MSOUP6.MSOUPO
  61. -INC SMCOORD
  62. -INC SMELEME
  63. -INC SMRIGID
  64. *
  65. CHARACTER*4 MOREFD, MOREFP
  66. CHARACTER*1 BLAN1
  67. SEGMENT ICCPR
  68. INTEGER ICPR(NNGOT)
  69. ENDSEGMENT
  70. SEGMENT IORGA
  71. REAL*8 VAA(ITES)
  72. ENDSEGMENT
  73. SEGMENT IPMAHY
  74. INTEGER MAHYBR(NSOUS)
  75. ENDSEGMENT
  76. SEGMENT MTRAV
  77. REAL*8 VECT(NBDDL),FLFOR(NBDDL),RFOR(NBDDL)
  78. ENDSEGMENT
  79. *
  80. BLAN1=' '
  81. *
  82. * Initialisations
  83. *
  84. IRET = 0
  85. IPT1 = IPFACE
  86. MRIGID = IPDARC
  87. RI1 = IPFORC
  88. MCHPOI = ITPN
  89. MCHELM = IPCK
  90. MCHEL1 = IORIE
  91. MCHPO1 = ISOUR
  92. MCHPO2 = ITP
  93. MCHPO3 = IH
  94. MCHPO4 = IFORC
  95. MCHPO5 = INORM
  96. MCHPO6 = ISURF
  97. VSOUR = 0.D0
  98. *
  99. * Initialisations se servant de IPMAHY et MRIGID
  100. *
  101. SEGACT IPMAHY
  102. NBMAIL = MAHYBR(/1)
  103. SEGACT MRIGID
  104. DO 5 IMAIL=1,NBMAIL
  105. IF (MAHYBR(IMAIL).NE.0) THEN
  106. DESCR = IRIGEL(3,IMAIL)
  107. SEGACT DESCR
  108. MOREFP = LISINC(1)
  109. SEGDES DESCR
  110. GOTO 10
  111. ENDIF
  112. 5 CONTINUE
  113. 10 CONTINUE
  114. NNGOT = XCOOR(/1)/(IDIM+1)
  115. SEGINI ICCPR
  116. *
  117. *= Creation du tableau ICPR pour le maillage IPT1
  118. *
  119. SEGACT IPT1
  120. N2 = IPT1.NUM(/2)
  121. IK = 0
  122. DO 15 I2=1,N2
  123. K = IPT1.NUM(1,I2)
  124. IF (ICPR(K).EQ.0) THEN
  125. IK = IK + 1
  126. ICPR(K) = IK
  127. ENDIF
  128. 15 CONTINUE
  129. SEGDES IPT1
  130. ITES = IK
  131. *
  132. *= Extraction des TH et réduction sur le support IPT1 du CHPO ITPN
  133. *= Les valeurs du CHPOINT en TH sont dans VAA avec l'ordre IPT1
  134. *
  135. SEGINI IORGA
  136. SEGACT MCHPOI
  137. NSOUPO = IPCHP(/1)
  138. DO 30 ISOU=1,NSOUPO
  139. MSOUPO = IPCHP(ISOU)
  140. SEGACT MSOUPO
  141. NC = NOCOMP(/2)
  142. NUMHAR = NOHARM(1)
  143. ITTP = 0
  144. DO 20 I=1,NC
  145. IF (NOCOMP(I).EQ.MOREFP) THEN
  146. ITTP = I
  147. ENDIF
  148. 20 CONTINUE
  149. IF (ITTP.NE.0) THEN
  150. MELEME = IGEOC
  151. SEGACT MELEME
  152. N2 = NUM(/2)
  153. MPOVAL = IPOVAL
  154. SEGACT MPOVAL
  155. DO 25 I2=1,N2
  156. K = NUM(1,I2)
  157. IF (ICPR(K).EQ.0) THEN
  158. INTERR(1) = K
  159. MOTERR(1:8) = 'FACE '
  160. CALL ERREUR(64)
  161. SEGDES MPOVAL, MELEME, MSOUPO
  162. SEGDES MCHPOI, MRIGID, IPMAHY
  163. SEGSUP ICCPR, IORGA
  164. RETURN
  165. ELSE
  166. VAA(ICPR(K)) = VPOCHA(I2,ITTP)
  167. ENDIF
  168. 25 CONTINUE
  169. SEGDES MPOVAL, MELEME
  170. ENDIF
  171. SEGDES MSOUPO
  172. 30 CONTINUE
  173. SEGDES MCHPOI
  174. *
  175. *= Création des chapeaux du CHPO centre P résultat
  176. *
  177. NSOUPO = 1
  178. NAT = 1
  179. SEGINI MCHPOI
  180. IRET = MCHPOI
  181. DO 31 ITIT=1,72
  182. MOCHDE(ITIT:ITIT)=BLAN1
  183. 31 CONTINUE
  184. MTYPOI = 'CENTRE '
  185. IFOPOI = IFOMOD
  186. JATTRI(1) = 1
  187. NC = 1
  188. SEGINI MSOUPO
  189. IPCHP(1) = MSOUPO
  190. SEGDES MCHPOI
  191. NOHARM(1) = NUMHAR
  192. IGEOC = ICENTR
  193. NOCOMP(1) = 'H '
  194. IPT2 = ICENTR
  195. SEGACT IPT2
  196. N = IPT2.NUM(/2)
  197. SEGDES IPT2
  198. SEGINI MPOVAL
  199. IPOVAL = MPOVAL
  200. SEGDES MSOUPO
  201. *
  202. * Activation du MPOVAL du CHPO centre SOUR
  203. *
  204. IF (ISOUR.NE.0) THEN
  205. SEGACT MCHPO1
  206. MSOUP1 = MCHPO1.IPCHP(1)
  207. SEGDES MCHPO1
  208. SEGACT MSOUP1
  209. MPOVA1 = MSOUP1.IPOVAL
  210. SEGDES MSOUP1
  211. SEGACT MPOVA1
  212. ENDIF
  213. *
  214. * Activation du MPOVAL du CHPO face trace de charge au temps (n-1)
  215. *
  216. IF (ITP.NE.0) THEN
  217. SEGACT MCHPO2
  218. MSOUP2 = MCHPO2.IPCHP(1)
  219. SEGDES MCHPO2
  220. SEGACT MSOUP2
  221. MPOVA2 = MSOUP2.IPOVAL
  222. SEGDES MSOUP2
  223. SEGACT MPOVA2
  224. ENDIF
  225. *
  226. * Activation du MPOVAL du CHPO centre charge au temps (n-1)
  227. *
  228. IF (IH.NE.0) THEN
  229. SEGACT MCHPO3
  230. MSOUP3 = MCHPO3.IPCHP(1)
  231. SEGDES MCHPO3
  232. SEGACT MSOUP3
  233. MPOVA3 = MSOUP3.IPOVAL
  234. SEGDES MSOUP3
  235. SEGACT MPOVA3
  236. ENDIF
  237. *
  238. * activation des objets liés à la présence d'une force volumique
  239. *
  240. IF (IFORC.NE.0) THEN
  241. *
  242. * Activation du MPOVAL du CHPO force appuyé au centre des éléments volumiques
  243. *
  244. SEGACT MCHPO4
  245. MSOUP4 = MCHPO4.IPCHP(1)
  246. SEGDES MCHPO4
  247. SEGACT MSOUP4
  248. MPOVA4 = MSOUP4.IPOVAL
  249. SEGDES MSOUP4
  250. SEGACT MPOVA4
  251. *
  252. * Activation du MPOVAL du CHPO des vecteurs normales appuyé au centre des faces
  253. *
  254. SEGACT MCHPO5
  255. MSOUP5 = MCHPO5.IPCHP(1)
  256. SEGDES MCHPO5
  257. SEGACT MSOUP5
  258. MPOVA5 = MSOUP5.IPOVAL
  259. SEGDES MSOUP5
  260. SEGACT MPOVA5
  261. *
  262. * Activation du MPOVAL du CHPO des surfaces appuyé au centre des faces
  263. *
  264. SEGACT MCHPO6
  265. MSOUP6 = MCHPO6.IPCHP(1)
  266. SEGDES MCHPO6
  267. SEGACT MSOUP6
  268. MPOVA6 = MSOUP6.IPOVAL
  269. SEGDES MSOUP6
  270. SEGACT MPOVA6
  271. *
  272. * Activation du MCHEL des orientations des faces
  273. *
  274. SEGACT MCHEL1
  275. *
  276. * Activation du MRIGI de la matrice masse hybride
  277. *
  278. SEGACT RI1
  279. ENDIF
  280. *
  281. * Activation de Ck|K|
  282. *
  283. IF (IPCK.NE.0) SEGACT MCHELM
  284. *
  285. *---------------------------------------
  286. *= Boucle sur les maillages elementaires
  287. *---------------------------------------
  288. *
  289. ITELEM = 0
  290. DO 110 IMAIL=1,NBMAIL
  291. IF (MAHYBR(IMAIL).EQ.0) GOTO 110
  292. *
  293. * Récupération des infos pour la zone IMAIL ou Darcy est définie
  294. *
  295. MELEME = IRIGEL(1,IMAIL)
  296. xMATRI = IRIGEL(4,IMAIL)
  297. SEGACT MELEME
  298. NBDDL = NUM(/1)
  299. NBELEM = NUM(/2)
  300. SEGACT xMATRI
  301. SEGINI MTRAV
  302. *
  303. * Activation des objets necessaires à la prise en compte des forces de volumes
  304. *
  305. IF (IFORC.NE.0) THEN
  306. MCHAM1 = MCHEL1.ICHAML(IMAIL)
  307. SEGACT MCHAM1
  308. MELVA1 = MCHAM1.IELVAL(1)
  309. SEGDES MCHAM1
  310. SEGACT MELVA1
  311. xMATR1 = RI1.IRIGEL(4,IMAIL)
  312. SEGACT xMATR1
  313. ELSE
  314. DO 35 I=1,NBDDL
  315. RFOR(I)=0.D0
  316. 35 CONTINUE
  317. ENDIF
  318. *
  319. C
  320. C Activation du MELVAL du MCHAML Ck|K|
  321. C
  322. IF (IPCK.NE.0) THEN
  323. MCHAML = ICHAML(IMAIL)
  324. SEGACT MCHAML
  325. MELVAL = IELVAL(1)
  326. SEGDES MCHAML
  327. SEGACT MELVAL
  328. ENDIF
  329. *
  330. *-------------------------------------------------------
  331. * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  332. *-------------------------------------------------------
  333. *
  334. DO 100 IEL=1,NBELEM
  335. ITELEM = ITELEM + 1
  336. *
  337. *- Calcul de VECT = DIV * RE , CONSD = VECT * DIV
  338. *
  339. * XMATRI = IMATTT(IEL)
  340. * SEGACT XMATRI
  341. CONSD = 0.D0
  342. DO 45 J=1,NBDDL
  343. VECT(J) = 0.D0
  344. DO 40 I=1,NBDDL
  345. VECT(J) = VECT(J) + RE(I,J,iel)
  346. 40 CONTINUE
  347. CONSD = CONSD + VECT(J)
  348. 45 CONTINUE
  349. * SEGDES XMATRI
  350. *
  351. *- Recuperation de SOUR et de H(n-1)
  352. *
  353. IF (ISOUR.NE.0) VSOUR=MPOVA1.VPOCHA(ITELEM,1)
  354. IF (IFORC.NE.0) THEN
  355. *
  356. *- calcul des flux de forces aux faces de l'element
  357. *
  358. DO 55 IDDL=1,NBDDL
  359. FLFOR(IDDL)= 0.D0
  360. IPOPTS = ICPR(NUM(IDDL,IEL))
  361. DO 50 I=1,IDIM
  362. FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I)
  363. $ *MELVA1.VELCHE(IDDL,IEL) * MPOVA4
  364. $ .VPOCHA(ITELEM,I) *MPOVA6.VPOCHA(IPOPTS,1)
  365. 50 CONTINUE
  366. 55 CONTINUE
  367. *
  368. *- Construction du tableau aux faces M.FORCE
  369. *
  370. * XMATR1 = IMATR1.IMATTT(IEL)
  371. * SEGACT XMATR1
  372. DO 65 I=1,NBDDL
  373. RFOR(I)=0.D0
  374. DO 60 J=1,NBDDL
  375. RFOR(I) = RFOR(I) + XMATR1.RE(I,J,iel)*FLFOR(J)
  376. 60 CONTINUE
  377. 65 CONTINUE
  378. * SEGDES XMATR1
  379. ENDIF
  380. *
  381. *- Calcul de p pour l'element IEL
  382. *
  383. VALP = 0.D0
  384. IF (IPCK.EQ.0) THEN
  385. DO 80 I=1,NBDDL
  386. NUMFA = ICPR(NUM(I,IEL))
  387. VALP = VALP + VECT(I) * (VAA(NUMFA) - RFOR(I))
  388. 80 CONTINUE
  389. VPOCHA(ITELEM,1) = (VALP + VSOUR) / CONSD
  390. ELSE
  391. BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT)
  392. DO 90 I=1,NBDDL
  393. NUMFA = ICPR(NUM(I,IEL))
  394. VALP = VALP + VECT(I) * ( THETA*VAA(NUMFA)
  395. S + (1.D0-THETA)*MPOVA2.VPOCHA(NUMFA,1))
  396. 90 CONTINUE
  397. VALP = (VALP + VSOUR) * BETA / CONSD
  398. VH = MPOVA3.VPOCHA(ITELEM,1)
  399. VPOCHA(ITELEM,1) = VALP + (1.D0-BETA)*VH
  400. ENDIF
  401. 100 CONTINUE
  402. IF (IFORC.NE.0) SEGDES xMATR1, MELVA1
  403. SEGDES xMATRI, MELEME
  404. IF (IPCK.NE.0) SEGDES MELVAL
  405. SEGSUP MTRAV
  406. 110 CONTINUE
  407. SEGDES MPOVAL
  408. SEGDES MRIGID, IPMAHY
  409. IF (IPCK.NE.0) SEGDES MCHELM, MPOVA2, MPOVA3
  410. IF (ISOUR.NE.0) SEGDES MPOVA1
  411. IF (IFORC.NE.0)
  412. S SEGDES MPOVA4, MPOVA5, MPOVA6,
  413. S MCHEL1, RI1
  414. SEGSUP ICCPR , IORGA
  415. *
  416. RETURN
  417. END
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  

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