Télécharger dmmu.eso

Retour à la liste

Numérotation des lignes :

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

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