1 : C PROFCH SOURCE CHAT 05/01/13 02:30:30 5004 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 SMMATRK1 59 : -INC SMCOORD 60 : -INC SMLENTI 61 : POINTEUR IZIPAD.MLENTI 62 : -INC SMELEME 63 : POINTEUR IZK.MELEME 64 : -INC SMTABLE 65 : POINTEUR MTABP.MTABLE 66 : PARAMETER (IBLSIZ=50000) 67 : 68 : SEGMENT IZTRAV 69 : INTEGER ITAB(NJA) 70 : ENDSEGMENT 71 : C 72 : C IVK = Tableau de travail 73 : PARAMETER (MAXIVK=100) 74 : INTEGER IVK(MAXIVK) 75 : C*** 76 : CHARACTER*8 TYPE 77 : CHARACTER*4 LMOT(1) 78 : DATA LMOT/'IMPR'/ 79 : DATA IOPTI/1/ 80 : C*** 81 : 82 : SEGACT MATRAK*MOD 83 : IF(KIZCL.NE.0)THEN 84 : SEGDES MATRAK 85 : RETURN 86 : ENDIF 87 : 88 : 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 108 : CALL KALNO0(MELEME,IZK,MLENTI,IZIPAD,IRET) 109 : IF(IRET.EQ.0)GO TO 90 110 : 111 : C--- Activation de tous les MELEMEs du domaine 112 : 113 : SEGACT MELEME,IZK,MLENTI 114 : NBSOUS=LISOUS(/1) 115 : IF(NBSOUS.EQ.0)NBSOUS=1 116 : DO 400 NS=1,NBSOUS 117 : IF(NBSOUS.EQ.1)IPT1=MELEME 118 : IF(NBSOUS.NE.1)IPT1=LISOUS(NS) 119 : SEGACT IPT1 120 : 400 CONTINUE 121 : 122 : C 123 : C************************************************************************ 124 : C C H O L E S K I 125 : C************************************************************************ 126 : C 127 : NJA=NBELC 128 : NJAN=NBELC 129 : NJAB=NBELC 130 : 131 : SEGINI IZL,IZTRAV 132 : KIZCL=IZL 133 : 134 : LONG=0 135 : 135 : CALL INITI(NUAN,NBELC,0) 135 : CALL INITI(NUNA,NBELC,0) 136 : NUAN(1)=1 137 : NUNA(1)=1 138 : 139 : C********* ON STOKE LES POINTEURS DES MELEME 140 : 141 : KI=1 142 : DO 40 NS=1,NBSOUS 143 : IF(NBSOUS.EQ.1)IPT1=MELEME 144 : IF(NBSOUS.NE.1)IPT1=LISOUS(NS) 145 : 146 : NEL=IPT1.NUM(/2) 147 : IPOINT=IPT1 147 : CALL INITI(IMEL(KI),NEL,IPOINT) 148 : 149 : DO 41 I=1,NEL 150 : IMJ(KI+I-1)=I 151 : 41 CONTINUE 152 : 153 : KI=KI+NEL 154 : 155 : 40 CONTINUE 156 : 157 : C********* 158 : 159 : IF(IOPTI.EQ.1)GO TO 33 160 : 161 : C********* ON NE RENUMEROTE PAS 162 : 163 : IF(IMPR.GE.2)WRITE(6,*)' PAS DE RENUMEROTATION DE LA PRESSION ' 164 : DO 32 I=1,NBELC 165 : NUAN(I)=I 166 : NUNA(I)=I 167 : 32 CONTINUE 168 : GO TO 34 169 : 170 : C********* OPTIMISATION DE LA NUMEROTAION DE LA MATRICE DE PRESSION 171 : 33 CONTINUE 172 : IF(IMPR.GE.2)WRITE(6,*)' ON RENUMEROTE LA MATRICE DE PRESSION ' 173 : C********* 174 : 175 : C write(6,*)' SUB PROFCH : IPADL :' 176 : C nbp=ipadl(/1) 177 : C write(6,1001)(IPADL(II),II=1,NBP) 178 : C im2=lect(/1) 179 : C im1=izk.num(/1) 180 : C write(6,*)' IM1=',IM1 ,' Taille IZK.NUM' 181 : C do 1873 i=1,im2 182 : C nbb=lect(i) 183 : C write(6,*)' I=->',i,' nb=',lect(i) 184 : C write(6,1001)(izk.num(ii,i),ii=1,nbb) 185 : C1873 continue 186 : 187 : KC=1 188 : 189 : DO 30 KK=1,NBELC 190 : 191 : K=NUNA(KK) 192 : 193 : IF(K.NE.0)GO TO 310 194 : C ON A PLUS D'ELEMENT NOUVELLEMENT NUMEROTE 195 : WRITE(6,*)' EST CE LA DISCONTINUITE ',K,KC,KK 196 : KC=KC+1 197 : DO 311 KKA=1,NBELC 198 : IF(NUAN(KKA).EQ.0)GO TO 312 199 : C SI NUAN(KKA)=0 C'EST QUE L'ELEMENT KKA N'A PAS ETE NUMEROTE 200 : 311 CONTINUE 201 : WRITE(6,*)' IL YA UN PROBLEME DE NUMEROTATION' 201 : CALL ARRET(0) 202 : 312 CONTINUE 203 : NUAN(KKA)=KC 204 : NUNA(KC)=KKA 205 : K=NUNA(KC) 206 : 310 CONTINUE 207 : 208 : IPT1=IMEL(K) 209 : NP=IPT1.NUM(/1) 210 : DO 31 I=1,NP 211 : NU=IPT1.NUM(I,IMJ(K)) 212 : C NU=IPT1.NUM(I,K) 213 : NU=IZIPAD.LECT(NU) 214 : NELV=LECT(NU) 215 : DO 31 NN=1,NELV 216 : KN=IZK.NUM(NN,NU) 217 : 218 : IF(NUAN(KN).NE.0)GO TO 31 219 : KC=KC+1 220 : NUAN(KN)=KC 221 : NUNA(KC)=KN 222 : 31 CONTINUE 223 : 224 : 30 CONTINUE 225 : 226 : DO 35 I=1,NBELC 227 : ITAB(I)=IMJ(NUNA(I)) 228 : 35 CONTINUE 228 : CALL RSETI(IMJ,ITAB,NBELC) 229 : 230 : DO 36 I=1,NBELC 231 : ITAB(I)=IMEL(NUNA(I)) 232 : 36 CONTINUE 232 : CALL RSETI(IMEL,ITAB,NBELC) 233 : 234 : SEGSUP IZTRAV 235 : C ********************************************** 236 : 34 CONTINUE 237 : 238 : IF(IMPR.GE.2)THEN 239 : WRITE(6,1929) 240 : 1929 FORMAT(/10X,' SUB PROFIM : RENUMEROTAION ZONE PRESSION ') 241 : WRITE(6,*)' CORRESPONDANCE NUMER. ANC. --> NUMER. NOUVELLE ' 242 : WRITE(6,1930)(N,NUAN(N),N=1,NBELC) 243 : WRITE(6,*)' CORRESPONDANCE NUMER. NOUVELLE --> NUMER. ANC. ' 244 : WRITE(6,1930)(N,NUNA(N),N=1,NBELC) 245 : 1930 FORMAT(7(1X,'**',I5,3X,I5,1X,'**')) 246 : ENDIF 247 : C 248 : DO 10 KK=1,NBELC 249 : K=IMJ(KK) 250 : IPT1=IMEL(KK) 251 : NP=IPT1.NUM(/1) 252 : K0=KK 253 : DO 11 I=1,NP 254 : NU=IPT1.NUM(I,K) 255 : NU=IZIPAD.LECT(NU) 256 : NELV=LECT(NU) 257 : DO 11 NN=1,NELV 258 : KN=IZK.NUM(NN,NU) 259 : KN=NUAN(KN) 260 : IF(KN.LT.K0)K0=KN 261 : 11 CONTINUE 262 : IF(KK.NE.NIMP.OR.KPIMP.NE.2)THEN 263 : KZA(KK)=KK-K0+1 264 : ELSE 265 : KZA(KK)=NBELC-1 266 : ENDIF 267 : LONG=LONG+KZA(KK) 268 : 10 CONTINUE 269 : 270 : C 271 : C La matrice est stockée par blocs de plusieurs lignes. En effet, il faut 272 : C l'équivalent de 900 opérations pour activer/désactiver un segment. Si 273 : C une ligne de la matrice comprend 50 éléments, on fait 100 opérations 274 : C utiles et 900 inutiles lors de la résolution. 275 : C La taille des blocs est fixée par le paramètre IBLSIZ. La diagonale est 276 : C stockée dans un segment à part. 277 : C 278 : C Il nous faut constituer un certain nombre de tableaux. Le premier segment 279 : C décrit la matrice. Il contient un tableau donnant les numéros des premières 280 : C lignes de chaque bloc et un tableau de pointeurs sur des segments 281 : C descripteurs de bloc. 282 : C 283 : C Il faut d'abord déterminer le nombre de blocs. 284 : C 285 : 286 : NBLK=1 287 : IF(LONG-NBELC.GT.IBLSIZ) THEN 288 : ILBL=0 289 : DO 12 KK=2,NBELC 290 : LA=KZA(KK)-1 291 : IF(LA.GT.IBLSIZ) THEN 292 : NBLK=NBLK+2 293 : ILBL=0 294 : GOTO 12 295 : ENDIF 296 : ILBL=ILBL+LA 297 : IF(ILBL.GT.IBLSIZ) THEN 298 : NBLK=NBLK+1 299 : ILBL=LA 300 : GOTO 12 301 : ENDIF 302 : 12 CONTINUE 303 : ENDIF 304 : C 305 : C On peut allouer le segment descripteur de la matrice. On va pouvoir 306 : C remplir des numéros de ligne des débuts de bloc. 307 : C 308 : SEGINI IDMAT 309 : IF(NBLK.EQ.1) THEN 310 : NLDBLK(1)=2 311 : NLDBLK(2)=NBELC+1 312 : ELSE 313 : IBLK=1 314 : ILBL=0 315 : NLDBLK(1)=2 316 : NLDBLK(NBLK+1)=NBELC+1 317 : DO 13 KK=2,NBELC 318 : LA=KZA(KK)-1 319 : IF(LA.GT.IBLSIZ) THEN 320 : NLDBLK(IBLK+1)=KK 321 : NLDBLK(IBLK+2)=KK+1 322 : IBLK=IBLK+2 323 : ILBL=0 324 : GOTO 13 325 : ENDIF 326 : ILBL=ILBL+LA 327 : IF(ILBL.GT.IBLSIZ) THEN 328 : NLDBLK(IBLK+1)=KK 329 : ILBL=LA 330 : IBLK=IBLK+1 331 : GOTO 13 332 : ENDIF 333 : 13 CONTINUE 334 : ENDIF 335 : C 336 : C On peut maintenant remplir le segment descripteur de chaque bloc. 337 : C 338 : DO 15 IBLK=1,NBLK 339 : NLBLK=NLDBLK(IBLK+1)-NLDBLK(IBLK) 340 : SEGINI IDBLK 341 : IDESCR(IBLK)=IDBLK 342 : ILON=0 343 : IDEBLK(1)=1 344 : DO 16 KK=NLDBLK(IBLK),NLDBLK(IBLK+1)-1 345 : LA=KZA(KK)-1 346 : ILON=ILON+LA 347 : IDEBLK(KK-NLDBLK(IBLK)+2)=IDEBLK(KK-NLDBLK(IBLK)+1)+LA 348 : 16 CONTINUE 349 : SEGDES IDBLK 350 : 15 CONTINUE 351 : SEGDES IDMAT 352 : C 353 : C On stocke dans KZA1 la référence au segment IDMAT à partir duquel tout se 354 : C retrouve. 355 : C 356 : C write(6,*)' sub PROFCH on met IDMAT -> KZA1 ',IDMAT 357 : KZA1=IDMAT 358 : 14 CONTINUE 359 : C 360 : IF(IMPR.GE.2)WRITE(6,*)(KZA(M),M=1,NBELC),LONG 361 : C 362 : IF(IMPR.GE.1)WRITE(6,200)LONG 363 : SEGDES IZL 364 : SEGSUP IZK,IZIPAD,MLENTI 365 : 366 : C--- D{sactivation de tous les MELEMEs 367 : 368 : DO 51 NS=1,NBSOUS 369 : IF(NBSOUS.EQ.1)IPT1=MELEME 370 : IF(NBSOUS.NE.1)IPT1=LISOUS(NS) 371 : SEGDES IPT1 372 : 51 CONTINUE 373 : SEGDES MELEME 374 : 375 : RETURN 376 : 90 CONTINUE 377 : WRITE(6,*)' Retour anormal de PROFCH' 378 : RETURN 379 : 380 : C************************************************************************ 381 : 100 FORMAT(40X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7, 382 : *' MOTS',13X,1H*/40X,1H*,72X,1H*) 383 : 200 FORMAT(10X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7, 384 : *' MOTS',13X,1H*) 385 : 1001 FORMAT(20(1X,I5)) 386 : END 387 : 388 : 389 :
© Cast3M 2003 - All rights reserved.
Disclaimer