1 : C RLECA1 SOURCE CHAT 05/01/13 03:01:02 5004 2 : SUBROUTINE RLECA1(MELFL,MELFP,MLEFSC,MACOE1) 3 : C************************************************************************ 4 : C 5 : C PROJET : CASTEM 2000 6 : C 7 : C NOM : RLECA1 8 : C 9 : C DESCRIPTION : Appelle par GRADIA 10 : C 11 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) 12 : C 13 : C AUTEUR : A. BECCANTINI 14 : C 15 : C************************************************************************ 16 : C 17 : C Inputs: 18 : C 19 : C MELFL : MELEME FACEL 20 : C 21 : C MELFP : MELEME FACEP 22 : C 23 : C Outputs: 24 : C 25 : C MLEFSC : list of FACES points and their neighbors (CENTRE and SOMMET 26 : C points) 27 : C 28 : C MACOE1 : coeff. for linear exact reconstruction of MLEFSC 29 : C 30 : C 31 : C**** Variables de COOPTIO 32 : C 33 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 34 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 35 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 36 : C & ,IECHO, IIMPI, IOSPI 37 : C & ,IDIM 38 : CC & ,MCOORD 39 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 40 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 41 : C & ,NORINC,NORVAL,NORIND,NORVAD 42 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI 43 : C 44 : IMPLICIT INTEGER(I-N) 45 : -INC CCOPTIO 46 : -INC SMCOORD 47 : -INC SMELEME 48 : POINTEUR MELFL.MELEME,MELFP.MELEME,MELFP1.MELEME 49 : INTEGER JG 50 : -INC SMLENTI 51 : POINTEUR MLEFP.MLENTI 52 : -INC SMLREEL 53 : POINTEUR MLRSIG.MLREEL, MLRTRA.MLREEL 54 : & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL 55 : & ,MLRIN3.MLREEL 56 : C 57 : INTEGER NBL, NBTPOI 58 : SEGMENT MLELEM 59 : INTEGER INDEX(NBL+1) 60 : INTEGER LESPOI(NBTPOI) 61 : ENDSEGMENT 62 : POINTEUR MLEFSC.MLELEM 63 : C 64 : INTEGER N1,N2 65 : SEGMENT MATRIX 66 : REAL*8 MAT(N1,N2) 67 : ENDSEGMENT 68 : POINTEUR MATA.MATRIX,MATU.MATRIX,MATV.MATRIX 69 : & ,MACOE1.MATRIX 70 : LOGICAL LOGSVD 71 : C 72 : INTEGER NBSOUS, NVMAX, NELT, NNOE, ISOUS, NVMIN, NVOI, IELEM 73 : & , IPOS, IELEM1,NGF,NGF1, NGN, IPCOOR, I1, I2, I3, IPVOI 74 : & , IERSVD, NGCD, NGCG, IVOI,IERR0,I4 75 : 76 : REAL*8 ERRTOL,SMAX,SMIN, XF, YF, ZF 77 : PARAMETER (ERRTOL=1.0D0-15) 78 : C 79 : NBL = 0 80 : NBTPOI = 0 81 : NVMAX = 0 82 : SEGACT MELFP 83 : NBSOUS=MELFP.LISOUS(/1) 84 : C NBSOUS=0 fais un peux chier! 85 : JG=MAX(NBSOUS,1) 86 : SEGINI MLEFP 87 : IF(NBSOUS .EQ. 0)THEN 88 : MLEFP.LECT(1)=MELFP 89 : ELSE 90 : DO ISOUS=1,NBSOUS,1 91 : MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS) 92 : ENDDO 93 : ENDIF 94 : SEGDES MELFP 95 : NBSOUS=JG 96 : C 97 : DO ISOUS=1,NBSOUS,1 98 : MELFP1=MLEFP.LECT(ISOUS) 99 : SEGACT MELFP1 100 : NNOE=MELFP1.NUM(/1)-1 101 : NELT=MELFP1.NUM(/2) 102 : NBTPOI=NBTPOI+((1+2+NNOE)*NELT) 103 : C 104 : C 1 = centre de face 105 : C 2 = centre d'elts "gauche" et "droite" (si differents) 106 : C NNOE = sommet voisins 107 : C 108 : NBL=NBL+NELT 109 : NVMAX=MAX(NVMAX,(2+NNOE)) 110 : ENDDO 111 : C 112 : C N.B.: NBTPOI surestimé 113 : C 114 : C**** MLEFSC 115 : C Structure Face-(Sommets-Centres) 116 : C 117 : SEGINI MLEFSC 118 : C 119 : C**** La matrice de coeff. 120 : C 121 : C MACOE1(*,I) = coefficients concernants MLEFSC.LESPOI(I) 122 : C 123 : NVMIN=IDIM+1 124 : C 125 : C**** NVMIN = nombre minimum de voisins au dessus duquel MATA est 126 : C singulier 127 : C 128 : N1=NVMIN 129 : N2=NBTPOI 130 : SEGINI MACOE1 131 : C 132 : C**** MATA = matrice à inverser (NVMAX,IDIM+1) 133 : C MATU = matice des vectors propres singuliers droite de MATA 134 : C (NVMAX,IDIM+1) 135 : C MATV = matrice des vectors propres singuliers gauche de MATA 136 : C (IDIM+1,IDIM+1) 137 : C Mais en INVSVD a dimension NVMAX,IDIM+1 138 : C MLRSIG =liste des valeurs singuliers de MATA 139 : C 140 : C N.B. MATA = MATU MATSIG t(MATV) 141 : C Si MATA non singulier 142 : C inv(MATA) = MATV inv(MATSIG) t(MATU) 143 : C 144 : N1=NVMAX 145 : N2=NVMIN 146 : SEGINI MATA 147 : SEGINI MATU 148 : SEGINI MATV 149 : JG=NVMIN 150 : SEGINI MLRSIG 151 : SEGINI MLRTRA 152 : C MLRTRA segment de travail de la subroutine invsvd 153 : C 154 : C MTAA : [t(MATA).MATA] 155 : C MINTAA : [inve(t(MATA) MATA)] 156 : C MLRIN1,2,3 : "temporary vectors" 157 : C 158 : N1=NVMIN 159 : N2=NVMIN 160 : SEGINI MTAA 161 : SEGINI MINTAA 162 : JG=NVMIN 163 : SEGINI MLRIN1 164 : SEGINI MLRIN2 165 : SEGINI MLRIN3 166 : C 167 : SEGACT MELFL 168 : MLEFSC.INDEX(1)=1 169 : IELEM=0 170 : DO ISOUS=1,NBSOUS,1 171 : MELFP1=MLEFP.LECT(ISOUS) 172 : NNOE=MELFP1.NUM(/1)-1 173 : NELT=MELFP1.NUM(/2) 174 : DO IELEM1=1,NELT,1 175 : IELEM=IELEM+1 176 : IPOS=MLEFSC.INDEX(IELEM) 177 : NGF=MELFP1.NUM(NNOE+1,IELEM1) 178 : NGF1=MELFL.NUM(2,IELEM) 179 : IF(NGF .NE. NGF1)THEN 180 : WRITE(IOIMP,*) 'subroutine rleca1.eso' 180 : CALL ERREUR(5) 181 : GOTO 9999 182 : ENDIF 183 : IPCOOR=(IDIM+1)*(NGF-1)+1 184 : XF=MCOORD.XCOOR(IPCOOR) 185 : YF=MCOORD.XCOOR(IPCOOR+1) 186 : IF(IDIM .EQ. 3) ZF=MCOORD.XCOOR(IPCOOR+2) 187 : MLEFSC.LESPOI(IPOS)=NGF 188 : DO I1 = 1, NVMIN, 1 189 : MACOE1.MAT(I1,IPOS)=0.0D0 190 : ENDDO 191 : C 192 : C********** Le centre "gauche" 193 : C 194 : NGCG=MELFL.NUM(1,IELEM) 195 : NGN=NGCG 196 : C NVOI=nombre de voisins 197 : NVOI=1 198 : MLEFSC.LESPOI(IPOS+1)=NGN 199 : MATA.MAT(1,1)=1.0D0 200 : IPCOOR=(IDIM+1)*(NGN-1)+1 201 : MATA.MAT(1,2)=MCOORD.XCOOR(IPCOOR)-XF 202 : MATA.MAT(1,3)=MCOORD.XCOOR(IPCOOR+1)-YF 203 : IF(IDIM.EQ.3) MATA.MAT(1,4)=MCOORD.XCOOR(IPCOOR+2)-ZF 204 : C 205 : C********** Le centre "droite" 206 : C 207 : NGCD=MELFL.NUM(3,IELEM) 208 : IF(NGCD .NE. NGCG)THEN 209 : NVOI=NVOI+1 210 : NGN=NGCD 211 : MLEFSC.LESPOI(IPOS+2)=NGN 212 : MATA.MAT(2,1)=1.0D0 213 : IPCOOR=(IDIM+1)*(NGN-1)+1 214 : MATA.MAT(2,2)=MCOORD.XCOOR(IPCOOR)-XF 215 : MATA.MAT(2,3)=MCOORD.XCOOR(IPCOOR+1)-YF 216 : IF(IDIM.EQ.3) MATA.MAT(2,4)=MCOORD.XCOOR(IPCOOR+2)-ZF 217 : ENDIF 218 : C 219 : C********** Les sommets 220 : C 221 : DO I1 = 1, NNOE, 1 222 : NGN=MELFP1.NUM(I1,IELEM1) 223 : MLEFSC.LESPOI(IPOS+NVOI+I1)=NGN 224 : MATA.MAT(NVOI+I1,1)=1.0D0 225 : IPCOOR=(IDIM+1)*(NGN-1)+1 226 : MATA.MAT(NVOI+I1,2)=MCOORD.XCOOR(IPCOOR)-XF 227 : MATA.MAT(NVOI+I1,3)=MCOORD.XCOOR(IPCOOR+1)-YF 228 : IF(IDIM.EQ.3) MATA.MAT(NVOI+I1,4)=MCOORD.XCOOR(IPCOOR+2) 229 : & -ZF 230 : ENDDO 231 : NVOI=NVOI+NNOE 232 : MLEFSC.INDEX(IELEM+1)=IPOS+NVOI+1 233 : C 234 : LOGSVD=.TRUE. 234 : CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT, 235 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD, 236 : & MLRTRA.PROG) 237 : IF(IERSVD.NE.0)THEN 238 : LOGSVD=.FALSE. 239 : ELSE 240 : SMAX=0.0D0 241 : DO I2=1,NVMIN,1 242 : SMAX=MAX(SMAX,MLRSIG.PROG(I2)) 243 : ENDDO 244 : SMIN=SMAX 245 : DO I2=1,NVMIN,1 246 : SMIN=MIN(SMIN,MLRSIG.PROG(I2)) 247 : ENDDO 248 : ENDIF 249 : IF((SMIN/SMAX) .LT. ERRTOL)THEN 250 : LOGSVD=.FALSE. 251 : ENDIF 252 : C 253 : C********** We have two different cases: 254 : C LOGSVD = .TRUE. -> we have the coeff. we are looking for 255 : C LOGSVD = .FALSE. -> we have not 256 : C 257 : C LOGSVD=.FALSE. 258 : IF(LOGSVD)THEN 259 : DO I4=1,NVMIN,1 260 : DO I3=1,NVOI,1 261 : IPVOI=IPOS+I3 262 : DO I2=1,NVMIN,1 263 : MACOE1.MAT(I2,IPVOI)=MACOE1.MAT(I2,IPVOI)+ 264 : & (MATV.MAT(I2,I4)*MATU.MAT(I3,I4)/ 265 : & MLRSIG.PROG(I4)) 266 : ENDDO 267 : ENDDO 268 : ENDDO 269 : ELSE 270 : WRITE (IOIMP,*) 'rleca1.eso' 271 : C 22 0 272 : C Opération malvenue. Résultat douteux 272 : CALL ERREUR(22) 273 : C 274 : C************* INVSVD does not worked 275 : C For each neighbor k we have to compute the solution 276 : C of 277 : C 278 : C t(MATA) MATA x = t(MATA) * b 279 : C 280 : C where b= \sum_l e_l \delta(k,l) = e_k 281 : C 282 : C To do that, we compute 283 : C 284 : C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 285 : C 286 : C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 287 : C [t(MATA) MATA] X_0] 288 : C 289 : C 290 : C************* We compute [t(MATA) MATA] 291 : C We store it in the upper triangle of MTAA(NVMIN,NVMIN) 292 : C 293 : DO I1=1,NVMIN,1 294 : DO I2=I1,NVMIN,1 295 : MTAA.MAT(I1,I2)=0.0D0 296 : ENDDO 297 : ENDDO 298 : C 299 : DO I1=1,NVMIN,1 300 : DO I2=I1,NVMIN,1 301 : DO I3=1,NVOI,1 302 : MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+ 303 : & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2)) 304 : ENDDO 305 : ENDDO 306 : ENDDO 307 : C 308 : C************* We compute [inve(t(MATA) MATA)] 309 : C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN) 310 : C 310 : CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT, 311 : & MLRIN2.PROG,IERR0) 312 : IF(IERR0.NE.0)THEN 313 : WRITE(IOIMP,*) 'subroutine rleca1.eso.' 314 : C 26 2 315 : C Tache impossible. Probablement données erronées 315 : CALL ERREUR(26) 316 : GOTO 9999 317 : ENDIF 318 : C 319 : C************* We complete MTAA and MINTAA 320 : C 321 : DO I1=1,NVMIN,1 322 : DO I2=I1+1,NVMIN,1 323 : MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2) 324 : MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2) 325 : ENDDO 326 : ENDDO 327 : C 328 : DO IVOI=1,NVOI,1 329 : C 330 : C**************** We compute [t(MATA) . b] and we store it in MLRIN1.PROG 331 : C 332 : DO I1=1,NVMIN,1 333 : MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1) 334 : MLRIN2.PROG(I1)=0.0D0 335 : MLRIN3.PROG(I1)=0.0D0 336 : ENDDO 337 : C 338 : C**************** X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 339 : C X_0(i1) into MLRIN2.PROG(I1) 340 : C 341 : DO I2=1,NVMIN,1 342 : DO I1=1,NVMIN,1 343 : MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+ 344 : & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2)) 345 : ENDDO 346 : ENDDO 347 : C 348 : C**************** X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 349 : C [t(MATA) MATA] X_0] 350 : C 351 : C [t(MATA) MATA] X_0 into MLRIN3.PROG 352 : C 353 : DO I2=1,NVMIN,1 354 : DO I1=1,NVMIN,1 355 : MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+ 356 : & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2)) 357 : ENDDO 358 : ENDDO 359 : C 360 : C**************** Now we have 361 : C [t(MATA) . b] in MLRIN1.PROG 362 : C X_0(i1) in MLRIN2.PROG(I1) 363 : C [t(MATA) MATA] X_0 in MLRIN3.PROG 364 : C 365 : C X_1(i1) in MLRCOE.MAT(i1,IVOI) 366 : C 367 : IPVOI=IPOS+IVOI 368 : DO I1=1,NVMIN,1 369 : DO I2=1,NVMIN,1 370 : MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+ 371 : & (MINTAA.MAT(I1,I2)* 372 : & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2))) 373 : ENDDO 374 : MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+ 375 : & MLRIN2.PROG(I1) 376 : ENDDO 377 : ENDDO 378 : ENDIF 379 : C 380 : C********** If we are here, we have the coeff. in the 'FACE' frame. 381 : C 382 : ENDDO 383 : SEGDES MELFP1 384 : ENDDO 385 : C 386 : NBTPOI=MLEFSC.INDEX(NBL+1)-1 387 : SEGADJ MLEFSC 388 : N1=NVMIN 389 : N2=NBTPOI 390 : SEGADJ MACOE1 391 : C 392 : SEGSUP MLEFP 393 : SEGDES MLEFSC 394 : SEGDES MACOE1 395 : SEGSUP MATA 396 : SEGSUP MATU 397 : SEGSUP MATV 398 : SEGSUP MLRSIG 399 : SEGSUP MLRTRA 400 : C 401 : SEGDES MELFL 402 : SEGSUP MTAA 403 : SEGSUP MINTAA 404 : SEGSUP MLRIN1 405 : SEGSUP MLRIN2 406 : SEGSUP MLRIN3 407 : C 408 : 9999 RETURN 409 : END 410 : 411 : 412 : 413 : 414 : 415 :
© Cast3M 2003 - All rights reserved.
Disclaimer