Télécharger profch.eso

Retour à la liste

Numérotation des lignes :

  1. C PROFCH SOURCE PV 16/11/17 22:01:15 9180
  2. SUBROUTINE PROFCH(MATRAK,IMPR,KPIMP,NIMP)
  3. C************************************************************************
  4. C
  5. C OBJET :
  6. C CE SOUS PROGRAMME PREPARE L'INVERSION DE LA MATRICE DE PRESSION
  7. C
  8. C KTYPI= 1 OU 5, METHODE DIRECTE
  9. C
  10. C 1/ IOPTI = 0 PAS DE RENUMEROTATION
  11. C IOPTI = 1 RENUMEROTATION DE LA MATRICE
  12. C DE PRESSION
  13. C 2/ CALCUL DU PROFIL ET DE LA TAILLE
  14. C DE LA DEMI-MATRICE
  15. C
  16. C MATRAK : objet de type MATRAK
  17. C Si KIZCL=0 On cree IZL
  18. C Sinon rien
  19. C IMPR 1 ou 2 niveau d'impressions
  20. C KPIMP = 0 pas de pression imposee
  21. C KPIMP = 1 pression imposée en un point (le dernier)
  22. C KPIMP = 2 pression moyenne imposée
  23. C NIMP numero de la pression imposée
  24. C
  25. C************************************************************************
  26. IMPLICIT INTEGER(I-N)
  27. IMPLICIT REAL*8 (A-H,O-Z)
  28.  
  29. C***
  30. C-INC SMMATRAKANC
  31. C*************************************************************************
  32. C
  33. C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees
  34. C
  35.  
  36. * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange
  37. * (points CENTRE ) pour chaque operateur de contrainte
  38. * KGEOC SPG pour la totalite des points CENTRE.
  39. * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse)
  40. * KLEMC Connectivites de l'ensemble des contraintes
  41. * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones
  42.  
  43. SEGMENT MATRAK
  44. INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP)
  45. INTEGER LIZAFM(NBSOUS)
  46. INTEGER IKAM0 (NBSOUS)
  47. INTEGER IMEM (NBELC)
  48. INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC
  49. ENDSEGMENT
  50.  
  51. SEGMENT IZAFM
  52. REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX)
  53. ENDSEGMENT
  54.  
  55. POINTEUR IPMJ.IZAFM,IPMK.IZAFM
  56.  
  57. C*******************************************************************
  58. -INC SMCOORD
  59. -INC SMLENTI
  60. POINTEUR IZIPAD.MLENTI
  61. -INC SMELEME
  62. POINTEUR IZK.MELEME
  63. -INC SMTABLE
  64. POINTEUR MTABP.MTABLE
  65. PARAMETER (IBLSIZ=50000)
  66.  
  67. SEGMENT IZTRAV
  68. INTEGER ITAB(NJA)
  69. ENDSEGMENT
  70. C
  71. C IVK = Tableau de travail
  72. PARAMETER (MAXIVK=100)
  73. INTEGER IVK(MAXIVK)
  74. C***
  75. CHARACTER*8 TYPE
  76. CHARACTER*4 LMOT(1)
  77. DATA LMOT/'IMPR'/
  78. DATA IOPTI/1/
  79. C***
  80.  
  81. SEGACT MATRAK*MOD
  82. IF(KIZCL.NE.0)THEN
  83. SEGDES MATRAK
  84. RETURN
  85. ENDIF
  86.  
  87. KGEO=KGEOC
  88. CALL KNBEL(KGEO,NBELC)
  89. MELEME=KLEMC
  90.  
  91. C
  92. C--- CALNOE renvoie dans IZTKB un tableau contenant pour chaque noeud
  93. C du domaine la liste des {l{ments poss{dant ce noeud. Les noeuds
  94. C sont rep{r{s dans la num{rotation locale(domaniale), les {l{ments
  95. C sont pris dans l'ordre fourni par la num{rotation naturelle.
  96. C Il est possible, partant de ce tableau, de construire pour chaque
  97. C {l{ment la liste des {l{ments voisins. Ceci est tr}s utile pour
  98. C le calcul de la matrice de pression ou Akl.ne.0 si k et l sont
  99. C deux {l{ments voisins.
  100. C
  101. C Cette technique n'est pas (encore) utilis{e si KTYPI = 1 ou 5.
  102. C Les maillages trait{s dans ce cas {tant "petits", ce n'est pas
  103. C trop grave. Par contre lorsque KTYPI=2,3,4, c'est @ dire pour les
  104. C tr}s gros maillages, son emploi est obligatoire.
  105. C
  106. C Cet appel ne doit etre fait que si on stocke la matrice de
  107. C pression (KTPI.LE.5).
  108. C
  109. CALL KALNO0(MELEME,IZK,MLENTI,IZIPAD,IRET)
  110. IF(IRET.EQ.0)GO TO 90
  111.  
  112. C--- Activation de tous les MELEMEs du domaine
  113.  
  114. SEGACT MELEME,IZK,MLENTI
  115. NBSOUS=LISOUS(/1)
  116. IF(NBSOUS.EQ.0)NBSOUS=1
  117. DO 400 NS=1,NBSOUS
  118. IF(NBSOUS.EQ.1)IPT1=MELEME
  119. IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  120. SEGACT IPT1
  121. 400 CONTINUE
  122.  
  123. C
  124. C************************************************************************
  125. C C H O L E S K I
  126. C************************************************************************
  127. C
  128. NJA=NBELC
  129. NJAN=NBELC
  130. NJAB=NBELC
  131.  
  132. SEGINI IZL,IZTRAV
  133. KIZCL=IZL
  134.  
  135. LONG=0
  136.  
  137. CALL INITI(NUAN,NBELC,0)
  138. CALL INITI(NUNA,NBELC,0)
  139. NUAN(1)=1
  140. NUNA(1)=1
  141.  
  142. C********* ON STOKE LES POINTEURS DES MELEME
  143.  
  144. KI=1
  145. DO 40 NS=1,NBSOUS
  146. IF(NBSOUS.EQ.1)IPT1=MELEME
  147. IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  148.  
  149. NEL=IPT1.NUM(/2)
  150. IPOINT=IPT1
  151. CALL INITI(IMEL(KI),NEL,IPOINT)
  152.  
  153. DO 41 I=1,NEL
  154. IMJ(KI+I-1)=I
  155. 41 CONTINUE
  156.  
  157. KI=KI+NEL
  158.  
  159. 40 CONTINUE
  160.  
  161. C*********
  162.  
  163. IF(IOPTI.EQ.1)GO TO 33
  164.  
  165. C********* ON NE RENUMEROTE PAS
  166.  
  167. IF(IMPR.GE.2)WRITE(6,*)' PAS DE RENUMEROTATION DE LA PRESSION '
  168. DO 32 I=1,NBELC
  169. NUAN(I)=I
  170. NUNA(I)=I
  171. 32 CONTINUE
  172. GO TO 34
  173.  
  174. C********* OPTIMISATION DE LA NUMEROTAION DE LA MATRICE DE PRESSION
  175. 33 CONTINUE
  176. IF(IMPR.GE.2)WRITE(6,*)' ON RENUMEROTE LA MATRICE DE PRESSION '
  177. C*********
  178.  
  179. C write(6,*)' SUB PROFCH : IPADL :'
  180. C nbp=ipadl(/1)
  181. C write(6,1001)(IPADL(II),II=1,NBP)
  182. C im2=lect(/1)
  183. C im1=izk.num(/1)
  184. C write(6,*)' IM1=',IM1 ,' Taille IZK.NUM'
  185. C do 1873 i=1,im2
  186. C nbb=lect(i)
  187. C write(6,*)' I=->',i,' nb=',lect(i)
  188. C write(6,1001)(izk.num(ii,i),ii=1,nbb)
  189. C1873 continue
  190.  
  191. KC=1
  192.  
  193. DO 30 KK=1,NBELC
  194.  
  195. K=NUNA(KK)
  196.  
  197. IF(K.NE.0)GO TO 310
  198. C ON A PLUS D'ELEMENT NOUVELLEMENT NUMEROTE
  199. WRITE(6,*)' EST CE LA DISCONTINUITE ',K,KC,KK
  200. KC=KC+1
  201. DO 311 KKA=1,NBELC
  202. IF(NUAN(KKA).EQ.0)GO TO 312
  203. C SI NUAN(KKA)=0 C'EST QUE L'ELEMENT KKA N'A PAS ETE NUMEROTE
  204. 311 CONTINUE
  205. WRITE(6,*)' IL YA UN PROBLEME DE NUMEROTATION'
  206. CALL ARRET(0)
  207. 312 CONTINUE
  208. NUAN(KKA)=KC
  209. NUNA(KC)=KKA
  210. K=NUNA(KC)
  211. 310 CONTINUE
  212.  
  213. IPT1=IMEL(K)
  214. NP=IPT1.NUM(/1)
  215. DO 31 I=1,NP
  216. NU=IPT1.NUM(I,IMJ(K))
  217. C NU=IPT1.NUM(I,K)
  218. NU=IZIPAD.LECT(NU)
  219. NELV=LECT(NU)
  220. DO 31 NN=1,NELV
  221. KN=IZK.NUM(NN,NU)
  222.  
  223. IF(NUAN(KN).NE.0)GO TO 31
  224. KC=KC+1
  225. NUAN(KN)=KC
  226. NUNA(KC)=KN
  227. 31 CONTINUE
  228.  
  229. 30 CONTINUE
  230.  
  231. DO 35 I=1,NBELC
  232. ITAB(I)=IMJ(NUNA(I))
  233. 35 CONTINUE
  234. CALL RSETI(IMJ,ITAB,NBELC)
  235.  
  236. DO 36 I=1,NBELC
  237. ITAB(I)=IMEL(NUNA(I))
  238. 36 CONTINUE
  239. CALL RSETI(IMEL,ITAB,NBELC)
  240.  
  241. SEGSUP IZTRAV
  242. C **********************************************
  243. 34 CONTINUE
  244.  
  245. IF(IMPR.GE.2)THEN
  246. WRITE(6,1929)
  247. 1929 FORMAT(/10X,' SUB PROFIM : RENUMEROTAION ZONE PRESSION ')
  248. WRITE(6,*)' CORRESPONDANCE NUMER. ANC. --> NUMER. NOUVELLE '
  249. WRITE(6,1930)(N,NUAN(N),N=1,NBELC)
  250. WRITE(6,*)' CORRESPONDANCE NUMER. NOUVELLE --> NUMER. ANC. '
  251. WRITE(6,1930)(N,NUNA(N),N=1,NBELC)
  252. 1930 FORMAT(7(1X,'**',I5,3X,I5,1X,'**'))
  253. ENDIF
  254. C
  255. DO 10 KK=1,NBELC
  256. K=IMJ(KK)
  257. IPT1=IMEL(KK)
  258. NP=IPT1.NUM(/1)
  259. K0=KK
  260. DO 11 I=1,NP
  261. NU=IPT1.NUM(I,K)
  262. NU=IZIPAD.LECT(NU)
  263. NELV=LECT(NU)
  264. DO 11 NN=1,NELV
  265. KN=IZK.NUM(NN,NU)
  266. KN=NUAN(KN)
  267. IF(KN.LT.K0)K0=KN
  268. 11 CONTINUE
  269. IF(KK.NE.NIMP.OR.KPIMP.NE.2)THEN
  270. KZA(KK)=KK-K0+1
  271. ELSE
  272. KZA(KK)=NBELC-1
  273. ENDIF
  274. LONG=LONG+KZA(KK)
  275. 10 CONTINUE
  276.  
  277. C
  278. C La matrice est stockée par blocs de plusieurs lignes. En effet, il faut
  279. C l'équivalent de 900 opérations pour activer/désactiver un segment. Si
  280. C une ligne de la matrice comprend 50 éléments, on fait 100 opérations
  281. C utiles et 900 inutiles lors de la résolution.
  282. C La taille des blocs est fixée par le paramètre IBLSIZ. La diagonale est
  283. C stockée dans un segment à part.
  284. C
  285. C Il nous faut constituer un certain nombre de tableaux. Le premier segment
  286. C décrit la matrice. Il contient un tableau donnant les numéros des premières
  287. C lignes de chaque bloc et un tableau de pointeurs sur des segments
  288. C descripteurs de bloc.
  289. C
  290. C Il faut d'abord déterminer le nombre de blocs.
  291. C
  292.  
  293. NBLK=1
  294. IF(LONG-NBELC.GT.IBLSIZ) THEN
  295. ILBL=0
  296. DO 12 KK=2,NBELC
  297. LA=KZA(KK)-1
  298. IF(LA.GT.IBLSIZ) THEN
  299. NBLK=NBLK+2
  300. ILBL=0
  301. GOTO 12
  302. ENDIF
  303. ILBL=ILBL+LA
  304. IF(ILBL.GT.IBLSIZ) THEN
  305. NBLK=NBLK+1
  306. ILBL=LA
  307. GOTO 12
  308. ENDIF
  309. 12 CONTINUE
  310. ENDIF
  311. C
  312. C On peut allouer le segment descripteur de la matrice. On va pouvoir
  313. C remplir des numéros de ligne des débuts de bloc.
  314. C
  315. SEGINI IDMAT
  316. IF(NBLK.EQ.1) THEN
  317. NLDBLK(1)=2
  318. NLDBLK(2)=NBELC+1
  319. ELSE
  320. IBLK=1
  321. ILBL=0
  322. NLDBLK(1)=2
  323. NLDBLK(NBLK+1)=NBELC+1
  324. DO 13 KK=2,NBELC
  325. LA=KZA(KK)-1
  326. IF(LA.GT.IBLSIZ) THEN
  327. NLDBLK(IBLK+1)=KK
  328. NLDBLK(IBLK+2)=KK+1
  329. IBLK=IBLK+2
  330. ILBL=0
  331. GOTO 13
  332. ENDIF
  333. ILBL=ILBL+LA
  334. IF(ILBL.GT.IBLSIZ) THEN
  335. NLDBLK(IBLK+1)=KK
  336. ILBL=LA
  337. IBLK=IBLK+1
  338. GOTO 13
  339. ENDIF
  340. 13 CONTINUE
  341. ENDIF
  342. C
  343. C On peut maintenant remplir le segment descripteur de chaque bloc.
  344. C
  345. DO 15 IBLK=1,NBLK
  346. NLBLK=NLDBLK(IBLK+1)-NLDBLK(IBLK)
  347. SEGINI IDBLK
  348. IDESCR(IBLK)=IDBLK
  349. ILON=0
  350. IDEBLK(1)=1
  351. DO 16 KK=NLDBLK(IBLK),NLDBLK(IBLK+1)-1
  352. LA=KZA(KK)-1
  353. ILON=ILON+LA
  354. IDEBLK(KK-NLDBLK(IBLK)+2)=IDEBLK(KK-NLDBLK(IBLK)+1)+LA
  355. 16 CONTINUE
  356. SEGDES IDBLK
  357. 15 CONTINUE
  358. SEGDES IDMAT
  359. C
  360. C On stocke dans KZA1 la référence au segment IDMAT à partir duquel tout se
  361. C retrouve.
  362. C
  363. C write(6,*)' sub PROFCH on met IDMAT -> KZA1 ',IDMAT
  364. KZA1=IDMAT
  365. 14 CONTINUE
  366. C
  367. IF(IMPR.GE.2)WRITE(6,*)(KZA(M),M=1,NBELC),LONG
  368. C
  369. IF(IMPR.GE.1)WRITE(6,200)LONG
  370. SEGDES IZL
  371. SEGSUP IZK,IZIPAD,MLENTI
  372.  
  373. C--- D{sactivation de tous les MELEMEs
  374.  
  375. DO 51 NS=1,NBSOUS
  376. IF(NBSOUS.EQ.1)IPT1=MELEME
  377. IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  378. SEGDES IPT1
  379. 51 CONTINUE
  380. SEGDES MELEME
  381.  
  382. RETURN
  383. 90 CONTINUE
  384. WRITE(6,*)' Retour anormal de PROFCH'
  385. RETURN
  386.  
  387. C************************************************************************
  388. 100 FORMAT(40X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7,
  389. *' MOTS',13X,1H*/40X,1H*,72X,1H*)
  390. 200 FORMAT(10X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7,
  391. *' MOTS',13X,1H*)
  392. 1001 FORMAT(20(1X,I5))
  393. END
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  

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