1 : C GRQUA2 SOURCE CHAT 05/01/13 00:21:26 5004 2 : SUBROUTINE GRQUA2(MCENVO,MCENEL,IGEOLI,INFOMC) 3 : C------------------------------------------------------------------ 4 : C PROJET : CASTEM 2000 5 : C 6 : C NOM : GRQUA2 7 : C 8 : C DESCRIPTION : Calcule les coefficients à appliquer 9 : C sur la fonction dont onv eut calculer le gradient et 10 : C le Hessien. 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 (Outils) : ERREUR, KRIPAD 19 : C 20 : C 21 : C APPELES (Calcul) : GRQUA3, GRQUA4 22 : C 23 : C----------------------------------------------------------------- 24 : C 25 : C ENTREES : 26 : C 27 : C MCENVO : maillage centre-voisins 28 : C 29 : C MCENEL : maillage de centre correspondant 30 : C à MCENVO 31 : C 32 : C IGEOLI : SPG de CHPOINT de C.L. 33 : C 34 : C SORTIES : 35 : C INFOMC : 36 : C segment qui contient les infos 37 : C pour créer le MCHAML des coefficients 38 : C 39 : C----------------------------------------------------------------- 40 : C 41 : C HISTORIQUE (Anomalies et modifications eventuelles) 42 : C 43 : C HISTORIQUE : 23.04.99, Creation 44 : C 45 : C----------------------------------------------------------------- 46 : C 47 : C 48 : IMPLICIT INTEGER(I-N) 49 : -INC CCOPTIO 50 : -INC SMCOORD 51 : -INC SMELEME 52 : -INC SMLENTI 53 : -INC SMLREEL 54 : C 55 : C----------------------------------------------------------------- 56 : C 57 : SEGMENT LISTPO 58 : INTEGER LISPO1(NPOILI),INDEX1(NBCEN+1) 59 : ENDSEGMENT 60 : C 61 : C**** LISTPO: list de points centres avec leur voisins 62 : C et position: 63 : C 64 : C LISPO1(k) t.q. 65 : C LISPO1(INDEX1(j)) = numero global du centre j 66 : C LISPO1(INDEX1(j)),...,LISPO1(INDEX1(j+1))-1 = voisins du centre 67 : C LISPO1(INDEX1(j)) 68 : C 69 : C**** Des tableuax a deux indices 70 : C 71 : SEGMENT MATRIX 72 : REAL*8 MAT(N1,N2) 73 : ENDSEGMENT 74 : C 75 : POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX,MATCEL.MATRIX, 76 : & COEF.MATRIX 77 : C 78 : C**** Segment INFOMC 79 : C 80 : SEGMENT INFOMC 81 : INTEGER NTOTSU,LAST(NLISOU,2), REST(NBCEN) 82 : INTEGER LISPOI(NPOILI),INDEX(NBCEN+1) 83 : REAL*8 COEFF(NPOILI,LDIM) 84 : ENDSEGMENT 85 : C 86 : C**** INFOMC segment qui contient les informations pour créer le 87 : C MCHAML pour le calcul du gradient et de l'hessien 88 : C 89 : C NTOTSU = nombre total de maillages elementaires 90 : C 91 : C LAST(i,2) = nombre de elts avec i-noueds 92 : C LAST(i,1) = si LAST(i,2) != 0, numero local du premier noeuds 93 : C avec i-1 elts; les autres sont dans REST(LAST(i,1)), 94 : C REST(REST(LAST(i,1))) si != 0 95 : C si j <= NBCEN, INDEX(j) contient la position du centre j dedans 96 : C LISPOI(k), i.e. 97 : C LISPOI(INDEX(j)) = numero global du centre j 98 : C LISPOI(INDEX(j)),...,LISPOI(INDEX(j+1))-1 = voisins du centre 99 : C LISPOI(INDEX(j)) 100 : C 101 : C N.B. Les informations ici sont redondantes, i.e. 102 : C si j appartient à la liste chaine qui parte de LAST(i,1), 103 : C INDEX(j+1)-INDEX(j) = i 104 : C 105 : C Mais ca nous permet de definir le MCHAML avec beaucoup de 106 : C facilité! 107 : C 108 : C LDIM = 5 en 2D, 9 en 3D 109 : C COEFF(j,*) = (d()/dx,d()/dy,d2()/dxdx,d2()/dydy,d2()/dxdy)_j en 2D 110 : C COEFF(j,*) = (d()/dx,d()/dy,d()/dz,d2()/dxdx,d2()/dydy,d2()/dzdz, 111 : C d2()/dxdy,d2()/dxdz,d2()/dydz)_j 112 : C 113 : POINTEUR MCENVO.MELEME,MCENEL.MELEME,KRILIM.MLENTI,KRICEN.MLENTI 114 : POINTEUR MLEVOI.MLENTI,MLEVO2.MLENTI,MLEVTO.MLENTI,MLRDIS.MLREEL, 115 : & MLRSIG.MLREEL, MLRCEL.MLREEL 116 : 117 : C 118 : C**** Variables de COOPTIO 119 : C 120 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 121 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 122 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 123 : C & ,IECHO, IIMPI, IOSPI 124 : C & ,IDIM 125 : CC & ,MCOORD 126 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 127 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 128 : C & ,NORINC,NORVAL,NORIND,NORVAD 129 : C & ,NUCROU, IPSAUV 130 : C 131 : INTEGER IGEOLI, NTSOUS, NTOTP, ICEN, ILISOU,NBELEM, NBNOEU, NTCEN2 132 : & , NLISOU, LDIM, NBMAX, NPOILI, NBCEN, IPOS, IVOIS 133 : & , IELEM, NOCENE, NOCEN1, IPOI 134 : & , NVOIMI 135 : & , IPOSNE, JG, INDER, IVECT, N1, N2, NTOTP0 136 : C 137 : REAL*8 CONTOL 138 : C 139 : C**** D'abord, on doit dimensioner le segment INFOMC 140 : C Boucle sur le LISOUS de MCENVO 141 : C 142 : SEGACT MCENVO 143 : SEGACT MCENEL 144 : C 145 : NTSOUS = MCENVO.LISOUS(/1) 146 : C 147 : C**** MCENEL est un maillage de POI1 148 : C 149 : IF((MCENEL.LISOUS(/1) .NE. 0) .OR. (MCENEL.ITYPEL .NE. 1))THEN 150 : C 151 : C******* 'MCENEL must contain only POI1 elements' 152 : C 153 : MOTERR(1:40) = 'Subroutine grqua2.eso. 1 ' 154 : WRITE(IOIMP,*) MOTERR(1:40) 155 : GOTO 9991 156 : ENDIF 157 : C 158 : C**** On calcule le nombre total du point dedans MCENVO 159 : C (les points doubles sont contés deux fois!) 160 : C 161 : NTOTP = 0 162 : NBMAX = 0 163 : ICEN = 0 164 : IF(NTSOUS .EQ. 0) NTSOUS = 1 165 : C 166 : DO ILISOU = 1, NTSOUS 167 : IF(NTSOUS .EQ. 1)THEN 168 : IPT1 = MCENVO 169 : ELSE 170 : IPT1 = MCENVO.LISOUS(ILISOU) 171 : SEGACT IPT1 172 : ENDIF 173 : NBNOEU = IPT1.NUM(/1) 174 : NBMAX = MAX(NBMAX, NBNOEU) 175 : C 176 : C******* NBMAX = nombre max de noeuds dans un elt 177 : C Son valeur va dimensioner INFOMC (NLISOU) 178 : C 179 : NBELEM = IPT1.NUM(/2) 180 : NTOTP = NTOTP + (NBNOEU * NBELEM) 181 : ICEN = ICEN + NBELEM 182 : ENDDO 183 : C 184 : NTCEN2 = MCENEL.NUM(/2) 185 : IF(NTCEN2 .NE. ICEN)THEN 186 : MOTERR(1:40) = 'Subroutine grqua2.eso. 2 ' 187 : WRITE(IOIMP,*) MOTERR(1:40) 188 : GOTO 9991 189 : ENDIF 190 : 191 : C 192 : NBCEN = NTCEN2 193 : NPOILI = NTOTP 194 : SEGINI LISTPO 195 : C 196 : C**** KRIPADs de MCENEL et de IGEOLI 197 : C 197 : IF(IGEOLI .GT. 0) CALL KRIPAD(IGEOLI,KRILIM) 197 : CALL KRIPAD(MCENEL,KRICEN) 198 : C 199 : C**** Dedans KRIPAD: 200 : C SEGACT KRILIM 201 : C SEGACT KRICEN 202 : C 203 : C 204 : C**** Boucle sur les LISOUS de MCENVO 205 : C 206 : ICEN = 0 207 : IPOS = 0 208 : C 209 : C**** IPOS: pour placer les noeuds dedans LISTPO.LISPO1 210 : C pour conter le nombre de points et dimensioner INFOMC 211 : C 212 : DO ILISOU = 1, NTSOUS 213 : IF(NTSOUS .EQ. 1)THEN 214 : IPT1 = MCENVO 215 : ELSE 216 : IPT1 = MCENVO.LISOUS(ILISOU) 217 : C 218 : C********** IPT1 déjà activés 219 : C 220 : ENDIF 221 : NBNOEU = IPT1.NUM(/1) 222 : NBELEM = IPT1.NUM(/2) 223 : DO IELEM = 1 , NBELEM 224 : NOCENE = IPT1.NUM(NBNOEU,IELEM) 225 : ICEN = ICEN + 1 226 : NOCEN1 = MCENEL.NUM(1,ICEN) 227 : IF(NOCENE .NE. NOCEN1)THEN 228 : MOTERR(1:40) = 'Subroutine grqua2.eso. 3 ' 229 : WRITE(IOIMP,*) MOTERR(1:40) 230 : GOTO 9991 231 : ENDIF 232 : IVOIS = 1 233 : IPOS = IPOS + 1 234 : LISTPO.INDEX1(ICEN)=IPOS 235 : LISTPO.LISPO1(IPOS) = NOCENE 236 : DO IPOI = 1, NBNOEU - 1 237 : NOCEN1 = IPT1.NUM(IPOI,IELEM) 238 : IF(KRICEN.LECT(NOCEN1) .NE. 0)THEN 239 : IPOS = IPOS + 1 240 : IVOIS = IVOIS + 1 241 : LISTPO.LISPO1(IPOS) = NOCEN1 242 : ELSEIF(IGEOLI .GT. 0)THEN 243 : IF(KRILIM.LECT(NOCEN1) .NE. 0)THEN 244 : IPOS = IPOS + 1 245 : IVOIS = IVOIS + 1 246 : LISTPO.LISPO1(IPOS) = NOCEN1 247 : ENDIF 248 : ENDIF 249 : ENDDO 250 : ENDDO 251 : IF(NTSOUS .GT. 1) SEGDES IPT1 252 : ENDDO 253 : LISTPO.INDEX1(NBCEN+1)=IPOS+1 254 : SEGDES MCENVO 255 : SEGDES MCENEL 256 : IF(IGEOLI .GT.0) SEGSUP KRILIM 257 : C 258 : C**** N.B: LISTPO est utilisé pour créer "les voisins des voisins"; 259 : C donc il faut le completer avant de créer cette liste. 260 : C 261 : C**** Maintenant on va s'amuser! 262 : C 263 : C 1) On doit faire une boucle sur les noeuds de 264 : C MCENEL 265 : C 266 : C 2) en 2D il faut avoir au moins 5 voisins (6 noeuds au total) 267 : C en 3D il faut avoir au moins 9 noeuds (NVOIMI) 268 : C 269 : IF(IDIM .EQ. 2)THEN 270 : NVOIMI = 6 271 : ELSE 272 : NVOIMI = 10 273 : ENDIF 274 : C 275 : C**** NTOTVO < NVOIMI: 276 : C 277 : C Il faut élargir le stencil 278 : C 279 : C Je prends les noeuds les plus proches 280 : C 281 : C En suit, il faut assembler A et la "pseudo-inverser" 282 : C Si rank(A) < NVOIMI, alors il faut repeter l'operation 283 : C d'élargissement du stencil. 284 : C 285 : C S'elle n'est pas inversible: sortir et obliger l'utilisateur a 286 : C refaire le maillage 287 : C 288 : C**** Si non (si NTOTVO >= NVOIMI) 289 : C 290 : C il faut assembler A et la "pseudo-inverser" . 291 : C 292 : C Si rank(A) < NVOIMI, alors il faut élargir le stencil ... 293 : C 294 : C**** Tous ca, on le stoke en INFOMC 295 : C 296 : C NBCEN déjà defini 297 : C 298 : NPOILI = IPOS 299 : NLISOU = NBMAX 300 : LDIM = NVOIMI - 1 301 : SEGINI INFOMC 302 : C 303 : C**** Object temporaires: 304 : C 305 : C MLEVOI.MLENTI : pour consider les voisins une seule fois! 306 : C liste chainé du noeud et ses voisins 307 : C 308 : C On dois aussi dimensioner: 309 : C 310 : C MLEVO2.MLENTI: liste de points "voisins des voisins" du 311 : C centre du quel on veut le gradient 312 : C 313 : C NBMAX-1 voisins qui ont NBMAX-1 voisins - le 314 : C point centre meme (NBMAX -1) * (NBMAX - 2) 315 : C 316 : C MLRDIS.MLREEL: distances des points MLEVO2 du point centre 317 : C 318 : C MLEVTO.MLENTI: le point centre et tous les voisins utilize 319 : C pour le calcul du gradient (NBMAX x (NBMAX-1))! 320 : C 321 : C Matrices pour la "pseudo-inversion" de A avec la "SVD" 322 : C (Singular Values Decomposition) 323 : C 324 : C MATA.MATRIX : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1) 325 : C 326 : C MATU.MATRIX : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1) 327 : C 328 : C MATV.MATRIX : matrice (NVOIMI-1) x (NVOIMI-1) 329 : C 330 : C MLRSIG.MLREEL : (NVOIMI-1) valeur singuliers de la matrice 331 : C MATA 332 : C 333 : C MLRCEL.MLREEL : segment de travail de longuer 334 : C (NBMAX x (NBMAX - 1)) 335 : C 336 : C MATCEL.MATRIX : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1) 337 : C qui contient la transpose de la 338 : C pseudo-inverse de MATA 339 : C 340 : C COEF.MATRIX : matrice (NBMAX x (NBMAX - 1)) x (NVOIMI-1) 341 : C qui contient la transpose de la 342 : C pseudo-inverse de MATA 343 : C 344 : C N.B.: MATA = MATU . diag(MLRSIG) . transpose(MATV) 345 : C 346 : C pseudoinverse(MATA) = 347 : C MATV . inv(diag(MLRSIG)) . transpose(MATU) 348 : C 349 : JG=MCOORD.XCOOR(/1)/(IDIM+1) 350 : SEGINI MLEVOI 351 : C 352 : JG = (NBMAX - 1) * (NBMAX - 2) 353 : SEGINI MLEVO2 354 : SEGINI MLRDIS 355 : C 356 : JG = (NBMAX - 1) * NBMAX 357 : SEGINI MLEVTO 358 : C 359 : N1 = (NBMAX - 1) * NBMAX 360 : N2 = NVOIMI - 1 361 : SEGINI MATA 362 : SEGINI MATCEL 363 : C 364 : N1 = (NBMAX - 1) * NBMAX 365 : N2 = NVOIMI - 1 366 : SEGINI MATU 367 : C 368 : N1 = (NBMAX - 1) * NBMAX 369 : N2 = NVOIMI - 1 370 : SEGINI MATV 371 : C 372 : JG = NVOIMI - 1 373 : SEGINI MLRSIG 374 : SEGINI MLRCEL 375 : C 376 : N1 = (NBMAX - 1) * NBMAX 377 : N2 = NVOIMI-1 378 : SEGINI COEF 379 : C 380 : IPOSNE = 0 381 : C 382 : C IPOSNE = position des noueuds dedans INFOMC.LISPOI 383 : C 384 : DO IELEM = 1, NBCEN 385 : C 386 : C******* D'abord, il faut remplir MLEVTO 387 : C definir NTOTP 388 : C 389 : IPOS = LISTPO.INDEX1(IELEM) 390 : NOCENE = LISTPO.LISPO1(IPOS) 391 : MLEVTO.LECT(1)=NOCENE 392 : NTOTP = LISTPO.INDEX1(IELEM+1) - IPOS 393 : DO IVOIS = 1, NTOTP - 1, 1 394 : NOCEN1 = LISTPO.LISPO1(IPOS+IVOIS) 395 : MLEVTO.LECT(1+IVOIS)=NOCEN1 396 : ENDDO 397 : IF(NTOTP .LT. NVOIMI)THEN 398 : INDER = 1 399 : ELSE 400 : CONTOL = 1.0D-4 400 : CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, 401 : & MATCEL,MLRSIG,MLRCEL,COEF,INDER) 402 : ENDIF 403 : IF(INDER .EQ. 2)THEN 404 : C 405 : C************* Erreur anormale 406 : C 407 : MOTERR(1:40) = 'Subroutine grqua2.eso. 4 ' 408 : WRITE(IOIMP,*) MOTERR(1:40) 409 : GOTO 9991 410 : ELSEIF(INDER .EQ. 1)THEN 411 : C 412 : NTOTP0 = NTOTP 413 : CONTOL = 1.0D-12 413 : CALL GRQUA4(LISTPO,MLEVOI,MLEVO2,MLRDIS, 414 : & KRICEN,CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO, 415 : & MATA,MATU,MATV,MATCEL,MLRSIG,MLRCEL,COEF,INDER) 416 : C 417 : IF(INDER .EQ. 1)THEN 418 : GOTO 9992 419 : ELSEIF(INDER .EQ. 2)THEN 420 : GOTO 9991 421 : ENDIF 422 : C 423 : IF(NTOTP .GT. (NBMAX - 1) * NBMAX)THEN 424 : C 425 : C************* Erreur anormale 426 : C 427 : MOTERR(1:40) = 'Subroutine grqua2.eso. 5 ' 428 : WRITE(IOIMP,*) MOTERR(1:40) 429 : GOTO 9991 430 : ENDIF 431 : C 432 : C************* N.B.: Il faut ajuster INFOMC 433 : C 434 : NLISOU = MAX(NLISOU, NTOTP) 435 : NPOILI = NPOILI + NTOTP - NTOTP0 436 : SEGADJ INFOMC 437 : ENDIF 438 : IPOSNE = IPOSNE + 1 439 : INFOMC.INDEX(IELEM)=IPOSNE 440 : INFOMC.LISPOI(IPOSNE)=NOCENE 441 : DO IVECT = 1, NVOIMI - 1 , 1 442 : INFOMC.COEFF(IPOSNE,IVECT) = COEF.MAT(1,IVECT) 443 : ENDDO 444 : DO IVOIS= 2, NTOTP,1 445 : IPOSNE = IPOSNE + 1 446 : NOCEN1 = MLEVTO.LECT(IVOIS) 447 : INFOMC.LISPOI(IPOSNE)=NOCEN1 448 : C 449 : C********** Cette boucle est au contraire. 450 : C Voir s'on peut faire mieux! 451 : C 452 : DO IVECT = 1, NVOIMI - 1 , 1 453 : INFOMC.COEFF(IPOSNE,IVECT) = COEF.MAT(IVOIS,IVECT) 454 : ENDDO 455 : ENDDO 456 : C 457 : IF(INFOMC.LAST(NTOTP,2) .EQ. 0) 458 : $ INFOMC.NTOTSU = INFOMC.NTOTSU + 1 459 : INFOMC.LAST(NTOTP,2) = INFOMC.LAST(NTOTP,2) + 1 460 : INFOMC.REST(IELEM) = INFOMC.LAST(NTOTP,1) 461 : C 462 : C******* Si INFOMC.LAST(NTOTP,2)=0, alors INFOMC.LAST(NTOTP,1)=0 463 : C 464 : INFOMC.LAST(NTOTP,1) = IELEM 465 : ENDDO 466 : INFOMC.INDEX(NBCEN+1)=IPOSNE+1 467 : SEGSUP LISTPO 468 : SEGDES INFOMC 469 : SEGSUP KRICEN 470 : SEGSUP MLEVOI 471 : SEGSUP MLEVO2 472 : SEGSUP MLRDIS 473 : SEGSUP MLEVTO 474 : SEGSUP MATA 475 : SEGSUP MATCEL 476 : SEGSUP MATU 477 : SEGSUP MATV 478 : SEGSUP MLRSIG 479 : SEGSUP MLRCEL 480 : SEGSUP COEF 481 : C 482 : 9999 RETURN 483 : C 484 : C**** Erreurs 485 : C 486 : 9991 CONTINUE 487 : C 488 : C**** Message d'erreur standard 489 : C 5 3 490 : C Erreur anormale.contactez votre support 491 : C 491 : CALL ERREUR(5) 492 : RETURN 493 : C 494 : C**** Donnés incompatibles 495 : C 496 : 9992 CONTINUE 497 : MOTERR(1:40) = 'Le maillage réalisé ne permet pas une ' 498 : WRITE(IOIMP,*) MOTERR(1:40) 499 : MOTERR(1:40) = 'reconstruction quadratique exacte. ' 500 : WRITE(IOIMP,*) MOTERR(1:40) 501 : C 502 : C**** Message d'erreur standard 503 : C 21 2 504 : C Données incompatibles 505 : C 505 : CALL ERREUR(21) 506 : RETURN 507 : C 508 : C**** Option non implementée 509 : C 510 : 9993 CONTINUE 511 : C 512 : C**** Message d'erreur standard 513 : C 251 2 514 : C Tentative d'utilisation d'une option non implémentée 515 : C 515 : CALL ERREUR(251) 516 : RETURN 517 : C 518 : END 519 : 520 : 521 :
© Cast3M 2003 - All rights reserved.
Disclaimer