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

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