1 : C RLENCF SOURCE CHAT 05/01/13 03:01:14 5004 2 : SUBROUTINE RLENCF(MELFL,MELFP,MLEPOI,MLECOE) 3 : C************************************************************************ 4 : C 5 : C PROJET : CASTEM 2000 6 : C 7 : C NOM : RLENCF 8 : C 9 : C DESCRIPTION : Appelle par GRADI2 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 MLESCF : list of SOMMET points and their CENTRE neighbors 20 : C 21 : C MATCOE : coeff. for linear exact reconstruction of MLESCF 22 : C 23 : C MLEFSC : list of FACES points and their neighbors (CENTRE and SOMMET 24 : C points) 25 : C 26 : C MACOE1 : coeff. for linear exact reconstruction of MLEFSC 27 : C 28 : C Output 29 : C 30 : C MLEFC : list of FACES points and their neighbors (CENTRE points only) 31 : C 32 : C MACOE2 : coeff. for linear exact reconstruction of MLEFC 33 : C 34 : C This subroutine computes the coefficients to compute the gradient 35 : C at intefaces with respect to the values on its neighbors. 36 : C The neighbors are 'CENTRE' points (in FACEL) ore 'SOMMET' points 37 : C (in FACEP). which allow to perform 38 : C A linear exact reconstruction is performed via a least square 39 : C method. 40 : C 41 : C**** Inputs: 42 : C 43 : C MELFL = the 'FACEL' meleme 44 : C MELFP = the 'FACEP' meleme 45 : C 46 : C**** Output: 47 : C 48 : C MLEPOI = list of integers. 49 : C MLEPOI.LECT(i) points to the list neighbors of 50 : C MELFL.NUM(2,i) 51 : C MLECOE = list of integers. 52 : C MLECOE.LECT(i) points to the matrix of coeffients to 53 : C compute the gradient 54 : C 55 : C**** Variables de COOPTIO 56 : C 57 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 58 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 59 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 60 : C & ,IECHO, IIMPI, IOSPI 61 : C & ,IDIM 62 : CC & ,MCOORD 63 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 64 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 65 : C & ,NORINC,NORVAL,NORIND,NORVAD 66 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI 67 : C 68 : IMPLICIT INTEGER(I-N) 69 : INTEGER NBSOUS,ISOUS,NBELEM,NBNO,NVMAX,NVMIN,NFAC,IFAC,IELEM 70 : & ,NGF,NGF1,NGC1,NGC2,INOEU,NGS,IPCOOR,NVOI,IVOI,NGV 71 : & ,IERSVD,IERR0 72 : & ,I1,I2,I3,I4 73 : 74 : REAL*8 XF,YF,ZF, XV, YV, ZV,SMAX,SMIN,ERRTOL 75 : PARAMETER(ERRTOL=1.0D-16) 76 : LOGICAL LOGSVD 77 : -INC CCOPTIO 78 : -INC SMCOORD 79 : -INC SMELEME 80 : -INC SMLREEL 81 : -INC SMLENTI 82 : -INC SMCHPOI 83 : INTEGER JG 84 : INTEGER N1,N2 85 : SEGMENT MATRIX 86 : REAL*8 MAT(N1,N2) 87 : ENDSEGMENT 88 : POINTEUR MELFL.MELEME, MELFP.MELEME, MLEPOI.MLENTI,MLECOE.MLENTI 89 : & , MATA.MATRIX,MATU.MATRIX,MATV.MATRIX 90 : & ,MATCOE.MATRIX,MLRSIG.MLREEL,MLRTRA.MLREEL 91 : & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL 92 : & ,MLRIN3.MLREEL,MLEFP.MLENTI,MELFP1.MELEME 93 : C 94 : SEGACT MELFP 95 : C 96 : C 97 : C**** En 2D MELFP est un maillage elementaire 98 : C En 3D pas à priori 99 : C -> MLEFP contains the list of LISOUS 100 : C 101 : NBSOUS=MELFP.LISOUS(/1) 102 : C NBSOUS=0 fais un peux chier! 103 : JG=MAX(NBSOUS,1) 104 : SEGINI MLEFP 105 : IF(NBSOUS .EQ. 0)THEN 106 : MLEFP.LECT(1)=MELFP 107 : ELSE 108 : DO ISOUS=1,NBSOUS,1 109 : MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS) 110 : ENDDO 111 : ENDIF 112 : SEGDES MELFP 113 : SEGACT MELFL 114 : NFAC=MELFL.NUM(/2) 115 : JG=NFAC 116 : SEGINI MLEPOI 117 : SEGINI MLECOE 118 : C 119 : C**** Loop on the faces to compute NVMAX 120 : C (maximum number of neighbors) 121 : C 122 : NBSOUS=MLEFP.LECT(/1) 123 : NVMAX=0 124 : DO ISOUS = 1, NBSOUS, 1 125 : MELFP1=MLEFP.LECT(ISOUS) 126 : SEGACT MELFP1 127 : NBELEM=MELFP1.NUM(/2) 128 : NBNO=MELFP1.NUM(/1) - 1 129 : C 130 : C The ISOUS elementary meleme of 'FACEP' has NBELEM elements 131 : C which contain NBNO sommets and one point face. Each face has 132 : C NBNO 'SOMMET' neighbors and 2 'CENTRE' neighbors. 133 : C 134 : NVMAX=MAX(NVMAX,NBNO+2) 135 : ENDDO 136 : C 137 : C MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1) 138 : C MATU = matrix of the singular right eigenvectors of MATA 139 : C (NVMAX,IDIM+1) 140 : C MATV = matrix of the singular left eigenvectors of MATA 141 : C (IDIM+1,IDIM+1) 142 : C But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1 143 : C MLRSIG = singular values of MATA (IDIM+1) 144 : C 145 : C N.B. MATA = MATU MATSIG t(MATV) 146 : C If MATA is non singular, 147 : C inv(MATA) = MATV inv(MATSIG) t(MATU) 148 : C 149 : C MLRTRA temporary vector of invsvd.eso 150 : C NVMIN = IDIM + 1 (the most little dimension of the matrices) 151 : C 152 : N1=NVMAX 153 : N2=IDIM+1 154 : SEGINI MATA 155 : SEGINI MATU 156 : SEGINI MATV 157 : JG=IDIM+1 158 : SEGINI MLRSIG 159 : SEGINI MLRTRA 160 : NVMIN=N2 161 : C 162 : C MTAA : [t(MATA).MATA] 163 : C MINTAA : [inve(t(MATA) MATA)] 164 : C MLRIN1,2,3 : "temporary vectors" 165 : C 166 : N1=NVMIN 167 : N2=NVMIN 168 : SEGINI MTAA 169 : SEGINI MINTAA 170 : JG=NVMIN 171 : SEGINI MLRIN1 172 : SEGINI MLRIN2 173 : SEGINI MLRIN3 174 : C 175 : C**** Loop on faces to compute the coefficient 176 : C 177 : NBSOUS=MLEFP.LECT(/1) 178 : IFAC=0 179 : DO ISOUS = 1, NBSOUS, 1 180 : MELFP1=MLEFP.LECT(ISOUS) 181 : NBELEM=MELFP1.NUM(/2) 182 : NBNO=MELFP1.NUM(/1) - 1 183 : DO IELEM=1,NBELEM,1 184 : IFAC=IFAC+1 185 : NGF=MELFP1.NUM(NBNO+1,IELEM) 186 : NGF1=MELFL.NUM(2,IFAC) 187 : NGC1=MELFL.NUM(1,IFAC) 188 : NGC2=MELFL.NUM(3,IFAC) 189 : IF(NGF .NE. NGF1)THEN 190 : WRITE(IOIMP,*) 'Subroutine rlencf.eso' 190 : CALL ERREUR(5) 191 : GOTO 9999 192 : ENDIF 193 : IF(NGC1 .EQ. NGC2)THEN 194 : JG=NBNO+1 195 : ELSE 196 : JG=NBNO+2 197 : ENDIF 198 : C 199 : C********** We create the list of neighbors. 200 : C 201 : SEGINI MLENTI 202 : MLEPOI.LECT(IFAC)=MLENTI 203 : C 204 : C********** We create the matrix of coefficients. 205 : C 206 : N2=JG 207 : N1=NVMIN 208 : SEGINI MATCOE 209 : MLECOE.LECT(IFAC)=MATCOE 210 : C 211 : C********** We fill the list of neighbors. 212 : C 213 : DO INOEU=1,NBNO,1 214 : NGS=MELFP1.NUM(INOEU,IELEM) 215 : MLENTI.LECT(INOEU)=NGS 216 : ENDDO 217 : MLENTI.LECT(NBNO+1)=NGC1 218 : IF(NGC1 .NE. NGC2) MLENTI.LECT(NBNO+2)=NGC2 219 : ENDDO 220 : SEGDES MELFP1 221 : ENDDO 222 : C 223 : C**** Test for MLEPOI 224 : C 225 : C do ifac=1,nfac,1 226 : C mlenti=mlepoi.lect(ifac) 227 : C nbno=mlenti.lect(/1) 228 : C write(*,*) 'ngf=',melfl.num(2,ifac) 229 : C write(*,*) 'nvoi=',nbno 230 : C write(*,*) (mlenti.lect(inoeu),inoeu=1,nbno,1) 231 : C enddo 232 : C 233 : C**** We have to fill the matrix coefficient 234 : C 235 : NFAC=MELFL.NUM(/2) 236 : DO IFAC=1,NFAC,1 237 : NGF=MELFL.NUM(2,IFAC) 238 : IPCOOR=((NGF-1)*(IDIM+1))+1 239 : XF=MCOORD.XCOOR(IPCOOR) 240 : YF=MCOORD.XCOOR(IPCOOR+1) 241 : IF(IDIM .EQ. 3) ZF=MCOORD.XCOOR(IPCOOR+2) 242 : C 243 : C******* The neighbors of NGF 244 : C 245 : MLENTI=MLEPOI.LECT(IFAC) 246 : NVOI=MLENTI.LECT(/1) 247 : C 248 : C******* The matrix where we store the coefficients 249 : C 250 : MATCOE=MLECOE.LECT(IFAC) 251 : C 252 : C******* Loop involving the neighbors 253 : C We create the matrix to "pseudoinvert" 254 : C 255 : DO IVOI=1,NVOI,1 256 : NGV=MLENTI.LECT(IVOI) 257 : IPCOOR=((NGV-1)*(IDIM+1))+1 258 : XV=MCOORD.XCOOR(IPCOOR) 259 : YV=MCOORD.XCOOR(IPCOOR+1) 260 : MATA.MAT(IVOI,1)=1 261 : MATA.MAT(IVOI,2)=XV-XF 262 : MATA.MAT(IVOI,3)=YV-YF 263 : IF(IDIM.EQ.3)THEN 264 : ZV=MCOORD.XCOOR(IPCOOR+2) 265 : MATA.MAT(IVOI,4)=ZV-ZF 266 : ENDIF 267 : ENDDO 268 : C 269 : C********** Now we have to invert this matrix 270 : C 271 : LOGSVD=.TRUE. 271 : CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT, 272 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD, 273 : & MLRTRA.PROG) 274 : IF(IERSVD.NE.0)THEN 275 : C 276 : C************* SVD decomposition of the matrix does not work 277 : C 278 : LOGSVD=.FALSE. 279 : ELSE 280 : C 281 : C************ We check the condition number of MATA 282 : C 283 : SMAX=0.0D0 284 : DO I1=1,NVMIN,1 285 : SMAX=MAX(SMAX,MLRSIG.PROG(I1)) 286 : ENDDO 287 : SMIN=SMAX 288 : DO I1=1,NVMIN,1 289 : SMIN=MIN(SMIN,MLRSIG.PROG(I1)) 290 : ENDDO 291 : IF((SMIN/SMAX) .LT. ERRTOL)THEN 292 : LOGSVD=.FALSE. 293 : ENDIF 294 : ENDIF 295 : C 296 : CC TEST 297 : C write(*,*) 'LOGSVD=.FALSE.' 298 : C LOGSVD=.FALSE. 299 : C 300 : IF(LOGSVD)THEN 301 : C 302 : C********** INVSVD worked 303 : C 304 : C MATA = MATU MATSIG t(MATV) 305 : C inv(MATA) = MATV inv(MATSIG) t(MATU) 306 : C 307 : DO I4=1,NVMIN,1 308 : DO IVOI=1,NVOI,1 309 : DO I2=1,IDIM+1,1 310 : C I2=1 is the only coefficient we are not interested 311 : C in. But we computed it to verify that 312 : C sum_ivoi MATCOE.MAT(ivoi,1) = 1 313 : C 314 : MATCOE.MAT(I2,IVOI)=MATCOE.MAT(I2,IVOI)+ 315 : & (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4) 316 : & /MLRSIG.PROG(I4)) 317 : ENDDO 318 : ENDDO 319 : ENDDO 320 : ELSE 321 : WRITE (IOIMP,*) 'rlencf.eso' 322 : C 22 0 323 : C Opération malvenue. Résultat douteux 323 : CALL ERREUR(22) 324 : C 325 : C********** INVSVD does not worked 326 : C For each neighbor k we have to compute the solution 327 : C of 328 : C 329 : C t(MATA) MATA x = t(MATA) * b 330 : C 331 : C where b= \sum_l e_l \delta(k,l) = e_k 332 : C 333 : C To do that, we compute 334 : C 335 : C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 336 : C 337 : C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 338 : C [t(MATA) MATA] X_0] 339 : C 340 : C 341 : C********** We compute [t(MATA) MATA] 342 : C We store it in the upper triangle of MTAA(NVMIN,NVMIN) 343 : C 344 : DO I1=1,NVMIN,1 345 : DO I2=I1,NVMIN,1 346 : MTAA.MAT(I1,I2)=0.0D0 347 : ENDDO 348 : ENDDO 349 : C 350 : DO I1=1,NVMIN,1 351 : DO I2=I1,NVMIN,1 352 : DO I3=1,NVOI,1 353 : MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+ 354 : & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2)) 355 : ENDDO 356 : ENDDO 357 : ENDDO 358 : C 359 : C********** We compute [inve(t(MATA) MATA)] 360 : C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN) 361 : C 361 : CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT, 362 : & MLRIN2.PROG,IERR0) 363 : IF(IERR0.NE.0)THEN 364 : WRITE(IOIMP,*) 'subroutine rlencf.eso.' 365 : C 26 2 366 : C Tache impossible. Probablement données erronées 366 : CALL ERREUR(26) 367 : GOTO 9999 368 : ENDIF 369 : C 370 : C********** We complete MTAA and MINTAA 371 : C 372 : DO I1=1,NVMIN,1 373 : DO I2=I1+1,NVMIN,1 374 : MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2) 375 : MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2) 376 : ENDDO 377 : ENDDO 378 : C 379 : DO IVOI=1,NVOI,1 380 : C 381 : C************* We compute [t(MATA) . b] and we store it in MLRIN1.PROG 382 : C 383 : DO I1=1,NVMIN,1 384 : MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1) 385 : MLRIN2.PROG(I1)=0.0D0 386 : MLRIN3.PROG(I1)=0.0D0 387 : ENDDO 388 : C 389 : C************* X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 390 : C X_0(i1) into MLRIN2.PROG(I1) 391 : C 392 : DO I2=1,NVMIN,1 393 : DO I1=1,NVMIN,1 394 : MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+ 395 : & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2)) 396 : ENDDO 397 : ENDDO 398 : C 399 : C************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 400 : C [t(MATA) MATA] X_0] 401 : C 402 : C [t(MATA) MATA] X_0 into MLRIN3.PROG 403 : C 404 : DO I2=1,NVMIN,1 405 : DO I1=1,NVMIN,1 406 : MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+ 407 : & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2)) 408 : ENDDO 409 : ENDDO 410 : C 411 : C************* Now we have 412 : C [t(MATA) . b] in MLRIN1.PROG 413 : C X_0(i1) in MLRIN2.PROG(I1) 414 : C [t(MATA) MATA] X_0 in MLRIN3.PROG 415 : C 416 : C X_1(i1) in MATCOE.MAT(IVOI,i1) 417 : C 418 : DO I1=1,IDIM+1,1 419 : C The only unuseful one is I1=1 420 : DO I2=1,NVMIN,1 421 : MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+ 422 : & (MINTAA.MAT(I1,I2)* 423 : & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2))) 424 : ENDDO 425 : MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+ 426 : & MLRIN2.PROG(I1) 427 : ENDDO 428 : ENDDO 429 : ENDIF 430 : CC 431 : CC TEST 432 : C write(*,*) 'ngf =', NGF 433 : C write(*,*) 'invide',LOGSVD 434 : C write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1) 435 : C write(*,*) 'coeff(1) =',(matcoe.mat(1,ivoi),ivoi=1,nvoi,1) 436 : C write(*,*) 'coeff(2) =',(matcoe.mat(2,ivoi),ivoi=1,nvoi,1) 437 : C write(*,*) 'coeff(3) =',(matcoe.mat(3,ivoi),ivoi=1,nvoi,1) 438 : C if(idim.eq.3) write(*,*) 'coeff(4)=', 439 : C & (matcoe.mat(4,ivoi),ivoi=1,nvoi,1) 440 : C xv=0.0D0 441 : C do ivoi=1,nvoi,1 442 : C xv=xv+matcoe.mat(1,ivoi) 443 : C enddo 444 : C write(*,*) 'sum=',xv 445 : CC 446 : MLENTI=MLEPOI.LECT(IFAC) 447 : SEGDES MLENTI 448 : MATCOE=MLECOE.LECT(IFAC) 449 : SEGDES MATCOE 450 : ENDDO 451 : C 452 : SEGDES MATCOE 453 : SEGDES MLEPOI 454 : C 455 : SEGDES MELFL 456 : SEGSUP MLEFP 457 : C 458 : SEGSUP MATA 459 : SEGSUP MATU 460 : SEGSUP MATV 461 : SEGSUP MLRSIG 462 : SEGSUP MLRTRA 463 : C 464 : SEGSUP MTAA 465 : SEGSUP MINTAA 466 : SEGSUP MLRIN1 467 : SEGSUP MLRIN2 468 : SEGSUP MLRIN3 469 : C 470 : 9999 CONTINUE 471 : RETURN 472 : END 473 : 474 : 475 :
© Cast3M 2003 - All rights reserved.
Disclaimer