1 : C RLEXVB SOURCE CHAT 05/01/13 03:03:11 5004 2 : SUBROUTINE RLEXVB(MELFL,MELFP,MELSOM,MLELEM) 3 : C 4 : C**** Variables de COOPTIO 5 : C 6 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 7 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 8 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 9 : C & ,IECHO, IIMPI, IOSPI 10 : C & ,IDIM 11 : CC & ,MCOORD 12 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 13 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 14 : C & ,NORINC,NORVAL,NORIND,NORVAD 15 : C & ,NUCROU, IPSAUV 16 : C 17 : IMPLICIT INTEGER(I-N) 18 : INTEGER NSOMM, LAST, NTOTN, NTOTE, NBSOUS, NBELEM, NBNO 19 : & , IELEM, NGF, NGF1, INOEU, NGS1, NLS1, ISOUS 20 : & , IPOS1, IPOSV1, NVOIS1, INOEU2, NGS2, NLS2 21 : & , IVOIS, IPOS, IELEMF, I1 22 : LOGICAL LOGCON 23 : C 24 : -INC SMCOORD 25 : -INC SMELEME 26 : INTEGER JG 27 : -INC SMLENTI 28 : -INC CCOPTIO 29 : C 30 : INTEGER NBL, NBTPOI 31 : SEGMENT MLELEM 32 : INTEGER INDEX(NBL+1) 33 : INTEGER LESPOI(NBTPOI) 34 : ENDSEGMENT 35 : C 36 : POINTEUR MELSOM.MELEME, MELFL.MELEME, MELFP.MELEME 37 : & ,MLESOM.MLENTI, MLEFP.MLENTI, MPOS.MLENTI, MTOUC.MLENTI 38 : & ,MLELE1.MLELEM, MELFP1.MELEME 39 : C 40 : C**** Le MELEME SOMMET 41 : C 41 : CALL KRIPAD(MELSOM,MLESOM) 42 : C 43 : C MLESOM: numerotation globale -> locale 44 : C 45 : C**** En KRIPAD 46 : C SEGACT MELSOM 47 : C SEGINI MLESOM 48 : C 49 : NSOMM = MELSOM.NUM(/2) 50 : JG=NSOMM 51 : SEGINI MTOUC 52 : C MTOUC.LECT(NLS1) = surestimation de nombre de voisins de NLS1 53 : SEGINI MPOS 54 : C MPOS.LECT + LAST : liste chaînée des sommets sur le bord 55 : LAST=-1 56 : C 57 : SEGACT MELFP 58 : C 59 : C**** En 2D MELFP est un maillage elementaire 60 : C En 3D pas à priori 61 : C 62 : C NTOTE : nombre total d'elts dans la même structure 63 : C 64 : NTOTN = 0 65 : NTOTE = 0 66 : NBSOUS=MELFP.LISOUS(/1) 67 : C NBSOUS=0 fais un peux chier! 68 : JG=MAX(NBSOUS,1) 69 : SEGINI MLEFP 70 : IF(NBSOUS .EQ. 0)THEN 71 : MLEFP.LECT(1)=MELFP 72 : ELSE 73 : DO ISOUS=1,NBSOUS,1 74 : MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS) 75 : ENDDO 76 : ENDIF 77 : SEGDES MELFP 78 : C 79 : SEGACT MELFL 80 : NBSOUS=MELFL.LISOUS(/1) 81 : IF(NBSOUS .NE. 0)THEN 82 : WRITE(IOIMP,*) 'FACEL = ???' 83 : WRITE(IOIMP,*) 'Subroutine rlexvb.eso' 83 : CALL ERREUR(5) 84 : GOTO 9999 85 : ENDIF 86 : C 87 : NBSOUS=JG 88 : IELEMF=0 89 : DO ISOUS = 1, NBSOUS, 1 90 : MELFP1=MLEFP.LECT(ISOUS) 91 : SEGACT MELFP1 92 : NBELEM=MELFP1.NUM(/2) 93 : NBNO=MELFP1.NUM(/1) - 1 94 : DO IELEM = 1, NBELEM,1 95 : IELEMF=IELEMF+1 96 : NGF=MELFP1.NUM(NBNO+1,IELEM) 97 : NGF1=MELFL.NUM(2,IELEMF) 98 : IF(NGF .NE. NGF1)THEN 99 : WRITE(IOIMP,*) 'FACEL = ???' 100 : WRITE(IOIMP,*) 'Subroutine rlexvb.eso' 100 : CALL ERREUR(5) 101 : GOTO 9999 102 : ENDIF 103 : IF(MELFL.NUM(1,IELEMF) .EQ. MELFL.NUM(3,IELEMF))THEN 104 : C 105 : C************* La face est sur le bord 106 : C 107 : DO INOEU = 1, NBNO, 1 108 : NGS1 = MELFP1.NUM(INOEU,IELEM) 109 : NLS1 = MLESOM.LECT(NGS1) 110 : IF(MPOS.LECT(NLS1).EQ.0)THEN 111 : NTOTE=NTOTE+1 112 : MPOS.LECT(NLS1)=LAST 113 : LAST=NLS1 114 : ENDIF 115 : ENDDO 116 : ENDIF 117 : ENDDO 118 : ENDDO 119 : C 120 : C NTOTN : nombre total de noeuds dans la structure sommets de bords 121 : C + sommets voisins 122 : C N.B. Les sommets voisins peuvent ne pas etre sur le bord 123 : C 124 : NTOTN = 0 125 : DO ISOUS = 1, NBSOUS, 1 126 : MELFP1=MLEFP.LECT(ISOUS) 127 : NBELEM=MELFP1.NUM(/2) 128 : NBNO=MELFP1.NUM(/1) - 1 129 : DO IELEM = 1, NBELEM,1 130 : DO INOEU = 1, NBNO, 1 131 : NGS1 = MELFP1.NUM(INOEU,IELEM) 132 : NLS1 = MLESOM.LECT(NGS1) 133 : IF(MPOS.LECT(NLS1).NE.0)THEN 134 : C 135 : C******************** Le sommet est sur le bord 136 : C 137 : MTOUC.LECT(NLS1)=MTOUC.LECT(NLS1)+NBNO -1 138 : NTOTN=NTOTN+(NBNO - 1) 139 : ENDIF 140 : ENDDO 141 : ENDDO 142 : ENDDO 143 : C 144 : C**** On stoke l'informaton dans une LISTE SEQUENTIELLE INDEXEE D'ELEMENTS 145 : C 146 : C NBL : NOMBRE D'ELEMENTS 147 : C NBTPOI : NOMBRE TOTAL DE POINTS REFERENCEES 148 : C INDEX(I) : INDICE DU 1ER POINT DU IEME ELEMENT 149 : C DANS LE TABLEAU LESPOI 150 : C LESPOI(INDEX(I) -> INDEX(I+1)-1) : NUMERO DES NOEUDS 151 : C DU IEME ELEMENT 152 : C 153 : C En plus, après le REPEAT UNTIL, 154 : C MPOS.LECT(NLS1) = position du sommet NLS1 i.e. 155 : C J=MPOS.LECT(NLS1) -> NLS1=LESPOI(INDEX(J)) 156 : C se voisins sont dedans LESPOI(INDEX(J)+1 -> INDEX(I+1)-1) 157 : NBL=NTOTE 158 : NBTPOI=NBL+NTOTN 159 : SEGINI MLELEM 160 : IPOS1=0 161 : MLELEM.INDEX(1)=1 162 : C 163 : C**** REPEAT UNTIL 164 : C 165 : 1 CONTINUE 166 : NLS1 = LAST 167 : IF(NLS1 .NE. -1)THEN 168 : IPOS1=IPOS1+1 169 : MLELEM.LESPOI(MLELEM.INDEX(IPOS1))=NLS1 170 : LAST = MPOS.LECT(NLS1) 171 : MPOS.LECT(NLS1)=IPOS1 172 : MLELEM.INDEX(IPOS1+1)=MLELEM.INDEX(IPOS1)+MTOUC.LECT(NLS1)+1 173 : MTOUC.LECT(NLS1)=0 174 : GOTO 1 175 : ENDIF 176 : C 177 : C**** N.B. MTOUC.LECT(I) = 0 \forall I 178 : C 179 : C Fin REPEAT UNTIL 180 : C On remplie LESPOI 181 : C 182 : DO ISOUS = 1, NBSOUS, 1 183 : MELFP1 = MLEFP.LECT(ISOUS) 184 : NBELEM=MELFP1.NUM(/2) 185 : NBNO=MELFP1.NUM(/1) - 1 186 : DO IELEM = 1, NBELEM, 1 187 : DO INOEU = 1, NBNO, 1 188 : NGS1 = MELFP1.NUM(INOEU,IELEM) 189 : NLS1 = MLESOM.LECT(NGS1) 190 : C IPOS1 = 0 -> NGS1 n'est pas sur le bords 191 : C IPOS1 = C > 0 192 : C C= position de NGS1 (NLS1) dedans MLELEM.INDEX 193 : IPOS1 = MPOS.LECT(NLS1) 194 : IF(IPOS1 .GT. 0)THEN 195 : C 196 : C IPOSV1 = position du premier voisin de NLS1 dedans 197 : C MLEMEM.LESPOI 198 : IPOSV1 = MLELEM.INDEX(IPOS1)+1 199 : C NVOIS1 = surestimation du nombre de voisins de NLS1 200 : NVOIS1 = MLELEM.INDEX(IPOS1+1)-IPOSV1 201 : DO INOEU2 = 1, NBNO, 1 202 : NGS2=MELFP1.NUM(INOEU2,IELEM) 203 : IF(NGS1 .NE. NGS2)THEN 204 : NLS2 = MLESOM.LECT(NGS2) 205 : IVOIS=0 206 : C 207 : C**************** REPEAT UNTIL 208 : C 209 : 10 CONTINUE 210 : IF(MLELEM.LESPOI(IPOSV1+IVOIS) .EQ. 0)THEN 211 : MTOUC.LECT(NLS1)=MTOUC.LECT(NLS1)+1 212 : MLELEM.LESPOI(IPOSV1+IVOIS) = NLS2 213 : ELSE 214 : LOGCON = (NLS2 .NE. MLELEM.LESPOI(IPOSV1 215 : $ +IVOIS)) 216 : IVOIS=IVOIS+1 217 : IF(LOGCON)THEN 218 : IF(IVOIS .GT. NVOIS1)THEN 219 : C 220 : C************************* Il y a un erreur de programmation (NVOIS1 n'a 221 : C pas bien été évalué) 222 : C 223 : WRITE(IOIMP,*) 'Subroutine rlexvb.eso' 223 : CALL ERREUR(5) 224 : GOTO 9999 225 : ELSE 226 : GOTO 10 227 : ENDIF 228 : ENDIF 229 : ENDIF 230 : ENDIF 231 : ENDDO 232 : ENDIF 233 : ENDDO 234 : ENDDO 235 : SEGDES MELFP1 236 : ENDDO 237 : SEGSUP MLEFP 238 : C 239 : C**** LESPOI SURDIMENSIONÉ 240 : C 241 : NBTPOI=0 242 : DO IELEM = 1 , NTOTE, 1 243 : NLS1 = MLELEM.LESPOI(MLELEM.INDEX(IELEM)) 244 : NBTPOI=NBTPOI+1+MTOUC.LECT(NLS1) 245 : ENDDO 246 : NBL=NTOTE 247 : SEGINI MLELE1 248 : C 249 : MLELE1.INDEX(1)=1 250 : DO IELEM = 1, NBL , 1 251 : IPOS=MLELEM.INDEX(IELEM) 252 : IPOS1=MLELE1.INDEX(IELEM) 253 : NLS1 = MLELEM.LESPOI(IPOS) 254 : NGS1 = MELSOM.NUM(1,NLS1) 255 : NVOIS1 = MTOUC.LECT(NLS1) 256 : MLELE1.LESPOI(IPOS1)=NGS1 257 : MLELE1.INDEX(IELEM+1)=MLELE1.INDEX(IELEM)+NVOIS1+1 258 : DO INOEU = 1, NVOIS1, 1 259 : NLS2 = MLELEM.LESPOI(IPOS+INOEU) 260 : IF(NLS2 .EQ.0)THEN 261 : WRITE(IOIMP,*) 'Subroutine rlexvb.eso' 261 : CALL ERREUR(5) 262 : GOTO 9999 263 : ELSE 264 : NGS2=MELSOM.NUM(1,NLS2) 265 : MLELE1.LESPOI(IPOS1+INOEU)=NGS2 266 : ENDIF 267 : ENDDO 268 : ENDDO 269 : C 270 : SEGSUP MLELEM 271 : MLELEM=MLELE1 272 : C 273 : SEGSUP MLESOM 274 : SEGSUP MPOS 275 : SEGSUP MTOUC 276 : SEGDES MLELEM 277 : SEGDES MELSOM 278 : SEGDES MELFL 279 : C 280 : C**** Problème détecté dans un cas test: mal-conditionnent 281 : C de la matrice associé à la méthode de moindres carres 282 : C aux bords. 283 : C Quand possible, on elimine de MLELEM tous les voisins 284 : C sommets qui sont sur le bord -> MLELE1 285 : C 286 : SEGACT MLELEM 287 : NBL=MLELEM.INDEX(/1)-1 288 : NBTPOI=MLELEM.LESPOI(/1) 289 : SEGINI MLELE1 290 : C 291 : C**** MLESOM.MLENTI contient les sommets aux bords 292 : C 293 : JG=MCOORD.XCOOR(/1)/(IDIM+1) 294 : SEGINI MLESOM 295 : DO IELEM=1,NBL,1 296 : IPOS=MLELEM.INDEX(IELEM) 297 : NGS1=MLELEM.LESPOI(IPOS) 298 : MLESOM.LECT(NGS1)=IELEM 299 : ENDDO 300 : C 301 : C 302 : C**** On rempli MLELE1 303 : C 304 : NBTPOI=0 305 : NVOIS1=0 306 : C 307 : C**** Pour chaque SOMMET NGS1, NVOIS1 = nombre des voisins qui ne 308 : C sont pas sur le bords 309 : C 310 : IPOS1=MLELEM.INDEX(1) 311 : DO IELEM=1,NBL,1 312 : IPOS=IPOS1 313 : IPOS1=MLELEM.INDEX(1+IELEM) 314 : NBTPOI=NBTPOI+1 315 : MLELE1.INDEX(IELEM)=NBTPOI 316 : NGS1=MLELEM.LESPOI(IPOS) 317 : MLELE1.LESPOI(NBTPOI)=NGS1 318 : DO I1=IPOS+1,IPOS1-1,1 319 : NGS2=MLELEM.LESPOI(I1) 320 : NLS2=MLESOM.LECT(NGS2) 321 : IF(NLS2.EQ.0)THEN 322 : C 323 : C************* NGS2 n'est pas sur le bord 324 : C 325 : NVOIS1=NVOIS1+1 326 : MLELE1.LESPOI(NBTPOI+NVOIS1)=NGS2 327 : ENDIF 328 : ENDDO 329 : IF(NVOIS1.EQ.0)THEN 330 : C 331 : C********** Tous les voisins de NGS1 sont sur le bords 332 : C 333 : DO I1=IPOS+1,IPOS1-1,1 334 : NGS2=MLELEM.LESPOI(I1) 335 : NVOIS1=NVOIS1+1 336 : MLELE1.LESPOI(NBTPOI+NVOIS1)=NGS2 337 : ENDDO 338 : ENDIF 339 : NBTPOI=NBTPOI+NVOIS1 340 : NVOIS1=0 341 : ENDDO 342 : MLELE1.INDEX(NBL+1)=NBTPOI+1 343 : C 344 : SEGSUP MLELEM 345 : SEGADJ MLELE1 346 : MLELEM=MLELE1 347 : SEGDES MLELEM 348 : SEGSUP MLESOM 349 : C 350 : 9999 RETURN 351 : END 352 : 353 :
© Cast3M 2003 - All rights reserved.
Disclaimer