1 : C GRQUA4 SOURCE CHAT 05/01/13 00:21:35 5004 2 : SUBROUTINE GRQUA4(LISTPO,MLEVOI,MLEVO2,MLRDIS, 3 : & KRICEN,CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO, 4 : & MATA,MATU,MATV,MATCEL,MLRSIG,MLRCEL,COEF,INDER) 5 : C------------------------------------------------------------------ 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : GRQUA4 9 : C 10 : C DESCRIPTION : Augmente le stencil pour le calcul du gradient. 11 : C 12 : C LANGAGE : FORTRAN 77 + ESOPE 2000 13 : C 14 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 15 : C 16 : C----------------------------------------------------------------- 17 : C 18 : C APPELES (E/S) : 19 : C 20 : C APPELES (Calcul): 21 : C 22 : C----------------------------------------------------------------- 23 : C 24 : C ENTREES : LISTPO = pointeur sur une liste de point contenant 25 : C les centres et ses voisins (voir la 26 : C description du segment) 27 : C 28 : C MLEVOI = pointeur sur une liste MLENTI 29 : C (dimensions = nombre total de points) 30 : C 31 : C MLEVO2 = pointeur sur une liste MLENTI qui va contenir 32 : C les points "voisins des voisins" de centre 33 : C 34 : C MLRDIS = pointeur sur une liste MLREEL qui va contenir 35 : C les distances des points MLEVO2 du point 36 : C centre 37 : C 38 : C KRICEN = pointeur sur une liste MLENTI qui contient 39 : C la correpondance LOCAL/GLOBAL du maillage 40 : C des points centres 41 : C 42 : C CONTOL = REAL 8 t.q. 1/cond(MATA) > CONTOL 43 : C Si non rank(MATA) < NVOIMI - 1 44 : C 45 : C NBMAX, NTOTP, NVOIMI : INTEGERs qui interviennent 46 : C dans les dimensions de 47 : C certains matrices 48 : C N.B. NTOTP = E/S 49 : C 50 : C 51 : C MATA = pointeur d'une matrice 52 : C (NBMAXx(NBMAX-1),NVOIMI-1) qui contient la 53 : C matrice (NTOTP-1,NVOIMI-1), a pseudoinverser 54 : C 55 : C MATU = pointeur d'une matrice 56 : C (NBMAXx(NBMAX-1),NVOIMI-1) on met la matrice 57 : C (NTOTP-1,NVOIMI-1) des vectors singuliers 58 : C droites de MATA 59 : C 60 : C MATV = pointeur de la matrice (NVOIMI-1,NVOIMI-1) 61 : C des vectors singuliers gauches de MATA 62 : C 63 : C MATCEL = pointeur sur une matrice de travail 64 : C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient 65 : C momentanément b 66 : C la matrice (NTOTP-1,NVOIMI-1), transposé 67 : C de la pseudo-inverse de MATA 68 : C 69 : C MLRSIG = pointeur sur une liste de reels (NVOIMI-1) 70 : C qui contient lesvaleurs singuliers de MATA 71 : C 72 : C MLRCEL = pointeur sur une liste de reels 73 : C (NBMAXx(NBMAX-1)) (vecteur de travail) 74 : C 75 : C E/S : NTOTP = nombre total des voisins du point centre 76 : C dont on veut augmenter le stencil 77 : C 78 : C MLEVTO = pointeur sur un list d'entiers 79 : C (NBMAX - 1) * NBMAX : 80 : C MLEVTO.LECT(1) = point ou je veux calculer 81 : C gradient et hessien 82 : C MLEVTO.LECT(i), i=2, NTOTP = voisins 83 : C 84 : C SORTIES : COEF = pointeur sur la matrice qui contient le 85 : C coefficient pour le calcul de l'hessien 86 : C 87 : C INDER = 1 -> probleme mal posé 88 : C = 2 -> erreur anormale 89 : C 90 : C----------------------------------------------------------------- 91 : C 92 : C HISTORIQUE (Anomalies et modifications eventuelles) 93 : C 94 : C HISTORIQUE : 11.05.2000, Creation 95 : C 96 : C----------------------------------------------------------------- 97 : C 98 : C 99 : IMPLICIT INTEGER(I-N) 100 : -INC CCOPTIO 101 : -INC SMCOORD 102 : -INC SMELEME 103 : -INC SMLENTI 104 : -INC SMLREEL 105 : C 106 : C----------------------------------------------------------------- 107 : C 108 : POINTEUR MLEVOI.MLENTI, MLEVO2.MLENTI, MLRDIS.MLREEL, 109 : & KRICEN.MLENTI, MLEVTO.MLENTI, MLRSIG.MLREEL, 110 : & MLRCEL.MLREEL 111 : C 112 : SEGMENT MATRIX 113 : REAL*8 MAT(N1,N2) 114 : ENDSEGMENT 115 : C 116 : POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX, MATCEL.MATRIX, 117 : & COEF.MATRIX 118 : C 119 : SEGMENT LISTPO 120 : INTEGER LISPO1(NPOILI),INDEX1(NBCEN+1) 121 : ENDSEGMENT 122 : C 123 : C**** LISTPO: list de points centres avec leur voisins 124 : C et position: 125 : C 126 : C LISPO1(k) t.q. 127 : C LISPO1(INDEX1(j)) = numero global du centre j 128 : C LISPO1(INDEX1(j)),...,LISPO1(INDEX1(j+1))-1 = voisins du centre 129 : C LISPO1(INDEX1(j)) 130 : 131 : C 132 : C 133 : C**** Variables de COOPTIO 134 : C 135 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 136 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 137 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 138 : C & ,IECHO, IIMPI, IOSPI 139 : C & ,IDIM 140 : CC & ,MCOORD 141 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 142 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 143 : C & ,NORINC,NORVAL,NORIND,NORVAD 144 : C & ,NUCROU, IPSAUV 145 : C 146 : INTEGER NVOIMI, NTOTP, INDER, INEUVO, NOCENE, LAST0, IVOIS, NOCEN1 147 : & , ICEN, IPOS1, NVOIS1, NOCEN2, ICOORX, IVOIS2, IVOIS3, IPOS0 148 : & , NBMAX 149 : 150 : C 151 : REAL*8 XCEN, YCEN, ZCEN, DIST, DIST1, CONTOL 152 : C 153 : C INEUVO = nombre de nouveaux voisins 154 : C 155 : INEUVO = 0 156 : C 157 : INDER = 0 158 : NOCENE = MLEVTO.LECT(1) 159 : C 160 : ICOORX = ((IDIM + 1) * (NOCENE - 1))+1 161 : XCEN = MCOORD.XCOOR(ICOORX) 162 : YCEN = MCOORD.XCOOR(ICOORX+1) 163 : IF(IDIM .EQ. 3) ZCEN = MCOORD.XCOOR(ICOORX+2) 164 : C 165 : LAST0 = NOCENE 166 : MLEVOI.LECT(NOCENE) = -1 167 : C 168 : C**** On mets -1 dans le dernier elt de la liste chaîné 169 : C -1 != 0 -> on sait que NOCENE fait partir de l'elt de 170 : C noeuds + voisins qu'on veut realizer! 171 : C 172 : DO IVOIS = 2, NTOTP, 1 173 : NOCEN1 = MLEVTO.LECT(IVOIS) 174 : MLEVOI.LECT(NOCEN1) = LAST0 175 : LAST0 = NOCEN1 176 : ENDDO 177 : C 178 : C**** Liste chainé du noeud et de ses principaux voisins remplie! 179 : C 180 : NOCEN1 = LAST0 181 : C 182 : C**** Si NOCEN1 = NOCENE, ça veut dire que le pouvre est tout 183 : C seul! C'est tres triste! 184 : C 185 : IF(NOCEN1 .EQ. NOCENE) THEN 186 : MOTERR(1:40) = 'Elt isolé???????????????????????????????' 187 : WRITE(IOIMP,*) MOTERR(1:40) 188 : INDER = 1 189 : GOTO 9999 190 : ENDIF 191 : C 192 : C*** On parcoure la liste chaîné dans un sens et on la rallonge 193 : C de nouveaux voisins dans l'autre 194 : C 195 : 2 CONTINUE 196 : ICEN = KRICEN.LECT(NOCEN1) 197 : IF(ICEN .NE. 0)THEN 198 : C 199 : C******* NOCEN1 est un point centre 200 : C 201 : IPOS1 = LISTPO.INDEX1(ICEN) 202 : NVOIS1 = LISTPO.INDEX1(ICEN+1) - IPOS1 203 : DO IVOIS = 1, NVOIS1 - 1, 1 204 : NOCEN2 = LISTPO.LISPO1(IPOS1+IVOIS) 205 : IF(MLEVOI.LECT(NOCEN2) .EQ. 0)THEN 206 : C 207 : C************** On peut le rajouter dedans MLEVOI.LECT 208 : C 209 : MLEVOI.LECT(NOCEN2)=LAST0 210 : LAST0 = NOCEN2 211 : ICOORX = ((IDIM + 1) * (NOCEN2 - 1))+1 212 : DIST = (MCOORD.XCOOR(ICOORX) - XCEN)**2 213 : DIST = DIST + ((MCOORD.XCOOR(ICOORX+1)- YCEN)**2) 214 : IF(IDIM .EQ. 3) 215 : & DIST = DIST + ((MCOORD.XCOOR(ICOORX+2)- ZCEN)**2) 216 : C 217 : C************* Je cherche la position ou il faut le mettre 218 : C 219 : DO IVOIS2 = 1, INEUVO, 1 220 : DIST1 = MLRDIS.PROG(IVOIS2) 221 : IF(DIST .LT. DIST1) GOTO 3 222 : ENDDO 223 : IVOIS2 = INEUVO + 1 224 : 3 CONTINUE 225 : DO IVOIS3 = INEUVO, IVOIS2, -1 226 : MLEVO2.LECT(IVOIS3+1) = MLEVO2.LECT(IVOIS3) 227 : MLRDIS.PROG(IVOIS3+1) = MLRDIS.PROG(IVOIS3) 228 : ENDDO 229 : MLEVO2.LECT(IVOIS2) = NOCEN2 230 : MLRDIS.PROG(IVOIS2) = DIST 231 : INEUVO = INEUVO + 1 232 : ENDIF 233 : ENDDO 234 : ENDIF 235 : NOCEN1 = MLEVOI.LECT(NOCEN1) 236 : C 237 : C**** Le dernier de la liste chaîné est NOCENE 238 : C 239 : IF(NOCEN1 .NE. NOCENE) GOTO 2 240 : C 241 : C 242 : C**** Il faut nettoyer MLEVOI 243 : C 244 : NOCEN1 = LAST0 245 : 5 CONTINUE 246 : NOCEN2 = MLEVOI.LECT(NOCEN1) 247 : MLEVOI.LECT(NOCEN1) = 0 248 : NOCEN1 = NOCEN2 249 : IF(NOCEN1 .NE. -1) GOTO 5 250 : C 251 : IF(INEUVO .LT. (NVOIMI - NTOTP))THEN 252 : MOTERR(1:40) = 'Pas assez de noueds voisins ' 253 : WRITE(IOIMP,*) MOTERR(1:40) 254 : INDER = 1 255 : GOTO 9999 256 : ENDIF 257 : C 258 : C**** Test 1: la liste ordonnée 259 : C 260 : C write(*,*) ' ' 261 : C write(*,*) 'NOCENE =', NOCENE 262 : C write(*,*) 'INEUVO =', INEUVO 263 : C do ivois3 = 1, ineuvo, 1 264 : C write(*,*) mlevo2.lect(ivois3), mlrdis.prog(ivois3) 265 : C enddo 266 : C 267 : C**** Maintenant il faut selectioner les noeuds au ajouter 268 : C 269 : C 270 : IPOS1 = MAX(NVOIMI - NTOTP +1,1) 271 : IF(IPOS1 .GT. INEUVO)THEN 272 : MOTERR(1:40) = 'Pas assez de noueds voisins ' 273 : WRITE(IOIMP,*) MOTERR(1:40) 274 : INDER = 1 275 : GOTO 9999 276 : ENDIF 277 : C 278 : C**** REPEAT UNTIL IVOIS < (INEUVO+1) .AND. INDER = 1 279 : C 280 : IPOS0 = 1 281 : 10 CONTINUE 282 : DIST = MLRDIS.PROG(IPOS1) * 1.05D0 283 : DO IVOIS = IPOS1+1, INEUVO, 1 284 : DIST1 = MLRDIS.PROG(IVOIS) 285 : IF(DIST1 .GT. DIST) GOTO 4 286 : ENDDO 287 : IVOIS = INEUVO + 1 288 : C 289 : C IVOIS-1 est le dernier qui je prends 290 : C 291 : 4 CONTINUE 292 : C 293 : C IPOS0 est le premier qui n'appartient pas 294 : C MLEVTO.LECT 295 : C 296 : DO IVOIS2 = IPOS0, IVOIS - 1, 1 297 : NTOTP = NTOTP + 1 298 : MLEVTO.LECT(NTOTP)=MLEVO2.LECT(IVOIS2) 299 : ENDDO 300 : IF(IVOIS .EQ. (INEUVO+1))THEN 301 : C 302 : C MLEVTO.LECT contient déjà tous les noueds 303 : C dedans MLEVO2.LECT 304 : C 304 : CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, 305 : & MATCEL,MLRSIG,MLRCEL,COEF,INDER) 306 : C 307 : IF(INDER .EQ. 1)THEN 308 : C 309 : C**** Test 2: Tous les voisins 310 : C 311 : MOTERR(1:40) = 'Probleme dans la inversion de matrices ' 312 : WRITE(IOIMP,*) MOTERR(1:40) 313 : ENDIF 314 : ELSEIF(IVOIS .LT. (INEUVO+1))THEN 314 : CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, 315 : & MATCEL,MLRSIG,MLRCEL,COEF,INDER) 316 : IF(INDER .EQ. 1)THEN 317 : IPOS0 = IVOIS 318 : IPOS1 = IVOIS 319 : GOTO 10 320 : ELSEIF(INDER .EQ. 2)THEN 321 : GOTO 9999 322 : ENDIF 323 : ELSEIF(IVOIS .GT. (INEUVO+1))THEN 324 : C 325 : C Erreur dans la subroutine 326 : C 326 : CALL ERREUR(5) 327 : ENDIF 328 : C 329 : 330 : 9999 RETURN 331 : END 332 : 333 : 334 : 335 : 336 : 337 : 338 :
© Cast3M 2003 - All rights reserved.
Disclaimer