Télécharger hybp1.eso

Retour à la liste

Numérotation des lignes :

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

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