Télécharger dmmu.eso

Retour à la liste

Numérotation des lignes :

dmmu
  1. C DMMU SOURCE CB215821 24/04/12 21:15:37 11897
  2. SUBROUTINE DMMU
  3. C-----------------------------------------------------------------------
  4. C -1 t
  5. C Calcul du produit D (M -U) B CHPface
  6. C -1 t
  7. C ou du produit D M B CHPface
  8. C le résultat est un CHPcentre
  9. C-----------------------------------------------------------------------
  10. C
  11. C---------------------------
  12. C Phrase d'appel (GIBIANE) :
  13. C---------------------------
  14. C
  15. C CHP3 = 'DMMU' MMODEL RIG1 (CHP1) CHP2 ;
  16. C
  17. C------------------------
  18. C Operandes et resultat :
  19. C------------------------
  20. C
  21. C MMODEL : modèle formulation DARCY.
  22. C RIG1 : Matrices hybrides elementaires de DARCY crees par MHYB.
  23. C CHP1 : flux aux faces
  24. C CHP2 : CHPO face à plusieur composantes.
  25. C CHP3 : CHPO centre à plusieur composantes dont le support géometrique
  26. C est le maillage centre du domaine les noms des composantes
  27. C sont ceux de CHP2
  28. C
  29. C-----------------------------------------------------------------------
  30. C
  31. C Langage : ESOPE + FORTRAN77
  32. C
  33. C
  34. C-----------------------------------------------------------------------
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8 (A-H,O-Z)
  37. *
  38.  
  39. -INC PPARAM
  40. -INC CCOPTIO
  41. -INC SMCOORD
  42. -INC SMCHPOI
  43. -INC SMCHAML
  44. -INC SMELEME
  45. -INC SMRIGID
  46. -INC SMTABLE
  47. -INC SMMODEL
  48. *
  49. SEGMENT ICCPR
  50. INTEGER ICPR(NNGOT)
  51. ENDSEGMENT
  52. SEGMENT IORGA
  53. INTEGER IAA(ITES)
  54. ENDSEGMENT
  55. *
  56. LOGICAL LOGRE,LOGIN
  57. CHARACTER*8 TAPIND,TYPOBJ,CHARRE,LETYPE
  58. CHARACTER*4 NOMTOT(1)
  59. *
  60. * Initialisations
  61. *
  62. IVALIN = 0
  63. XVALIN = 0.D0
  64. LOGIN = .TRUE.
  65. IOBIN = 0
  66. TAPIND = 'MOT '
  67. *
  68. *
  69. * Lecture du MMODEL
  70. *
  71. CALL LIROBJ('MMODEL',IPMODE,1,IRET)
  72. IF(IERR.NE.0)RETURN
  73. MMODEL=IPMODE
  74. SEGACT MMODEL
  75. N1=KMODEL(/1)
  76. DO 7 I=1,N1
  77. IMODEL=KMODEL(I)
  78. SEGACT IMODEL
  79. IF(FORMOD(1).NE.'DARCY')THEN
  80. MOTERR(1:16) = 'DARCY '
  81. CALL ERREUR(719)
  82. RETURN
  83. ENDIF
  84. 7 CONTINUE
  85. C on récupère la table DOMAINE
  86. IPTABL = 0
  87. CALL LEKMOD(MMODEL,IPTABL,IRET)
  88. IF (IERR.NE.0) RETURN
  89. TYPOBJ = 'MAILLAGE'
  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. IFACE = 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. *
  107. *
  108. * Test du sous-type de la matrice de rigiditée récupérée
  109. *
  110. SEGACT MRIGID
  111. LETYPE = MTYMAT
  112. IF (LETYPE.NE.'DARCY') THEN
  113. MOTERR(1:8) = 'RIGIDITE'
  114. MOTERR(9:16) = 'DARCY '
  115. CALL ERREUR(79)
  116. SEGDES MRIGID
  117. GOTO 100
  118. ENDIF
  119. *
  120. * Controle des pointeurs de MELEME de la rigidité
  121. *
  122. NRIGEL=IRIGEL(/2)
  123. C write(6,*)' NRIGEL----- ' ,nrigel
  124. MELEME=IELTFA
  125. SEGACT MELEME
  126. NBSOUS=LISOUS(/1)
  127. IF(NBSOUS.EQ.0)THEN
  128. IF((NRIGEL.NE.1).OR.(IRIGEL(1,1).NE.MELEME))THEN
  129. MOTERR(1:8) = 'DARCY '
  130. MOTERR(9:16) = 'ELTFA '
  131. INTERR(1) = 1
  132. CALL ERREUR(698)
  133. SEGDES MRIGID
  134. GOTO 100
  135. ENDIF
  136. ELSE
  137. IF(NRIGEL.NE.NBSOUS)THEN
  138. MOTERR(1:8) = 'DARCY '
  139. MOTERR(9:16) = 'ELTFA '
  140. INTERR(1) = 1
  141. CALL ERREUR(698)
  142. SEGDES MRIGID
  143. GOTO 100
  144. ENDIF
  145. DO 10 ISOUS=1,NBSOUS
  146. IF (LISOUS(ISOUS).NE.IRIGEL(1,ISOUS)) THEN
  147. MOTERR(1:8) = 'DARCY '
  148. MOTERR(9:16) = 'ELTFA '
  149. INTERR(1) = ISOUS
  150. CALL ERREUR(698)
  151. SEGDES MRIGID
  152. GOTO 100
  153. ENDIF
  154. 10 CONTINUE
  155. ENDIF
  156. *
  157. * LECTURE DES FLUX OU
  158. * DU CHPOIN FACE A MULTIPLIER
  159. *
  160. CALL LIROBJ('CHPOINT',IPCHF,1,IRET)
  161. IF(IRET.NE.1) RETURN
  162. MCHPO1=IPCHF
  163. NOMTOT(1)=' '
  164. NBCOMP=-1
  165. CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT)
  166. IF(IERR.NE.0)RETURN
  167. SEGACT MCHPO1
  168. MSOUP1=MCHPO1.IPCHP(1)
  169. SEGACT MSOUP1
  170. IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN
  171. MPOVA1=MSOUP1.IPOVAL
  172. SEGACT MPOVA1
  173. *
  174. *
  175. CALL LIROBJ('CHPOINT',IPCHF,1,IRET)
  176. IF(IRET.NE.1) RETURN
  177. MCHPO2=IPCHF
  178. NOMTOT(1)=' '
  179. NBCOMP=-1
  180. CALL QUEPOI(MCHPO2,IFACE,INDIC,NBCOMP,NOMTOT)
  181. IF(IERR.NE.0)RETURN
  182. SEGACT MCHPO2
  183. MSOUP2=MCHPO2.IPCHP(1)
  184. SEGACT MSOUP2
  185. MPOVA2=MSOUP2.IPOVAL
  186. SEGACT MPOVA2
  187. ELSE
  188. MCHPO2=MCHPO1
  189. MCHPO1=0
  190. MSOUP2=MSOUP1
  191. MPOVA2=MSOUP2.IPOVAL
  192. SEGACT MPOVA2
  193. CALL LIROBJ('CHPOINT',IPCHF,0,IRET)
  194. IF(IRET.EQ.1) THEN
  195. MCHPO1=IPCHF
  196. NOMTOT(1)=' '
  197. NBCOMP=-1
  198. CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT)
  199. IF(IERR.NE.0)RETURN
  200. SEGACT MCHPO1
  201. MSOUP1=MCHPO1.IPCHP(1)
  202. SEGACT MSOUP1
  203. IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN
  204. MPOVA1=MSOUP1.IPOVAL
  205. SEGACT MPOVA1
  206. ELSE
  207. CALL ERREUR(21)
  208. RETURN
  209. ENDIF
  210. ENDIF
  211. ENDIF
  212. *
  213. * Construction de MCHPOI
  214. *
  215. *
  216. IPT1=ICENTR
  217. SEGACT IPT1
  218. NPN=IPT1.NUM(/2)
  219. NSOUPO=1
  220. NAT=1
  221. SEGINI MCHPOI
  222. MTYPOI=' '
  223. MOCHDE=' CHPOIN CREE PAR DMMU '
  224. IFOPOI=IFOUR
  225. JATTRI(1)=2
  226. NC=NBCOMP
  227. SEGINI MSOUPO
  228. IPCHP(1)=MSOUPO
  229. DO 5 L=1,NBCOMP
  230. NOCOMP(L)=MSOUP2.NOCOMP(L)
  231. NOHARM(L)=MSOUP2.NOHARM(L)
  232. 5 CONTINUE
  233. IGEOC=ICENTR
  234. N=NPN
  235. SEGINI MPOVAL
  236. IPOVAL=MPOVAL
  237. NB=N*NC
  238. CALL INITD(VPOCHA,NB,0.D0)
  239. *
  240. * Creation du tableau ICPR
  241. *
  242. IK = 0
  243. NNGOT = nbpts
  244. SEGINI ICCPR
  245. MELEME = IFACE
  246. SEGACT MELEME
  247. N2 = NUM(/2)
  248. DO 15 I2=1,N2
  249. K = NUM(1,I2)
  250. IF (ICPR(K).EQ.0) THEN
  251. IK = IK + 1
  252. ICPR(K) = IK
  253. ENDIF
  254. 15 CONTINUE
  255. ITES=IK
  256. SEGDES MELEME
  257. SEGINI IORGA
  258. MELEME = MSOUP2.IGEOC
  259. SEGACT MELEME
  260. N2 = NUM(/2)
  261. DO 25 I2=1,N2
  262. K = NUM(1,I2)
  263. IF (ICPR(K).EQ.0) THEN
  264. INTERR(1) = K
  265. MOTERR(1:8) = 'FACE '
  266. CALL ERREUR(64)
  267. SEGSUP ICCPR, IORGA
  268. RETURN
  269. ELSE
  270. IAA(ICPR(K)) = I2
  271. ENDIF
  272. 25 CONTINUE
  273. SEGDES MELEME
  274. CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  275. MCHELM=IORIE
  276. SEGACT MCHELM
  277. C
  278. C Calcul du produit
  279. C
  280. MELEME=IELTFA
  281. SEGACT MELEME
  282. IPT1=MELEME
  283. ITELEM=0
  284. IF(MCHPO1.NE.0)THEN
  285. C
  286. C Cas avec vitesses
  287. C
  288. DO 50 ISOUS=1,NRIGEL
  289. MCHAML= ICHAML(ISOUS)
  290. SEGACT MCHAML
  291. MELVAL=IELVAL(1)
  292. SEGACT MELVAL
  293. IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS)
  294. SEGACT IPT1
  295. xMATRI=IRIGEL(4,ISOUS)
  296. SEGACT xMATRI
  297. NELRIG=re(/3)
  298. C write(6,*)' ISOUS NELRIG ',ISOUS,NELRIG
  299. DO 60 IEL=1,NELRIG
  300. ITELEM=ITELEM+1
  301. * XMATRI=IMATTT(IEL)
  302. * SEGACT XMATRI
  303. NLIGRD=RE(/1)
  304. NLIGRP=RE(/2)
  305. C write(6,*) ' NLIGRD NLIGRP ',NLIGRD,NLIGRP
  306. DO 40 I=1,NLIGRD
  307. RLIGN=0.D0
  308. DO 30 J=1,NLIGRP
  309. RLIGN=RLIGN+RE(I,J,iel)
  310. 30 CONTINUE
  311. IPOPTS=ICPR(IPT1.NUM(I,IEL))
  312. III=IAA(IPOPTS)
  313. C write(6,*)' dmmu ',itelem,i,rlign,MPOVA1.VPOCHA(IPOPTS,1)
  314. C * , VELCHE(I,IEL)
  315. RLIGN=RLIGN-MPOVA1.VPOCHA(IPOPTS,1)*VELCHE(I,IEL)
  316. DO 20 K=1,NBCOMP
  317. VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN*
  318. * MPOVA2.VPOCHA(III,K)
  319. 20 CONTINUE
  320.  
  321. 40 CONTINUE
  322. * SEGDES XMATRI
  323. 60 CONTINUE
  324. SEGDES xMATRI,MCHAML,MELVAL
  325. 50 CONTINUE
  326. ELSE
  327. C
  328. C Cas ou il n'y a pas de vitesse
  329. C On a dupliqué toute la boucle 50 pour ne pas faire
  330. C NRIGEL*NELRIG*NLIGRD fois le meme test
  331. C
  332. DO 51 ISOUS=1,NRIGEL
  333. MCHAML= ICHAML(ISOUS)
  334. SEGACT MCHAML
  335. MELVAL=IELVAL(1)
  336. SEGACT MELVAL
  337. IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS)
  338. SEGACT IPT1
  339. xMATRI=IRIGEL(4,ISOUS)
  340. SEGACT xMATRI
  341. NELRIG=re(/3)
  342. DO 61 IEL=1,NELRIG
  343. ITELEM=ITELEM+1
  344. * XMATRI=IMATTT(IEL)
  345. * SEGACT XMATRI
  346. NLIGRD=RE(/1)
  347. NLIGRP=RE(/2)
  348. DO 41 I=1,NLIGRD
  349. RLIGN=0.D0
  350. DO 31 J=1,NLIGRP
  351. RLIGN=RLIGN+RE(I,J,iel)
  352. 31 CONTINUE
  353. IPOPTS=ICPR(IPT1.NUM(I,IEL))
  354. III=IAA(IPOPTS)
  355. DO 21 K=1,NBCOMP
  356. VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN*
  357. * MPOVA2.VPOCHA(III,K)
  358. 21 CONTINUE
  359.  
  360. 41 CONTINUE
  361. * SEGDES XMATRI
  362. 61 CONTINUE
  363. SEGDES xMATRI,MCHAML,MELVAL
  364. 51 CONTINUE
  365. ENDIF
  366. SEGDES MRIGID
  367. CALL ECROBJ('CHPOINT',MCHPOI)
  368. *
  369. * Ménage
  370. *
  371. SEGDES MELEME
  372. SEGSUP ICCPR ,IORGA
  373. SEGDES MPOVAL,MSOUPO,MCHPOI
  374. SEGDES MPOVA2,MSOUP2,MCHPO2
  375. IF(MCHPO1.NE.0) THEN
  376. SEGDES MPOVA1,MSOUP1,MCHPO1
  377. ENDIF
  378. 100 CONTINUE
  379. RETURN
  380. END
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  

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