1 : C PROGCS SOURCE CHAT 05/01/13 02:31:16 5004 2 : SUBROUTINE PROGCS 3 : C************************************************************************ 4 : C 5 : C CE SOUS PROGRAMME RANGE DANS LES TABLES MATAP, MATAP.METHODE ET 6 : C MATAP.MATRIS LES POINTEURS SUR LES LISTENTI D'ADRESSAGE. 7 : C 8 : C KTYPI= 2,3,4 METHODE DE GRADIENT CONJUGUE AVEC CALCUL 9 : C EXPLICITE DE LA MATRICE DE PRESSION. DANS CE CAS 10 : C PROGCS CONSTRUIT LES TABLEAUX D'INDEXAGE EN 11 : C FONCTION DU MODE DE STOCKAGE DEMANDE : MORSE OU 12 : C COMPRESSE. 13 : C 14 : C SYNTAXE : MTABM = PROGCS MTABPIMPR 15 : C 16 : C************************************************************************ 17 : IMPLICIT INTEGER(I-N) 18 : IMPLICIT REAL*8 (A-H,O-Z) 19 : CHARACTER*8 TYPE 20 : 21 : C*** 22 : 23 : -INC CCOPTIO 24 : C************************************************************************* 25 : C 26 : C REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES puis assemblees 27 : C 28 : 29 : * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange 30 : * (points CENTRE ) pour chaque operateur de contrainte 31 : * KGEOC SPG pour la totalite des points CENTRE. 32 : * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse) 33 : * KLEMC Connectivites de l'ensemble des contraintes 34 : * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones 35 : 36 : SEGMENT MATRAK 37 : INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP) 38 : INTEGER LIZAFM(NBSOUS) 39 : INTEGER IKAM0 (NBSOUS) 40 : INTEGER IMEM (NBELC) 41 : INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC 42 : ENDSEGMENT 43 : 44 : SEGMENT IZAFM 45 : REAL*8 AM(NNELP,NP,IESP),RPGI(NELAX) 46 : ENDSEGMENT 47 : 48 : POINTEUR IPMJ.IZAFM,IPMK.IZAFM 49 : 50 : C******************************************************************* 51 : -INC SMCOORD 52 : -INC SMTABLE 53 : POINTEUR MTABP.MTABLE 54 : POINTEUR MTABM.MTABLE 55 : -INC SMLENTI 56 : POINTEUR IZIPAD.MLENTI 57 : POINTEUR KA.MLENTI 58 : POINTEUR IA.MLENTI 59 : POINTEUR ITAB.MLENTI 60 : -INC SMELEME 61 : POINTEUR IZK.MELEME 62 : 63 : C 64 : C IVK = Tableau de travail 65 : PARAMETER (MAXIVK=100) 66 : INTEGER IVK(MAXIVK) 67 : C 68 : C MTABP : TABLE MATAP 69 : C MTABM : SOUS-TABLE MATRICE 70 : C MTAB2 : SOUS-TABLE METHODE 71 : C 72 : CHARACTER*4 LMOT(1) 73 : 74 : PARAMETER (NTB=1) 75 : CHARACTER*8 LTAB(NTB) 76 : DIMENSION KTAB(NTB) 77 : DATA LTAB/'EQPR '/ 78 : DATA LMOT/'IMPR'/ 79 : C*** 80 : 81 : NTO=1 81 : CALL LITABS(LTAB,KTAB,NTB,NTO,IRET) 82 : IF(IRET.EQ.0)THEN 83 : WRITE(6,*)' On attend une table soustype EQPR' 84 : RETURN 85 : ENDIF 86 : 87 : IMPR=0 87 : CALL LIRMOT(LMOT,1,IP,0) 88 : IF(IP.NE.0)THEN 88 : CALL LIRENT(IMPR,1,IRET) 89 : IF(IRET.EQ.0)RETURN 90 : ENDIF 91 : 92 : MTABP=KTAB(1) 93 : SEGACT MTABP 94 : 95 : TYPE=' ' 95 : CALL ACMO(MTABP,'MATC',TYPE,MATRAK) 96 : IF(TYPE.NE.'MATRAK')THEN 97 : WRITE(6,*)' Il n''y a pas d''entree MATRAK dans la table' 98 : RETURN 99 : ENDIF 100 : 101 : SEGACT MATRAK 102 : KGEO=KGEOC 102 : CALL KNBEL(KGEO,NBELC) 103 : MELEME=KLEMC 104 : 105 : C--- CALNOE renvoie dans IZTKB un tableau contenant pour chaque noeud 106 : C du domaine la liste des {l{ments poss{dant ce noeud. Les noeuds 107 : C sont rep{r{s dans la num{rotation locale(domaniale), les {l{ments 108 : C sont pris dans l'ordre fourni par la num{rotation naturelle. 109 : C Il est possible, partant de ce tableau, de construire pour chaque 110 : C {l{ment la liste des {l{ments voisins. Ceci est tr}s utile pour 111 : C le calcul de la matrice de pression ou Akl.ne.0 si k et l sont 112 : C deux {l{ments voisins. 113 : C 114 : C Cette technique n'est pas (encore) utilis{e si KTYPI = 1 ou 5. 115 : C Les maillages trait{s dans ce cas {tant "petits", ce n'est pas 116 : C trop grave. Par contre lorsque KTYPI=2,3,4, c'est @ dire pour les 117 : C tr}s gros maillages, son emploi est obligatoire. 118 : C 119 : C Cet appel ne doit etre fait que si on stocke la matrice de 120 : C pression (KTPI.LE.5). 121 : C 121 : CALL KALNO0(MELEME,IZK,MLENTI,IZIPAD,IRET) 122 : IF(IRET.EQ.0)GO TO 90 123 : C 124 : C Creation de la table MATRIS 125 : C 125 : CALL CRTABL(MTABM) 125 : CALL ECMM(MTABM,'SOUSTYPE','MATRIS') 126 : 127 : C--- Activation de tous les MELEMEs du domaine 128 : 129 : SEGACT MELEME,IZK,MLENTI 130 : NBSOUS=LISOUS(/1) 131 : IF(NBSOUS.EQ.0)NBSOUS=1 132 : DO 400 NS=1,NBSOUS 133 : IF(NBSOUS.EQ.1)IPT1=MELEME 134 : IF(NBSOUS.NE.1)IPT1=LISOUS(NS) 135 : SEGACT IPT1 136 : 400 CONTINUE 137 : 138 : * METHODE 139 : TYPE=' ' 139 : CALL ACMO(MTABP,'METHODE',TYPE,MTAB2) 140 : SEGACT MTAB2 141 : * KTYPI 141 : CALL ACME(MTAB2,'KTYPI',KTPI) 142 : * KSTOCK 142 : CALL ACME(MTAB2,'KSTOCK',KSTOCK) 143 : C 144 : C---- CALCUL DES TABLEAUX DE VOISINAGE 145 : C 146 : IF(KSTOCK.EQ.0) THEN 147 : JG=NBELC+1 148 : SEGINI IA 149 : ENDIF 150 : C 151 : C---- Le premier travail consiste a fabriquer la liste des voisins de 152 : C chaque element a partir du travail effectue par CALNOE. 153 : C On a la liste des elements auquel appartient chaque noeud. Il faut 154 : C consolider ces informations pour tous les noeuds. Si un noeud 155 : C appartient a l'element k et a l'element l, k et l sont voisins. 156 : C On remplit une table IVK par element. Les tableaux IVK sont 157 : C reunis dans un segment de travail ITAB. 158 : C 159 : 160 : IF(IDIM.EQ.2) MAJOR=20 161 : IF(IDIM.EQ.3) MAJOR=40 162 : JG=MAJOR*NBELC 163 : SEGINI ITAB 164 : 165 : MAXVOI=0 166 : LA=0 167 : K=0 168 : IF(KSTOCK.EQ.0) IA.LECT(1)=1 169 : DO 401 NS=1,NBSOUS 170 : IF(NBSOUS.EQ.1)IPT1=MELEME 171 : IF(NBSOUS.NE.1)IPT1=LISOUS(NS) 172 : NP=IPT1.NUM(/1) 173 : NEL=IPT1.NUM(/2) 174 : DO 1000 KK=1,NEL 175 : K=K+1 176 : IVK(1)=1 177 : IVK(2)=K 178 : DO 101 I=1,NP 179 : IU=IZIPAD.LECT(IPT1.NUM(I,KK)) 180 : NELV=LECT(IU) 181 : DO 102 KAP=1,NELV 182 : KN=IZK.NUM(KAP,IU) 183 : C 184 : C KN appartient il deja a la liste ? Si oui sauter au KN suivant 185 : C Si non l'ajouter 186 : DO 103 J=1,IVK(1) 187 : IF(IVK(J+1).EQ.KN) GOTO 102 188 : 103 CONTINUE 189 : IVK(1)=IVK(1)+1 190 : IVK(IVK(1)+1)=KN 191 : 102 CONTINUE 192 : 101 CONTINUE 193 : 194 : C--- LA contient le nombre d'elements de la matrice en morse 195 : C--- MAXVOI est le nombre maxi d'elements voisins d'un element donne 196 : C--- On les remet e jour. 197 : 198 : IF(IVK(1).GT.MAJOR) THEN 199 : PRINT *,'PROGCS : L''element ',KK,' du sous objet ',NS 200 : PRINT *,' i.e. ',K ,' dans la num. naturelle' 201 : PRINT *,'a un nombre de voisins',ivk(1), 202 : & 'superieur au maxi prevu' 203 : PRINT *,'Est-ce bien raisonnable ?' 204 : PRINT *,'Revoyez votre maillage |' 205 : STOP 206 : END IF 207 : IF(IVK(1).GT.MAXVOI) MAXVOI=IVK(1) 208 : LA=LA+IVK(1) 209 : 210 : C--- Le tableau IVK contient les connectivites. Il faut le reordonner 211 : C en fonction du mode de stockage demande et stocker dans IA et 212 : C ITAB pour le morse et ITAB seulement pour le mode compresse. 213 : C Pour le morse, comme pour le mode compresse on met en premier l'elt 214 : C lui-meme, puis on range dans l'ordre croissant On complete la table 215 : C jusqu'a MAJOR avec le numero de l'element dans le cas compresse. 216 : 217 : IF(KSTOCK.EQ.0) THEN 218 : C--- Cas du stockage morse 219 : IA.LECT(K+1)=IA.LECT(K)+IVK(1) 219 : C CALL ITRI(IVK(3),IVK(1)-1) 219 : CALL ORDOTA(IVK(3),IVK(1)-1) 220 : C--- On positionne dans ITAB le tableau IVK 220 : C CALL MOVMEM(IVK(2),IVK(1),ITAB.LECT(IA.LECT(K))) 220 : CALL RSETI(ITAB.LECT(IA.LECT(K)),IVK(2),IVK(1)) 221 : 222 : ELSEIF(KSTOCK.EQ.1) THEN 223 : C--- Cas du stockage compresse 224 : C--- On place en premier dans IVK le terme diagonal 225 : DO 121 I=1,IVK(1) 226 : IF(IVK(I+1).EQ.K) THEN 227 : ID=I+1 228 : GOTO 122 229 : ENDIF 230 : 121 CONTINUE 231 : 122 CONTINUE 232 : C--- On echange 232 : CALL ISWAP(IVK,2,ID,IVK(1)+1) 233 : C--- On r{ordonne le reste du tableau IVK 233 : C CALL ITRI(IVK(3),IVK(1)-1) 233 : CALL ORDOTA(IVK(3),IVK(1)-1) 234 : C--- On comble jusqu'a MAJOR avec le premier numero 235 : IF(IVK(1).LT.MAJOR) THEN 236 : DO 124 I=IVK(1)+1,MAJOR 237 : IVK(I+1)=IVK(2) 238 : 124 CONTINUE 239 : ENDIF 240 : C--- On positionne IVK dans le tableau ITAB (Indexage). ITAB 241 : C est stocke comme s'il etait declare ITAB(NBELC,MAJOR) 242 : DO 125 I=1,MAJOR 243 : IOFF=K+(I-1)*NBELC 244 : ITAB.LECT(IOFF)=IVK(I+1) 245 : 125 CONTINUE 246 : ENDIF 247 : 248 : 1000 CONTINUE 249 : SEGDES IPT1 250 : 401 CONTINUE 251 : SEGDES MELEME 252 : 253 : C--- Nous pouvons proceder a l'ajustement de la longueur du segment 254 : C ITAB : LA pour le morse, NEL*MAXVOI pour le mode compresse. 255 : 256 : IF(KSTOCK.EQ.0) JG=LA 257 : IF(KSTOCK.EQ.1) JG=NBELC*MAXVOI 258 : SEGADJ ITAB 259 : 260 : IF(KSTOCK.EQ.0) JA=ITAB 261 : IF(KSTOCK.EQ.1) KA=ITAB 262 : 263 : C! 264 : C Impression de controle 265 : C 266 : IF(IMPR.GE.2) THEN 267 : IF(KSTOCK.EQ.0) THEN 268 : WRITE(6,555) (IA.LECT(K),K=1,NBELC+1) 269 : DO 2345 K=1,NBELC 270 : WRITE(6,555) (ITAB.LECT(I),I=IA.LECT(K),IA.LECT(K+1)-1) 271 : 2345 CONTINUE 272 : ELSE 273 : DO 1234 K=1,NBELC 274 : WRITE(6,555) (ITAB.LECT(K+(I-1)*NBELC),I=1,MAXVOI) 275 : 555 FORMAT(1X,12(1X,i4)) 276 : 1234 CONTINUE 277 : ENDIF 278 : ENDIF 279 : C! 280 : SEGDES ITAB 281 : IF(KSTOCK.EQ.0) SEGDES IA 282 : IF(KSTOCK.EQ.0) LONG=LA 283 : IF(KSTOCK.EQ.1) LONG=NBELC*MAXVOI 284 : IF(IMPR.GE.1) WRITE(6,200)LONG 285 : C 285 : CALL ECME(MTABM,'NL',NBELC) 285 : CALL ECME(MTAB2,'NL',NBELC) 285 : CALL ECME(MTAB2,'NPTC',NPTC) 286 : TYPE=' ' 287 : 288 : IF(KSTOCK.EQ.0) THEN 288 : CALL ECMO(MTABM,'IA','LISTENTI',IA) 288 : CALL ECMO(MTABM,'JA','LISTENTI',JA) 288 : CALL ECME(MTABM,'LA',LA) 288 : CALL ECME(MTABM,'ILG',LA) 289 : ELSE 290 : SEGACT KA 291 : NNZ=KA.LECT(/1)/NBELC 291 : CALL ECMO(MTABM,'KA','LISTENTI',KA) 291 : CALL ECME(MTABM,'NNZ',NNZ) 291 : CALL ECME(MTABM,'ILG',NNZ) 292 : SEGDES KA 293 : ENDIF 294 : 295 : SEGSUP IZK,IZIPAD,MLENTI 296 : 296 : CALL ECROBJ('TABLE',MTABM) 297 : 298 : RETURN 299 : 90 CONTINUE 300 : WRITE(6,*)' Retour anormal de PROGCS' 301 : RETURN 302 : 200 FORMAT(10X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7, 303 : *' MOTS',13X,1H*) 304 : END 305 : 306 : 307 : 308 : 309 :
© Cast3M 2003 - All rights reserved.
Disclaimer