Télécharger kress.eso

Retour à la liste

Numérotation des lignes :

  1. C KRESS SOURCE PV 16/11/17 22:00:27 9180
  2. SUBROUTINE KRESS(MTABP,MCHB,MCHR,IMPR,BETA,KDPDQ,KPIMP,PIMP,NIMP)
  3. C ***********************************************************************
  4. C
  5. C RESOLUTION DE LA PRESSION DANS LE CAS D'UNE METHODE ITERATIVE
  6. C AVEC MATRICE STOCKEE. ON ENCHAINE, SI BESOIN EST : PROGCS,
  7. C CGPR ET RESCGS.
  8. C
  9. C DANS LE CAS D'UNE METHODE ITERATIVE AVEC STOCKAGE DE LA MATRICE
  10. C DE PRESSION ....
  11. C
  12. C CE SOUS PROGRAMME EFFECTUE LES OPERATIONS SUIVANTES :
  13. C
  14. C 1. STOCKAGE DANS LA TABLE METHODE DES POINTEURS ET DONNEES NECESSAIRES
  15. C
  16. C 2. CALCUL ET STOCKAGE DE LA MATRICE C D C CONNAISSANT C
  17. C EN FONCTION DU MODE DE STOCKAGE DEMANDE
  18. C
  19. C DEUX MODES DE STOCKAGE SONT DISPONIBLES
  20. C - MORSE CLASSIQUE (KSTOCK.EQ.0)
  21. C - MODE COMPRESSE (KSTOCK.EQ.1)
  22. C
  23. C 3. CALCUL D'UN PRECONDITIONNEMENT SI NECESSAIRE, I.E. LORSQUE
  24. C KTYPI=3 OU 4
  25. C
  26. C
  27. C
  28. C
  29. C POINTEURS :
  30. C
  31. C MTABP CONTIENT ENTRE AUTRES LES TABLEAUX D'INDIRECTION
  32. C (CF PROFIM)
  33. C MELEME OBJET MAILLAGE SUR LEQUEL EST FAIT LE CALCUL
  34. C MATRAK MATRICES ELEMENTAIRES DE LA DIVERGENCE (ALIAS"C")
  35. C
  36. C IZIPAD CORRESPONDANCE NUMER. GLOBALE --> NUMER. LOCALE
  37. C (DOMAINE SUR LEQUEL PORTE AP ET NON LA PRESSION)
  38. C IZTAB POROSITE (VECT DOMA) OU (VECT NOEU) PAR DEFAUT 1.,1.,1.
  39. C
  40. C EN SORTIE :
  41. C
  42. C LA MATRICE EST RANGEE DANS UN SEUL SEGMENT DE MEMOIRE
  43. C
  44. C MORSE : PAR LIGNES, LE PREMIER DE LA LIGNE EST L'ELEMENT LUI-MEME
  45. C COMPRESSE : COMME UN TABLEAU A(NEL,NNZ) AVEC NNZ = NOMBRE MAXI
  46. C DE TERMES NON NULS PAR LIGNE. LA DIAGONALE EST STOCKEE
  47. C EN PREMIER.
  48. C
  49. C LE CALCUL DE LA MATRICE EST FAIT SOUS-OBJET PAR SOUS-OBJET CE QUI
  50. C LIMITE A DEUX LE NOMBRE DE MELEMEs ET DE MATRICES DE LA DIVERGENCE
  51. C ACTIVES AU MEME INSTANT (EN PLUS DE LA MATRICE ET DE SON TABLEAU
  52. C D'INDEXAGE)
  53. C
  54. C ***********************************************************************
  55. IMPLICIT INTEGER(I-N)
  56. IMPLICIT REAL*8 (A-H,O-Z)
  57. C
  58. C***
  59.  
  60.  
  61. -INC PPARAM
  62. -INC CCOPTIO
  63. -INC SMELEME
  64. POINTEUR IPTJ.MELEME,IPTK.MELEME,MELEM1.MELEME
  65. POINTEUR IGEOM.MELEME,IGEOMI.MELEME
  66. POINTEUR MELEMK.MELEME,MELEMC.MELEME,IGEOMC.MELEME
  67. POINTEUR MELSTB.MELEME
  68. -INC SMLENTI
  69. POINTEUR IZIPAD.MLENTI,MLENTP.MLENTI,IZIPAK.MLENTI
  70. POINTEUR IA.MLENTI,JA.MLENTI,KA.MLENTI
  71. -INC SMCOORD
  72. -INC SMCHPOI
  73. POINTEUR IZDV.MCHPOI
  74. POINTEUR IZDP.MCHPOI,IZDDP.MPOVAL
  75. POINTEUR MCHB.MCHPOI,MCHR.MCHPOI
  76. POINTEUR IZR.MPOVAL,IZPOC.MPOVAL
  77. POINTEUR IZP.MPOVAL,IZF.MPOVAL
  78. POINTEUR IPRESU.MCHPOI
  79. POINTEUR IPRESV.MPOVAL
  80. C-INC SMMATRAK
  81. C*************************************************************************
  82. C
  83. C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees
  84. C
  85.  
  86. * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange
  87. * (points CENTRE ) pour chaque operateur de contrainte
  88. * KGEOC SPG pour la totalite des points CENTRE.
  89. * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse)
  90. * KLEMC Connectivites de l'ensemble des contraintes
  91. * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones
  92.  
  93. SEGMENT MATRAK
  94. INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP)
  95. INTEGER LIZAFM(NBSOUS)
  96. INTEGER IKAM0 (NBSOUS)
  97. INTEGER IMEM (NBELC)
  98. INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC
  99. ENDSEGMENT
  100.  
  101. SEGMENT IZAFM
  102. REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX)
  103. ENDSEGMENT
  104.  
  105. POINTEUR IPMJ.IZAFM,IPMK.IZAFM
  106.  
  107. C*******************************************************************
  108. -INC SMLREEL
  109. POINTEUR IZB.MLREEL
  110. POINTEUR IZDPR.MLREEL,IZPIM.MLREEL
  111. POINTEUR IZICCG.MLREEL,IZA.MLREEL
  112. POINTEUR IZD.MLREEL
  113. C
  114. CHARACTER*8 TYPE,TYPC
  115. C
  116. C Pour verifier ou non les numeros des elements sur
  117. C une condition de sortie fluide.
  118. C
  119. LOGICAL SOMTER
  120. DATA SOMTER/.FALSE./
  121.  
  122. C***
  123. IZL=0
  124. KDPDQ=0
  125. C write(6,*)' Entree dans kress KDPDQ=',kdpdq
  126.  
  127. C
  128. C*** LECTURE OU CALCUL DES TABLEAUX D'INDIRECTION ***
  129. C
  130. C*** On essaie de lire les tableaux d'indirection crees par PROGCS
  131. C Si la table MATRIS n'existe pas, alors LEKTAB fera appel
  132. C a PROGCS.
  133.  
  134. CALL LEKTAB(MTABP,'MATRIS',MTAB2)
  135. IF(MTAB2.EQ.0)GO TO 90
  136. TYPE=' '
  137. CALL ACMO(MTABP,'METHODE',TYPE,MTAB3)
  138. IF(MTAB3.EQ.0)GO TO 90
  139. CALL ACME(MTAB3,'KTYPI',KTYPI)
  140. CALL ACME(MTAB3,'KSTOCK',KSTO)
  141.  
  142. C
  143. C*** CALCUL DES COEFFICIENTS DE LA MATRICE ***
  144. C
  145.  
  146. TYPE=' '
  147. CALL ACMO(MTABP,'MATC',TYPE,MATRAK)
  148. IF(TYPE.NE.'MATRAK')THEN
  149. WRITE(6,*)' Il n''y a pas d''entree MATC dans la table '
  150. GO TO 90
  151. ENDIF
  152.  
  153. TYPE=' '
  154. CALL ACMO(MTABP,'DIAGV',TYPE,IZDV)
  155. IF(TYPE.NE.'CHPOINT ')THEN
  156. WRITE(6,*)' Il n''y a pas d''entree DIAGV dans la table '
  157. GO TO 90
  158. ENDIF
  159.  
  160. SEGACT MATRAK*MOD
  161.  
  162. IF(KDIAG.EQ.IZDV.AND.KCAC.EQ.1)THEN
  163. GO TO 500
  164. ENDIF
  165.  
  166. IF(KDPDQ.NE.0)THEN
  167. C DPDQ
  168. TYPE=' '
  169. CALL ACMO(MTABP,'DOMAINE',TYPE,MTABD)
  170. IF(TYPE.NE.'TABLE')GO TO 90
  171. TYPE=' '
  172. CALL ACMO(MTABD,'MELSTB',TYPE,MELSTB)
  173. IF(TYPE.NE.'MAILLAGE') GO TO 90
  174. CALL LEKTAB(MTABD,'CENTRE',MELEMC)
  175.  
  176. TYPE=' '
  177. CALL ACMO(MTABD,'MCHPOC',TYPE,MCHPOC)
  178. IF(TYPE.NE.'CHPOINT') GO TO 90
  179. CALL LICHT(MCHPOC,IZPOC,TYPC,IGEOMC)
  180. NBK=IZPOC.VPOCHA(/2)
  181. SEGACT MELEMC,IGEOMC
  182. C DPDQ
  183. ENDIF
  184.  
  185. KDIAG=IZDV
  186. CALL LICHT(IZDV,MPOVAL,TYPC,IGEOM)
  187.  
  188.  
  189.  
  190. WRITE(6,*)' *------------------------------------------------*'
  191. WRITE(6,*)' * KTYPI=2,3,4 RESOLUTION DE L''EQUATION DE *'
  192. WRITE(6,*)' * PRESSION PAR UNE METHODE ITERATIVE. *'
  193. WRITE(6,*)' * LA MATRICE EST CALCULEE ET STOCKEE UNE FOIS *'
  194. WRITE(6,*)' * POUR TOUTES. L''INVERSION UTILISE UNE *'
  195. WRITE(6,*)' * METHODE DE GRADIENT CONJUGUE PRECONDITIONNE *'
  196. WRITE(6,*)' * PAR LA DIAGONALE OU UNE FACTORISATION *'
  197. WRITE(6,*)' * INCOMPLETE. *'
  198. C
  199. C*********************************************************************
  200. C--- Initialisation des segments IZDPR IZPIM
  201. C et stockage de leurs pointeurs dans la table METHODE.
  202. C*********************************************************************
  203. C
  204. C--- NL nombre total d'elements
  205. C
  206. CALL ACME(MTAB3,'NL',NL)
  207. JG=NL
  208. SEGINI IZDPR
  209. CALL ECMO(MTAB3,'DPR','LISTREEL',IZDPR)
  210. SEGDES IZDPR
  211. IF(KTYPI.EQ.4) THEN
  212. IF(KSTO.EQ.0) CALL ACME(MTAB2,'LA',NICCG)
  213. IF(KSTO.EQ.1) CALL ACME(MTAB2,'NNZ',NNZ)
  214. IF(KSTO.EQ.1) NICCG=NL*NNZ
  215. JG=NICCG
  216. SEGINI IZICCG
  217. CALL ECMO(MTAB3,'ICCG','LISTREEL',IZICCG)
  218. SEGDES IZICCG
  219. ENDIF
  220. C
  221. SEGINI IZPIM
  222. CALL INITD (IZPIM.PROG,NL,0.D0)
  223. CALL ECMO(MTAB3,'PIM','LISTREEL',IZPIM)
  224. SEGDES IZPIM
  225. C
  226. C*********************************************************************
  227. C---- P A R A M E T R E S D E C O N V E R G E N C E
  228. C*********************************************************************
  229. C
  230. CALL ACME(MTAB3,'NITMAX',NITMAX)
  231. CALL ACMF(MTAB3,'EPSI',EPI)
  232. CALL ACME(MTAB3,'NPITE',NPITE)
  233. CALL ACME(MTAB3,'NFIMPR',NFIMPR)
  234. WRITE(6,*)' * *'
  235. WRITE(6,*)' * PARAMETRES DE CONTROLE CHOISIS : *'
  236. WRITE(6,*)' * *'
  237. WRITE(6,*)' * KTYPI : ',KTYPI
  238. WRITE(6,*)' * KSTOCK : ',KSTO
  239. WRITE(6,*)' * NITMAX : ',NITMAX
  240. WRITE(6,666) EPI
  241. WRITE(6,*)' * JTPA : ',NPITE
  242. WRITE(6,*)' * FIMPR : ',NFIMPR
  243. WRITE(6,*)' * *'
  244. 666 FORMAT(1X,' * JTEPS : ',3X,E8.2)
  245. C
  246. C*********************************************************************
  247. C---- I N I T I A L I S A T I O N D E L A P R E S S I O N
  248. C*********************************************************************
  249. C
  250.  
  251. C
  252. C AU CAS ON ON FOURNIT UNE PRESSION INITIALE, ON SAUVE LE CONTENU
  253. C DANS PIM.
  254. C
  255. TYPE=' '
  256. CALL ACMO(MTAB3,'PINITIAL',TYPE,IPRESU)
  257. C
  258. IF (IPRESU.EQ.0) THEN
  259. SEGACT IZPIM*MOD
  260. DO II=1,NL
  261. IZPIM.PROG(II)=0.D0
  262. ENDDO
  263. SEGDES IZPIM
  264. ELSE
  265. CALL LICHT(IPRESU,IPRESV,TYPC,IGEOM)
  266. SEGACT IPRESV*MOD
  267. LL1=IPRESV.VPOCHA(/2)
  268. LL2=IPRESV.VPOCHA(/1)
  269. IF (LL1.NE.1.OR.LL2.NE.NL) THEN
  270. WRITE(6,*)'PB DANS CGPR : DIM DE LA PRESSION'
  271. WRITE(6,*)'DIM ATTENDUE : ',NL
  272. WRITE(6,*)'DIM RECUE : ',LL2
  273. IF(NL.LT.LL2) THEN
  274. WRITE(6,*) 'NOMBRES DE CONTRAINTES INCOMPATIBLES'
  275. STOP
  276. ELSE
  277. WRITE(6,*) 'IL Y A DES CONTRAINTES SUPPLEMENTAIRES'
  278. SEGACT IZPIM*MOD
  279. DO II=1,LL2
  280. IZPIM.PROG(II)=IPRESV.VPOCHA(II,1)
  281. ENDDO
  282. DO II=LL2+1,NL
  283. IZPIM.PROG(II)=0.D0
  284. ENDDO
  285. SEGDES IZPIM,IPRESV
  286. ENDIF
  287. ENDIF
  288. SEGACT IZPIM*MOD
  289. DO 15 II=1,NL
  290. IZPIM.PROG(II)=IPRESV.VPOCHA(II,1)
  291. 15 CONTINUE
  292. SEGDES IZPIM,IPRESU,IPRESV
  293. ENDIF
  294. C
  295. C********************************************************************
  296. C--- A L L O C A T I O N D E L A M A T R I C E
  297. C********************************************************************
  298. C
  299.  
  300. C--- Allocation de memoire pour la matrice A et la diagonale inverse
  301.  
  302. IF(KSTO.EQ.0) CALL ACME(MTAB2,'LA',LA)
  303. IF(KSTO.EQ.1) CALL ACME(MTAB2,'NNZ',NNZ)
  304. IF(KSTO.EQ.1) LA=NL*NNZ
  305. JG=LA
  306. SEGINI IZA
  307. CALL ECMO(MTAB2,'IMAT','LISTREEL',IZA)
  308. SEGINI IZD
  309. C
  310. C********************************************************************
  311. C--- C A L C U L D E L A M A T R I C E D E P R E S S I O N
  312. C********************************************************************
  313. C
  314. MELEME=KLEMC
  315. C
  316. C ON COMMENCE PAR ACTIVER
  317. C
  318.  
  319. MELEM1=KGEOS
  320. CALL KRIPAD(MELEM1,IZIPAD)
  321.  
  322. SEGACT MELEME
  323.  
  324. C
  325. C Nombre de sous-objets
  326. C
  327. NBSOUS=LISOUS(/1)
  328. IF(NBSOUS.EQ.0)NBSOUS=1
  329.  
  330. IES=VPOCHA(/2)
  331. NPT=VPOCHA(/1)
  332. C
  333.  
  334. C--- S T O C K A G E M O R S E
  335. IF(KSTO.EQ.0) THEN
  336. TYPE=' '
  337. CALL ACMO(MTAB2,'IA',TYPE,IA)
  338. TYPE=' '
  339. CALL ACMO(MTAB2,'JA',TYPE,JA)
  340. SEGACT IA*MOD,JA*MOD
  341. KKA=0
  342. C
  343. C On procede sous-objet par sous-objet
  344. C
  345.  
  346. C? write(6,*)' KDPDQ=',KDPDQ
  347. IF(KDPDQ.NE.0)THEN
  348. C
  349. C OPERATEUR DPDQ
  350. C
  351. CALL KRIPAD(IGEOMC,MLENTI)
  352. SEGACT MELSTB*MOD,IGEOMC
  353. MELEMK=KGEOC
  354. CALL KRIPAD(MELEMK,IZIPAK)
  355. SEGACT MELEMK*MOD
  356. ENDIF
  357.  
  358. KJ=0
  359. DO 10 JIS=1,NBSOUS
  360. IF(NBSOUS.EQ.1) IPTJ=MELEME
  361. IF(NBSOUS.NE.1) IPTJ=LISOUS(JIS)
  362. SEGACT IPTJ*MOD
  363. IPMJ=LIZAFM(JIS)
  364. JKAM=IKAM0 (JIS)
  365. SEGACT IPMJ*MOD
  366. NPJ=IPTJ.NUM(/1)
  367. NEL=IPTJ.NUM(/2)
  368. C Un sous-objet est active, on traite tous ses elements
  369. DO 1 KJJ=1,NEL
  370. KJ=KJ+1
  371. C
  372. C--- Boucle sur le nombe de voisins, les numeros de ceux-ci sont
  373. C dans JA
  374. C
  375. C On v‚rifie que la somme des termes d'une ligne de la matrice est nulle
  376. C en accumulant dans la variable STERM.
  377. C
  378. STERM=0.D0
  379. DO 2 KV=IA.LECT(KJ),IA.LECT(KJ+1)-1
  380. KKA=KKA+1
  381. KK=JA.LECT(KV)
  382. C De deux choses l'une, KK appartient au meme sous-objet
  383. C que KJ ou non. Par defaut on dit oui puis on teste.
  384. IPTK=IPTJ
  385. IPMK=IPMJ
  386. KKAM=JKAM
  387. NPK=NPJ
  388. NSOBJ=IMEM(KK)
  389. IF(NSOBJ.NE.JIS) THEN
  390. C KK et KJ sont dans deux sous-objets differents, il faut
  391. C activer celui de KK
  392. IPTK=LISOUS(NSOBJ)
  393. IPMK=LIZAFM(NSOBJ)
  394. KKAM=IKAM0 (NSOBJ)
  395. SEGACT IPTK*MOD
  396. SEGACT IPMK*MOD
  397. NPK=IPTK.NUM(/1)
  398. ENDIF
  399. KKK=KK-KKAM+1
  400.  
  401. AUX=0.D0
  402.  
  403. DO 31 I=1,NPJ
  404. DO 32 N=1,NPK
  405. IF(IPTJ.NUM(I,KJJ).NE.IPTK.NUM(N,KKK))GO TO 3
  406. JD=IZIPAD.LECT(IPTK.NUM(N,KKK))
  407. * IPR=IC*(JD-1)+1
  408. AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1)
  409. & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2)
  410. IF(IES.EQ.2)GO TO 3
  411. AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3)
  412. 3 CONTINUE
  413. 32 CONTINUE
  414. 31 CONTINUE
  415. IZA.PROG(KKA)=AUX
  416. IF(KK.EQ.KJ) IZD.PROG(KJ)=AUX
  417.  
  418. C? write(6,*)' KDPDQ 2 = ',KDPDQ
  419. IF(KDPDQ.NE.0)THEN
  420. C
  421. C OPERATEUR DPDQ
  422. C
  423. C PV nuna est dans le segment izl. izl n'est pas initialisé
  424. if (izl.eq.0) call erreur(5)
  425. KJQ=NUNA(KJ)
  426. LPJ=LECT(MELEMK.NUM(1,KJQ))
  427. IF(LPJ.NE.0)THEN
  428. P1=BETA*IZPOC.VPOCHA(LPJ,1)
  429.  
  430. C?! D(KJ)=D(KJ)+ABS(P1)
  431. IZD.PROG(KJ)=IZD.PROG(KJ)+ABS(P1)
  432.  
  433. DO 3819 JC=2,NBK
  434. NPC=MELSTB.NUM(JC,LPJ)
  435. P=BETA*IZPOC.VPOCHA(LPJ,JC)
  436. KPI=IZIPAK.LECT(NPC)
  437. KPJN=NUAN(KPI)
  438. IF(KPJN.GT.KJ)GO TO 3819
  439. IF(KPJN.GT.KJ)GO TO 3819
  440.  
  441. C# MC : IDEBLK non initialisé + IDECI et LLI uniquement dans KJAA
  442. C KJAA non utilisé => on vire !!!
  443. C IDECI=IDEBLK(KJ-KJD+1)
  444. C LLI=IDEBLK(KJ-KJD+2)-IDEBLK(KJ-KJD+1)
  445. C KJAA=LLI-KJ+KPJN+IDECI
  446.  
  447. C?! A(KJAA)=A(KJAA)+P
  448. IZA.PROG(KKA)=IZA.PROG(KKA)+P
  449.  
  450. 3819 CONTINUE
  451. ENDIF
  452. ENDIF
  453. C FIN DPDQ
  454.  
  455. C
  456. STERM=STERM+AUX
  457. C
  458. C Il faut desactiver IPMK et IPTK sauf bien sur si on n'a
  459. C qu'un seul sous-objet
  460. C
  461. IF(NSOBJ.NE.JIS) THEN
  462. SEGDES IPMK
  463. SEGDES IPTK
  464. ENDIF
  465. 2 CONTINUE
  466. C
  467. C Controle sur la somme des termes d'une ligne
  468. C
  469. IF(SOMTER.AND.ABS(STERM).GT.1.D-10) THEN
  470. WRITE(6,*) 'L''ELEMENT KJ = ',KJ,' DOIT ETRE SUR UNE SORTIE'
  471. ENDIF
  472. 1 CONTINUE
  473. SEGDES IPTJ
  474. SEGDES IPMJ
  475. 10 CONTINUE
  476. C
  477. C Sortie
  478. C dimension ts(9)
  479. C do 1001 kj=1,nelp
  480. C write(61,60) (itab(ka),ka=iat(kj),iat(kj+1)-1)
  481. C1001 continue
  482. C do 1000 kj=1,nelp
  483. C call initd(ts,9,0.d0)
  484. C do 1003 ka=iat(kj),iat(kj+1)-1
  485. C i=ka-iat(kj)+1
  486. C ts(i)=a(ka)
  487. C1003 continue
  488. C write(61,61) (a(ka),ka=iat(kj),iat(kj+1)-1)
  489. C write(61,61) (ts(i),i=1,9)
  490. C1000 continue
  491. 60 format(9(1x,i4))
  492. 61 format(9(1x,f7.4))
  493. C
  494. SEGDES IA,JA,IZA
  495. SEGDES MELEME
  496.  
  497. IF(KDPDQ.NE.0)THEN
  498. SEGSUP IZIPAK,MLENTI
  499. ENDIF
  500.  
  501.  
  502. C--- S T O C K A G E C O M P R E S S E
  503.  
  504. ELSEIF(KSTO.EQ.1) THEN
  505. TYPE=' '
  506. CALL ACMO(MTAB2,'KA',TYPE,KA)
  507. SEGACT KA*MOD
  508. KJ=0
  509. DO 20 JIS=1,NBSOUS
  510. IF(NBSOUS.EQ.1) IPTJ=MELEME
  511. IF(NBSOUS.NE.1) IPTJ=LISOUS(JIS)
  512. SEGACT IPTJ*MOD
  513. IPMJ=LIZAFM(JIS)
  514. JKAM=IKAM0 (JIS)
  515. SEGACT IPMJ*MOD
  516. NPJ=IPTJ.NUM(/1)
  517. NEL=IPTJ.NUM(/2)
  518. C Un sous-objet est active, on traite tous ses elements
  519. DO 101 KJJ=1,NEL
  520. KJ=KJ+1
  521. C
  522. C--- Boucle sur le nombre de voisins, la liste se termine lorsque
  523. C le numero du voisin est egal au numero de l'element.
  524. C KA et A sont ranges comme si ils etaient declares KA(NL,NNZ)
  525. C et A(NL,NNZ).
  526. C
  527. STERM=0.D0
  528. DO 102 KV=1,NNZ
  529. KKA=KJ+(KV-1)*NL
  530. KK=KA.LECT(KKA)
  531. IF(KV.GT.1.AND.KK.EQ.KJ) THEN
  532. C--- On est arrive a la fin de la liste des voisins
  533. IZA.PROG(KKA)=0.D0
  534. GOTO 102
  535. ENDIF
  536. C De deux choses l'une, KK appartient au meme sous-objet
  537. C que KJ ou non. Par defaut on dit oui puis on teste.
  538. IPTK=IPTJ
  539. IPMK=IPMJ
  540. KKAM=JKAM
  541. NPK=NPJ
  542. NSOBJ=IMEM(KK)
  543. IF(NSOBJ.NE.JIS) THEN
  544. C KK et KJ sont dans deux sous-objets differents, il faut
  545. C activer celui de KK
  546. IPTK=LISOUS(NSOBJ)
  547. IPMK=LIZAFM(NSOBJ)
  548. KKAM=IKAM0 (NSOBJ)
  549. SEGACT IPTK*MOD
  550. SEGACT IPMK*MOD
  551. NPK=IPTK.NUM(/1)
  552. ENDIF
  553. KKK=KK-KKAM+1
  554.  
  555. AUX=0.D0
  556. DO 131 I=1,NPJ
  557. DO 132 N=1,NPK
  558. IF(IPTJ.NUM(I,KJJ).NE.IPTK.NUM(N,KKK))GO TO 103
  559. JD=IZIPAD.LECT(IPTK.NUM(N,KKK))
  560. * IPR=IC*(JD-1)+1
  561. AUX=AUX+IPMJ.AM(KJJ,I,1)*IPMK.AM(KKK,N,1)/VPOCHA(JD,1)
  562. & +IPMJ.AM(KJJ,I,2)*IPMK.AM(KKK,N,2)/VPOCHA(JD,2)
  563. IF(IES.EQ.2)GO TO 103
  564. AUX=AUX+IPMJ.AM(KJJ,I,3)*IPMK.AM(KKK,N,3)/VPOCHA(JD,3)
  565. 103 CONTINUE
  566. 132 CONTINUE
  567. 131 CONTINUE
  568. IZA.PROG(KKA)=AUX
  569. IF(KV.EQ.1) IZD.PROG(KJ)=AUX
  570. C
  571. C Il faut desactiver IPMK et IPTK sauf bien sur si on n'a
  572. C qu'un seul sous-objet
  573. C
  574. STERM=STERM+AUX
  575. IF(NSOBJ.NE.JIS) THEN
  576. SEGDES IPMK
  577. SEGDES IPTK
  578. ENDIF
  579. 102 CONTINUE
  580. IF(SOMTER.AND.ABS(STERM).GT.1.D-10) THEN
  581. WRITE(6,*) 'L''ELEMENT KJ = ',KJ,' DOIT ETRE SUR UNE SORTIE'
  582. ENDIF
  583. 101 CONTINUE
  584. SEGDES IPTJ
  585. SEGDES IPMJ
  586. 20 CONTINUE
  587. C
  588. C Sortie
  589. C do 2001 kj=1,nelp
  590. C write(60,60) (itab(kj+(kv-1)*nelp),kv=1,nnz)
  591. C2001 continue
  592. C do 2000 kj=1,nelp
  593. C write(60,61) ( a (kj+(kv-1)*nelp),kv=1,nnz)
  594. C2000 continue
  595. SEGDES KA,IZA
  596. SEGDES MELEME
  597. ELSE
  598. WRITE(6,*) 'CGPR Mode de stockage inconnu |'
  599. ENDIF
  600. C
  601. C*** RAJOUTER LA DIAGONALE DE LA PRESSION ***
  602. C
  603. TYPE=' '
  604. CALL ACMO(MTABP,'DIAGP',TYPE,IZDP)
  605. IF(TYPE.EQ.'CHPOINT ')THEN
  606. CALL LICHT(IZDP,IZDDP,TYPE,IGEOMI)
  607. IGEOM=KGEOC
  608.  
  609. SEGACT IZD*MOD
  610. SEGACT IGEOM,IGEOMI
  611.  
  612. NPT=IZD.PROG(/1)
  613. NPTI=IZDDP.VPOCHA(/1)
  614.  
  615. IF(KSTO.EQ.0) THEN
  616. SEGACT IA*MOD,JA*MOD,IZA*MOD
  617. ELSE
  618. SEGACT KA*MOD,IZA*MOD
  619. ENDIF
  620. DO 51 I=1,NPTI
  621. XVALI=IZDDP.VPOCHA(I,1)
  622. IPOI=IGEOMI.NUM(1,I)
  623. DO 52 J=1,NPT
  624. IF(IGEOM.NUM(1,J).EQ.IPOI)THEN
  625. IZD.PROG(J)=IZD.PROG(J)+XVALI
  626. C
  627. C*** Il faut aussi mettre a jour la matrice
  628. C
  629. IF(KSTO.EQ.0) THEN
  630. IZA.PROG(IA.LECT(J))=IZA.PROG(IA.LECT(J))+XVALI
  631. ELSE
  632. IZA.PROG(J)=IZA.PROG(J)+XVALI
  633. ENDIF
  634.  
  635. GO TO 53
  636. ENDIF
  637. 52 CONTINUE
  638. 53 CONTINUE
  639. 51 CONTINUE
  640. IF(KSTO.EQ.0) THEN
  641. SEGDES IA,JA,IZA
  642. ELSE
  643. SEGDES KA,IZA
  644. ENDIF
  645. SEGDES IZDP,IZDDP,IGEOM,IGEOMI
  646. ENDIF
  647. C
  648. C********************************************************************
  649. C--- P R E C O N D I T I O N N E M E N T
  650. C********************************************************************
  651. C
  652.  
  653. C write(6,*)' P R E C O N'
  654. IF(KTYPI.EQ.2) THEN
  655. SEGACT IZDPR*MOD
  656. CALL INITD(IZDPR.PROG,NL,1.D0)
  657. CALL ECMO(MTAB3,'IPRC','LISTREEL',IZDPR)
  658. SEGDES IZDPR
  659. ELSEIF(KTYPI.EQ.3) THEN
  660. SEGACT IZDPR*MOD
  661. DO 23 K=1,NL
  662. IZDPR.PROG(K)=1./SQRT(IZD.PROG(K))
  663. 23 CONTINUE
  664. CALL ECMO(MTAB3,'IPRC','LISTREEL',IZDPR)
  665. SEGDES IZDPR
  666. ELSEIF(KTYPI.EQ.4) THEN
  667. IF(KSTO.EQ.0) THEN
  668. SEGACT IZICCG*MOD,IZA*MOD,IA*MOD,JA*MOD
  669. CALL ICGFAM(IZICCG.PROG,IZA.PROG,IA.LECT,JA.LECT,NL)
  670. SEGDES IZICCG,IZA,IA,JA
  671. ENDIF
  672. IF(KSTO.EQ.1) THEN
  673. SEGACT IZICCG*MOD,IZA*MOD,KA*MOD
  674. CALL ICGFAC(IZICCG.PROG,IZA.PROG,KA.LECT,NL,NNZ)
  675. SEGDES IZICCG,IZA,KA
  676. ENDIF
  677. CALL ECMO(MTAB3,'IPRC','LISTREEL',IZICCG)
  678. ELSE
  679. WRITE(6,*) 'CGPR Option inconnue'
  680. END IF
  681. SEGSUP IZD
  682. SEGSUP IZIPAD
  683. SEGDES MPOVAL,IZDV
  684. KCAC=1
  685. SEGDES MATRAK
  686. SEGDES MELEME
  687. C
  688. C*** RESOLUTION ***
  689. C
  690. 500 CONTINUE
  691. CALL LICHT(MCHB,IZB,TYPC,IGEOM)
  692. CALL LICHT(MCHR,IZR,TYPC,IGEOM)
  693. C write(6,*)' AVT KAL2P '
  694. CALL KAL2P(MTABP,MCHB)
  695. C write(6,*)' AVT RESCGS'
  696. CALL RESCGS(IZB,MTABP,IZR)
  697. C write(6,*)' APR RESCGS'
  698. segact melstb
  699. RETURN
  700. 90 CONTINUE
  701. WRITE(6,*)' Arret anormal de KRESS'
  702. RETURN
  703. END
  704.  
  705.  
  706.  
  707.  
  708.  
  709.  
  710.  
  711.  
  712.  
  713.  
  714.  

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