Télécharger hybp1.eso

Retour à la liste

Numérotation des lignes :

  1. C HYBP1 SOURCE CB215821 19/08/20 21:18:21 10287
  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. ITES = IK
  130. *
  131. *= Extraction des TH et réduction sur le support IPT1 du CHPO ITPN
  132. *= Les valeurs du CHPOINT en TH sont dans VAA avec l'ordre IPT1
  133. *
  134. SEGINI IORGA
  135. SEGACT MCHPOI
  136. NSOUPO = IPCHP(/1)
  137. DO 30 ISOU=1,NSOUPO
  138. MSOUPO = IPCHP(ISOU)
  139. SEGACT MSOUPO
  140. NC = NOCOMP(/2)
  141. NUMHAR = NOHARM(1)
  142. ITTP = 0
  143. DO 20 I=1,NC
  144. IF (NOCOMP(I).EQ.MOREFP) THEN
  145. ITTP = I
  146. ENDIF
  147. 20 CONTINUE
  148. IF (ITTP.NE.0) THEN
  149. MELEME = IGEOC
  150. SEGACT MELEME
  151. N2 = NUM(/2)
  152. MPOVAL = IPOVAL
  153. SEGACT MPOVAL
  154. DO 25 I2=1,N2
  155. K = NUM(1,I2)
  156. IF (ICPR(K).EQ.0) THEN
  157. INTERR(1) = K
  158. MOTERR(1:8) = 'FACE '
  159. CALL ERREUR(64)
  160. SEGDES MRIGID, IPMAHY
  161. SEGSUP ICCPR, IORGA
  162. RETURN
  163. ELSE
  164. VAA(ICPR(K)) = VPOCHA(I2,ITTP)
  165. ENDIF
  166. 25 CONTINUE
  167. ENDIF
  168. 30 CONTINUE
  169. *
  170. *= Création des chapeaux du CHPO centre P résultat
  171. *
  172. NSOUPO = 1
  173. NAT = 1
  174. SEGINI MCHPOI
  175. IRET = MCHPOI
  176. DO 31 ITIT=1,72
  177. MOCHDE(ITIT:ITIT)=BLAN1
  178. 31 CONTINUE
  179. MTYPOI = 'CENTRE '
  180. IFOPOI = IFOMOD
  181. JATTRI(1) = 1
  182. NC = 1
  183. SEGINI MSOUPO
  184. IPCHP(1) = MSOUPO
  185. NOHARM(1) = NUMHAR
  186. IGEOC = ICENTR
  187. NOCOMP(1) = 'H '
  188. IPT2 = ICENTR
  189. SEGACT IPT2
  190. N = IPT2.NUM(/2)
  191. SEGINI MPOVAL
  192. IPOVAL = MPOVAL
  193. *
  194. * Activation du MPOVAL du CHPO centre SOUR
  195. *
  196. IF (ISOUR.NE.0) THEN
  197. SEGACT MCHPO1
  198. MSOUP1 = MCHPO1.IPCHP(1)
  199. SEGACT MSOUP1
  200. MPOVA1 = MSOUP1.IPOVAL
  201. SEGACT MPOVA1
  202. ENDIF
  203. *
  204. * Activation du MPOVAL du CHPO face trace de charge au temps (n-1)
  205. *
  206. IF (ITP.NE.0) THEN
  207. SEGACT MCHPO2
  208. MSOUP2 = MCHPO2.IPCHP(1)
  209. SEGACT MSOUP2
  210. MPOVA2 = MSOUP2.IPOVAL
  211. SEGACT MPOVA2
  212. ENDIF
  213. *
  214. * Activation du MPOVAL du CHPO centre charge au temps (n-1)
  215. *
  216. IF (IH.NE.0) THEN
  217. SEGACT MCHPO3
  218. MSOUP3 = MCHPO3.IPCHP(1)
  219. SEGACT MSOUP3
  220. MPOVA3 = MSOUP3.IPOVAL
  221. SEGACT MPOVA3
  222. ENDIF
  223. *
  224. * activation des objets liés à la présence d'une force volumique
  225. *
  226. IF (IFORC.NE.0) THEN
  227. *
  228. * Activation du MPOVAL du CHPO force appuyé au centre des éléments volumiques
  229. *
  230. SEGACT MCHPO4
  231. MSOUP4 = MCHPO4.IPCHP(1)
  232. SEGACT MSOUP4
  233. MPOVA4 = MSOUP4.IPOVAL
  234. SEGACT MPOVA4
  235. *
  236. * Activation du MPOVAL du CHPO des vecteurs normales appuyé au centre des faces
  237. *
  238. SEGACT MCHPO5
  239. MSOUP5 = MCHPO5.IPCHP(1)
  240. SEGACT MSOUP5
  241. MPOVA5 = MSOUP5.IPOVAL
  242. SEGACT MPOVA5
  243. *
  244. * Activation du MPOVAL du CHPO des surfaces appuyé au centre des faces
  245. *
  246. SEGACT MCHPO6
  247. MSOUP6 = MCHPO6.IPCHP(1)
  248. SEGACT MSOUP6
  249. MPOVA6 = MSOUP6.IPOVAL
  250. SEGACT MPOVA6
  251. *
  252. * Activation du MCHEL des orientations des faces
  253. *
  254. SEGACT MCHEL1
  255. *
  256. * Activation du MRIGI de la matrice masse hybride
  257. *
  258. SEGACT RI1
  259. ENDIF
  260. *
  261. * Activation de Ck|K|
  262. *
  263. IF (IPCK.NE.0) SEGACT MCHELM
  264. *
  265. *---------------------------------------
  266. *= Boucle sur les maillages elementaires
  267. *---------------------------------------
  268. *
  269. ITELEM = 0
  270. DO 110 IMAIL=1,NBMAIL
  271. IF (MAHYBR(IMAIL).EQ.0) GOTO 110
  272. *
  273. * Récupération des infos pour la zone IMAIL ou Darcy est définie
  274. *
  275. MELEME = IRIGEL(1,IMAIL)
  276. xMATRI = IRIGEL(4,IMAIL)
  277. SEGACT MELEME
  278. NBDDL = NUM(/1)
  279. NBELEM = NUM(/2)
  280. SEGACT xMATRI
  281. SEGINI MTRAV
  282. *
  283. * Activation des objets necessaires à la prise en compte des forces de volumes
  284. *
  285. IF (IFORC.NE.0) THEN
  286. MCHAM1 = MCHEL1.ICHAML(IMAIL)
  287. SEGACT MCHAM1
  288. MELVA1 = MCHAM1.IELVAL(1)
  289. SEGACT MELVA1
  290. xMATR1 = RI1.IRIGEL(4,IMAIL)
  291. SEGACT xMATR1
  292. ELSE
  293. DO 35 I=1,NBDDL
  294. RFOR(I)=0.D0
  295. 35 CONTINUE
  296. ENDIF
  297. *
  298. C
  299. C Activation du MELVAL du MCHAML Ck|K|
  300. C
  301. IF (IPCK.NE.0) THEN
  302. MCHAML = ICHAML(IMAIL)
  303. SEGACT MCHAML
  304. MELVAL = IELVAL(1)
  305. SEGACT MELVAL
  306. ENDIF
  307. *
  308. *-------------------------------------------------------
  309. * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  310. *-------------------------------------------------------
  311. *
  312. DO 100 IEL=1,NBELEM
  313. ITELEM = ITELEM + 1
  314. *
  315. *- Calcul de VECT = DIV * RE , CONSD = VECT * DIV
  316. *
  317. * XMATRI = IMATTT(IEL)
  318. * SEGACT XMATRI
  319. CONSD = 0.D0
  320. DO 45 J=1,NBDDL
  321. VECT(J) = 0.D0
  322. DO 40 I=1,NBDDL
  323. VECT(J) = VECT(J) + RE(I,J,iel)
  324. 40 CONTINUE
  325. CONSD = CONSD + VECT(J)
  326. 45 CONTINUE
  327. * SEGDES XMATRI
  328. *
  329. *- Recuperation de SOUR et de H(n-1)
  330. *
  331. IF (ISOUR.NE.0) VSOUR=MPOVA1.VPOCHA(ITELEM,1)
  332. IF (IFORC.NE.0) THEN
  333. *
  334. *- calcul des flux de forces aux faces de l'element
  335. *
  336. DO 55 IDDL=1,NBDDL
  337. FLFOR(IDDL)= 0.D0
  338. IPOPTS = ICPR(NUM(IDDL,IEL))
  339. DO 50 I=1,IDIM
  340. FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I)
  341. $ *MELVA1.VELCHE(IDDL,IEL) * MPOVA4
  342. $ .VPOCHA(ITELEM,I) *MPOVA6.VPOCHA(IPOPTS,1)
  343. 50 CONTINUE
  344. 55 CONTINUE
  345. *
  346. *- Construction du tableau aux faces M.FORCE
  347. *
  348. * XMATR1 = IMATR1.IMATTT(IEL)
  349. * SEGACT XMATR1
  350. DO 65 I=1,NBDDL
  351. RFOR(I)=0.D0
  352. DO 60 J=1,NBDDL
  353. RFOR(I) = RFOR(I) + XMATR1.RE(I,J,iel)*FLFOR(J)
  354. 60 CONTINUE
  355. 65 CONTINUE
  356. * SEGDES XMATR1
  357. ENDIF
  358. *
  359. *- Calcul de p pour l'element IEL
  360. *
  361. VALP = 0.D0
  362. IF (IPCK.EQ.0) THEN
  363. DO 80 I=1,NBDDL
  364. NUMFA = ICPR(NUM(I,IEL))
  365. VALP = VALP + VECT(I) * (VAA(NUMFA) - RFOR(I))
  366. 80 CONTINUE
  367. VPOCHA(ITELEM,1) = (VALP + VSOUR) / CONSD
  368. ELSE
  369. BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT)
  370. DO 90 I=1,NBDDL
  371. NUMFA = ICPR(NUM(I,IEL))
  372. VALP = VALP + VECT(I) * ( THETA*VAA(NUMFA)
  373. S + (1.D0-THETA)*MPOVA2.VPOCHA(NUMFA,1))
  374. 90 CONTINUE
  375. VALP = (VALP + VSOUR) * BETA / CONSD
  376. VH = MPOVA3.VPOCHA(ITELEM,1)
  377. VPOCHA(ITELEM,1) = VALP + (1.D0-BETA)*VH
  378. ENDIF
  379. 100 CONTINUE
  380. IF (IFORC.NE.0) SEGDES xMATR1
  381. SEGDES xMATRI
  382. SEGSUP MTRAV
  383. 110 CONTINUE
  384. SEGDES MRIGID, IPMAHY
  385. IF (IFORC.NE.0) SEGDES RI1
  386. SEGSUP ICCPR , IORGA
  387. *
  388. END
  389.  
  390.  
  391.  

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