Télécharger hybp.eso

Retour à la liste

Numérotation des lignes :

  1. C HYBP SOURCE CHAT 11/03/16 21:25:04 6902
  2. SUBROUTINE HYBP
  3. C-----------------------------------------------------------------------
  4. C Cet operateur permet de calculer la charge en fonction des traces de
  5. C charge dans le cadre d'une formulation variationnelle mixte hybride
  6. C pour les equations de DARCY.
  7. C
  8. C (n) (n-1) -1 t -1
  9. C H = (1 - BETA) * H + (DIV * M1 * DIV) * BETA *
  10. C -1 (n) (n-1)
  11. C ( DIV * M1 * (THETA * TH + (1-THETA) * TH - M * F) + SOUR )
  12. C
  13. C---------------------------
  14. C Phrase d'appel (GIBIANE) :
  15. C---------------------------
  16. C
  17. C CHPO1 = HYBP MMODEL MRIGI1 CHPO2 (TABLE2) (CHPO3) (MRIGI2 CHPO4);
  18. C
  19. C------------------------
  20. C Operandes et resultat :
  21. C------------------------
  22. C
  23. C CHPO1 : CHPO centre de nom de composante H.
  24. C MRIGI1 : RIGIDITE masses hybrides élémentaires de sous-type DARCY.
  25. C MRIGI2 : Matrices masses hybrides elementaires de sous-type MASSE.
  26. C MMODEL : MODELE spécifiant la formulation.
  27. C CHPO2 : CHPO face de nom de composante TH.
  28. C TABLE2 : TABLE DARCY_TRANSITOIRE a inclure en transitoire; contient
  29. C Theta, le pas de temps, Ck|K|, TH(n-1) et H(n-1).
  30. C CHPO3 : CHPO centre de nom de composante SOUR.
  31. C CHPO4 : CHPO centre de composante FX,FY(,FZ),
  32. C densite de force moyenne par élément.
  33. C
  34. C-----------------------------------------------------------------------
  35. C
  36. C Langage : ESOPE + FORTRAN77
  37. C
  38. C Auteurs : 09/93 F.DABBENE - Cas permanent
  39. C 09/94 X.NOUVELLON - Extension au cas transitoire
  40. C 02/96 L.V.BENET - Prise en compte de forces de volume
  41. C
  42. C-----------------------------------------------------------------------
  43. IMPLICIT INTEGER(I-N)
  44. IMPLICIT REAL*8 (A-H,O-Z)
  45. *
  46. -INC CCOPTIO
  47. -INC SMCHAML
  48. -INC SMCHPOI
  49. -INC SMELEME
  50. -INC SMMODEL
  51. -INC SMRIGID
  52. -INC SMTABLE
  53. *
  54. SEGMENT IPMAHY
  55. INTEGER MAHYBR(NSOUS)
  56. ENDSEGMENT
  57. *
  58. LOGICAL LOGRE,LOGIN
  59. CHARACTER*4 NOMTOT(1)
  60. CHARACTER*4 NOMTO3(3)
  61. CHARACTER*6 CHAR6
  62. CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  63. *
  64. * Initialisations
  65. *
  66. IVALIN = 0
  67. XVALIN = 0.D0
  68. LOGIN = .TRUE.
  69. IOBIN = 0
  70. TAPIND = 'MOT '
  71. *
  72. * Lecture du MMODEL et test de la formulation
  73. *
  74. CALL LIROBJ('MMODEL',IPMODE,1,IRET)
  75. IF (IERR.NE.0) RETURN
  76. MMODEL = IPMODE
  77. *
  78. * Lecture de la TABLE domaine
  79. *
  80. IPTABL = 0
  81. C CALL LIRTAB('DOMAINE',IPTABL,1,IRET)
  82. CALL LEKMOD(MMODEL,IPTABL,IRET)
  83. IF (IERR.NE.0) RETURN
  84. CHARIN = 'MAILLAGE'
  85. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  86. IF (IERR.NE.0) RETURN
  87. IPGEOM = IOBRE
  88. CHARIN = 'ELTFA'
  89. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  90. IF (IERR.NE.0) RETURN
  91. IELTFA = IOBRE
  92. CHARIN = 'CENTRE'
  93. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  94. IF (IERR.NE.0) RETURN
  95. ICENTR = IOBRE
  96. CHARIN = 'FACE'
  97. CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  98. IF (IERR.NE.0) RETURN
  99. IPFACE = IOBRE
  100. *
  101. * récup. MCHAML orientation normale
  102. *
  103. CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  104. IF (IERR.NE.0) RETURN
  105. *
  106. * récup. CHPO orientation normale
  107. *
  108. CALL LEKTAB(IPTABL,'XXNORMAF',INORM)
  109. IF (IERR.NE.0) RETURN
  110. *
  111. * récup. CHPO surface des faces
  112. *
  113. CALL LEKTAB(IPTABL,'XXSURFAC',ISURF)
  114. IF (IERR.NE.0) RETURN
  115. *
  116. * Lecture de RIGIDITE obligatoire
  117. *
  118. CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  119. IF (IERR.NE.0) RETURN
  120. RI1 = IPRIGI
  121. *
  122. * Lecture du CHAMPOINT obligatoire
  123. *
  124. CALL LIROBJ('CHPOINT',ITPN,1,IRET)
  125. IF (IERR.NE.0) RETURN
  126. *
  127. * Lecture eventuelle du 1er CHAMPOINT optionnel
  128. *
  129. ICHP1 = 0
  130. CALL LIROBJ('CHPOINT',ICHP1,0,IRET1)
  131. IF (IERR.NE.0) RETURN
  132. *
  133. * Lecture eventuelle du 2me CHAMPOINT optionnel
  134. *
  135. ICHP2 = 0
  136. CALL LIROBJ('CHPOINT',ICHP2,0,IRET2)
  137. IF (IERR.NE.0) RETURN
  138. *
  139. * Lecture de la RIGIDITE complementaire
  140. *
  141. IPRIGI = 0
  142. CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR)
  143. IF (IERR.NE.0) RETURN
  144. RI2 = IPRIGI
  145. *
  146. * Lecture éventuelle de la TABLE Darcy_transitoire
  147. *
  148. IPTABL = 0
  149. CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET)
  150. IF (IERR.NE.0) RETURN
  151. IF (IRET.EQ.1) THEN
  152. TYPOBJ = 'FLOTTANT'
  153. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  154. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  155. IF (IERR.NE.0) RETURN
  156. THETA = XVALRE
  157. THEMIN = -1.D-12
  158. THEMAX = 1.D0 - THEMIN
  159. IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  160. REAERR(1) = REAL(THETA)
  161. REAERR(2) = REAL(0.D0)
  162. REAERR(3) = REAL(1.D0)
  163. MOTERR(1:8) = ' THETA '
  164. CALL ERREUR(42)
  165. RETURN
  166. ENDIF
  167. IF (THETA.LT.0.D0) THETA=0.D0
  168. IF (THETA.GT.1.D0) THETA=1.D0
  169. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  170. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  171. IF (IERR.NE.0) RETURN
  172. DELTAT = XVALRE
  173. IF (DELTAT.LT.0.D0) THEN
  174. REAERR(1) = REAL(DELTAT)
  175. REAERR(2) = REAL(0.D0)
  176. MOTERR(1:8) = ' DELTAT '
  177. CALL ERREUR(41)
  178. RETURN
  179. ENDIF
  180. TYPOBJ = 'MCHAML '
  181. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  182. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  183. IF (IERR.NE.0) RETURN
  184. IPCK = IOBRE
  185. TYPOBJ = 'CHPOINT '
  186. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN,
  187. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  188. IF (IERR.NE.0) RETURN
  189. ITP = IOBRE
  190. CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN,
  191. S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  192. IF (IERR.NE.0) RETURN
  193. IH = IOBRE
  194. ELSE
  195. THETA = 0.D0
  196. DELTAT = 0.D0
  197. IPCK = 0
  198. ITP = 0
  199. IH = 0
  200. ENDIF
  201. *
  202. * Test CHAMPOINT centre H
  203. *
  204. IF (IH.NE.0) THEN
  205. NBCOMP = 1
  206. INDIC = 1
  207. NOMTOT(1) = 'H'
  208. CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT)
  209. IF (IERR.NE.0) RETURN
  210. ENDIF
  211. * (n-1)
  212. * Test du CHAMPOINT face TH
  213. *
  214. IF (ITP.NE.0) THEN
  215. NBCOMP = 1
  216. INDIC = 1
  217. NOMTOT(1) = 'TH'
  218. CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT)
  219. IF (IERR.NE.0) RETURN
  220. ENDIF
  221. *
  222. * Test de la formulation contenu dans MMODEL
  223. *
  224. SEGACT MMODEL
  225. NSOUS = KMODEL(/1)
  226. SEGINI IPMAHY
  227. IDARCY = 0
  228. IPT1=IPGEOM
  229. IPT2=IPGEOM
  230. DO 10 ISOUS=1,NSOUS
  231. IF(NSOUS.GT.1)THEN
  232. SEGACT IPT2
  233. IPT1= IPT2.LISOUS(ISOUS)
  234. SEGDES IPT2
  235. ENDIF
  236. IMODEL = KMODEL(ISOUS)
  237. SEGACT IMODEL
  238. LETYPE = FORMOD(1)
  239. IF (LETYPE.EQ.'DARCY') THEN
  240. IDARCY = IDARCY + 1
  241. MAHYBR(ISOUS) = IPT1
  242. ENDIF
  243. SEGDES IMODEL
  244. 10 CONTINUE
  245. SEGDES MMODEL
  246. IF (IDARCY.EQ.0) THEN
  247. MOTERR = LETYPE
  248. CALL ERREUR(193)
  249. GOTO 100
  250. ENDIF
  251. *
  252. * Tri et test des CHAMPOINTs ARGUMENTS
  253. *
  254. ISOUR=0
  255. IFORC=0
  256. IF(IRET1.EQ.1)THEN
  257. NBCOMP = 0
  258. INDIC = 1
  259. NOMTO3(1) = ' '
  260. NOMTO3(2) = ' '
  261. NOMTO3(3) = ' '
  262. CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTO3)
  263. IF (IERR.NE.0) RETURN
  264. IF(NOMTO3(1).EQ.'SOUR')THEN
  265. ISOUR=ICHP1
  266. ELSE IF(NOMTO3(1).EQ.'FX')THEN
  267. IFORC=ICHP1
  268. ENDIF
  269. IF(IRET2.EQ.1)THEN
  270. NBCOMP = 0
  271. INDIC = 1
  272. NOMTO3(1) = ' '
  273. NOMTO3(2) = ' '
  274. NOMTO3(3) = ' '
  275. CALL QUEPOI(ICHP2,ICENTR,INDIC,NBCOMP,NOMTO3)
  276. IF (IERR.NE.0) RETURN
  277. IF(NOMTO3(1).EQ.'SOUR')THEN
  278. ISOUR=ICHP2
  279. ELSE IF(NOMTO3(1).EQ.'FX')THEN
  280. IFORC=ICHP2
  281. ENDIF
  282. ENDIF
  283. ENDIF
  284. *
  285. * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  286. *
  287. MELEME = IELTFA
  288. SEGACT MELEME
  289. LZONES = LISOUS(/1)
  290. IF (LZONES.EQ.0) LZONES = 1
  291. IPT4 = IPGEOM
  292. SEGACT IPT4
  293. DO 30 ISOUS=1,NSOUS
  294. IMAMEL = MAHYBR(ISOUS)
  295. IF (IMAMEL.NE.0) THEN
  296. DO 20 ISZ=1,LZONES
  297. IF (LZONES.EQ.1) THEN
  298. IPT2 = IPT4
  299. IPT3 = MELEME
  300. ELSE
  301. IPT2 = IPT4.LISOUS(ISZ)
  302. IPT3 = LISOUS(ISZ)
  303. ENDIF
  304. IF (IPT2.EQ.IMAMEL) THEN
  305. MAHYBR(ISOUS) = IPT3
  306. GOTO 30
  307. ENDIF
  308. 20 CONTINUE
  309. IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  310. MOTERR(1:8) = ' MODELE '
  311. MOTERR(9:16) = ' ELTFA '
  312. INTERR(1) = ISOUS
  313. CALL ERREUR(698)
  314. SEGDES IPT4
  315. SEGDES MELEME
  316. GOTO 100
  317. ENDIF
  318. ENDIF
  319. 30 CONTINUE
  320. SEGDES IPT4
  321. SEGDES MELEME
  322. *
  323. * Test et tri du sous-type des matrices de rigiditée récupérées
  324. *
  325. IPFORC=0
  326. IPDARC=0
  327. SEGACT RI1
  328. LETYPE = RI1.MTYMAT
  329. IF(LETYPE.EQ.'DARCY')THEN
  330. IPDARC=RI1
  331. ELSE IF(LETYPE.EQ.'MASSE')THEN
  332. IPFORC=RI1
  333. ENDIF
  334. IF(RI2.NE.0)THEN
  335. SEGACT RI2
  336. LETYPE = RI2.MTYMAT
  337. IF(LETYPE.EQ.'MASSE')THEN
  338. IPFORC=RI2
  339. ELSE IF(LETYPE.EQ.'DARCY')THEN
  340. IPDARC=RI2
  341. ENDIF
  342. IF(IPFORC.EQ.0)THEN
  343. MOTERR(1:8) = 'RIGIDITE'
  344. MOTERR(9:16) = 'MASSE'
  345. CALL ERREUR(79)
  346. SEGDES RI1
  347. SEGDES RI2
  348. GOTO 100
  349. ELSE IF(IFORC.EQ.0)THEN
  350. MOTERR(1:8) = 'CHPOINT'
  351. MOTERR(9:16) = 'FORCE'
  352. CALL ERREUR(79)
  353. SEGDES RI1
  354. SEGDES RI2
  355. GOTO 100
  356. ENDIF
  357. ENDIF
  358. IF(IPDARC.EQ.0)THEN
  359. MOTERR(1:8) = 'RIGIDITE'
  360. MOTERR(9:16) = 'DARCY'
  361. CALL ERREUR(79)
  362. SEGDES RI1
  363. IF(IPFORC.NE.0)SEGDES RI2
  364. GOTO 100
  365. ENDIF
  366. *
  367. * Controle des pointeurs de MELEME de la rigidité
  368. *
  369. DO 40 ISOUS=1,NSOUS
  370. IMAMEL = MAHYBR(ISOUS)
  371. IF (IMAMEL.NE.0) THEN
  372. RI1=IPDARC
  373. IPTEST = RI1.IRIGEL(1,ISOUS)
  374. CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  375. IF (.NOT. LOGRE) THEN
  376. MOTERR(1:8) = 'DARCY'
  377. MOTERR(9:16) = ' ELTFA '
  378. INTERR(1) = ISOUS
  379. CALL ERREUR(698)
  380. SEGDES RI1
  381. GOTO 100
  382. ENDIF
  383. IF(IPFORC.NE.0)THEN
  384. RI2=IPFORC
  385. IPTEST = RI2.IRIGEL(1,ISOUS)
  386. CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  387. IF (.NOT. LOGRE) THEN
  388. MOTERR(1:8) = 'MASSE'
  389. MOTERR(9:16) = ' ELTFA '
  390. INTERR(1) = ISOUS
  391. CALL ERREUR(698)
  392. SEGDES RI1
  393. SEGDES RI2
  394. GOTO 100
  395. ENDIF
  396. ENDIF
  397. ENDIF
  398. 40 CONTINUE
  399. SEGDES RI1
  400. IF(IPFORC.NE.0)SEGDES RI2
  401. *
  402. * Vérification du support du MCHAML Ck|K| % au MMODEL
  403. * Controle du nombre de composantes reelles
  404. *
  405. IF (IPCK.NE.0) THEN
  406. IF(IPFORC.NE.0)CALL ERREUR(34)
  407. ISUP = 2
  408. ICOND = 1
  409. CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  410. IF (IRET.GT.1) GOTO 100
  411. MCHELM = IPCK
  412. SEGACT MCHELM
  413. DO 50 ISOUS=1,NSOUS
  414. IF (MAHYBR(ISOUS).NE.0) THEN
  415. MCHAML = ICHAML(ISOUS)
  416. SEGACT MCHAML
  417. N2 = TYPCHE(/2)
  418. IF (N2.GT.1) THEN
  419. CALL ERREUR(320)
  420. SEGDES MCHELM
  421. SEGDES MCHAML
  422. GOTO 100
  423. ENDIF
  424. CHAR6 = TYPCHE(1)(1:6)
  425. IF (CHAR6.NE.'REAL*8') THEN
  426. MOTERR(1:8) = NOMCHE(1)
  427. CALL ERREUR(679)
  428. SEGDES MCHELM
  429. SEGDES MCHAML
  430. GOTO 100
  431. ENDIF
  432. SEGDES MCHAML
  433. ENDIF
  434. 50 CONTINUE
  435. SEGDES MCHELM
  436. ENDIF
  437. * (n)
  438. * Test du CHAMPOINT face TH
  439. * On ne teste que la présence d'une composante de nom TH
  440. *
  441. MCHPO1 = ITPN
  442. SEGACT MCHPO1
  443. NSOUP1 = MCHPO1.IPCHP(/1)
  444. DO 70 ISOUS=1,NSOUP1
  445. MSOUPO = MCHPO1.IPCHP(ISOUS)
  446. SEGACT MSOUPO
  447. NC = NOCOMP(/2)
  448. DO 60 J=1,NC
  449. IF (NOCOMP(J).EQ.'TH ') GOTO 80
  450. 60 CONTINUE
  451. SEGDES MSOUPO
  452. 70 CONTINUE
  453. SEGDES MCHPO1
  454. MOTERR(1:4) = ' TH '
  455. CALL ERREUR(181)
  456. GOTO 100
  457. 80 CONTINUE
  458. SEGDES MSOUPO
  459. SEGDES MCHPO1
  460. *
  461. * Construction du CHPO des charges
  462. *
  463. SEGDES IPMAHY
  464. CALL HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE,
  465. S ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,ICHP1)
  466. IF (ICHP1.NE.0) CALL ECROBJ('CHPOINT',ICHP1)
  467. *
  468. * Ménage
  469. *
  470. 100 CONTINUE
  471. SEGSUP IPMAHY
  472. RETURN
  473. END
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  

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