1 : C ALEA1 SOURCE CHAT 11/03/16 21:15:20 6902 2 : SUBROUTINE ALEA1 3 : & (MMODEL,IPOSI,LADIM,XDIR1,XDIR2,XDIR3,ZSIG,CLAMD1,CLAMD2,CLAMD3, 4 : & VALMOY,DELZET,OMMAX) 5 : C 6 : C======================================================================= 7 : C 8 : C Subroutine ALEA1 : génération d'un CHAMELEM aléatoire 9 : C 10 : C Appellée par ALEA 11 : C 12 : C--------------------- 13 : C Variables internes : 14 : C--------------------- 15 : C 16 : C NBEL : nombre d'éléments 17 : C NBPT : nombre de noeuds par éléments 18 : C NBPGAU : nombre de valeurs à stocker par élément 19 : C 20 : C IPTABL : pointeur sur la table domaine, si elle existe 21 : C 22 : C-------------------------------------- 23 : C Développements ultérieurs possibles : 24 : C-------------------------------------- 25 : C 26 : C Génération aux points de gauss pour la MECA-FLUX (pression et vitesse) 27 : C 28 : C======================================================================= 29 : C 30 : IMPLICIT INTEGER(I-N) 31 : IMPLICIT REAL*8(A-H,O-Z) 32 : REAL*8 ZDIST, VALMOY, OMMAX 33 : C 34 : -INC CCOPTIO 35 : -INC CCREEL 36 : -INC CCGEOME 37 : -INC SMELEME 38 : -INC SMCOORD 39 : -INC SMCHAML 40 : -INC SMMODEL 41 : -INC SMTABLE 42 : -INC SMINTE 43 : 44 : C Table des coordonnées des points supports 45 : SEGMENT TABCOR 46 : REAL*8 COR(NBL,3) 47 : ENDSEGMENT 48 : 49 : C Table des valeurs en ces points 50 : SEGMENT TABVAL 51 : REAL*8 VAL(NBL) 52 : ENDSEGMENT 53 : 54 : C Table des partitions ok pour la génération 55 : SEGMENT TPART 56 : LOGICAL LPART(NBPART) 57 : ENDSEGMENT 58 : 59 : C Objet de type 'LISTENTI' utilise par ELQUOI. 60 : SEGMENT/INFO/(INFELL(JG)),INFO1.INFO,INFO2.INFO 61 : 62 : C Coordonnées reduites et fonctions de forme pour trouver 63 : C les coordonnées réelles des points de gauss 64 : SEGMENT WRK2 65 : REAL*8 XE(3,NBBB),SHP(6,NBBB) 66 : ENDSEGMENT 67 : 68 : DIMENSION XDIR1(3),XDIR2(3),XDIR3(3) 69 : CHARACTER*16 NOMFOR 70 : POINTEUR MELCHP.MELEME, MELENT.MELEME, MELSUP.MELEME 71 : 72 : INTEGER IPOSI 73 : C 74 : C epsilon servant à différents tests 75 : EPS = 1.D-12 76 : * 77 : * ------------------------------------------------------------------ 78 : * 2. On construit les coordonnées des points supports 79 : * on initialise le champ résultat en parallèle (CHAMELEM) 80 : * ------------------------------------------------------------------ 81 : 82 : SEGACT MMODEL 83 : * Nombre de partitions : 84 : NBPART = KMODEL(/1) 85 : 86 : * En-tête du CHAMELEM 87 : * Le nombre de partitions est fixé à zéro pour l'instant, et 88 : * il va augmenter à chaque partition invalide du modèle. 89 : L1 = 40 90 : N1 = 0 91 : N3 = 6 92 : SEGINI MCHELM 93 : TITCHE=' CHAMELEM cree par l operateur ALEA ' 94 : IFCHE = IFOUR 95 : * 96 : * Boucle sur les partitions du modèle. 97 : * Le chamelem résultat aura autant de partitions qu'il y a de 98 : * sous-maillages supports correspondant. 99 : * 100 : NBL = 1 101 : SEGINI TABCOR 102 : SEGINI TPART 103 : 104 : NBL = 0 105 : DO 4 IPART=1,NBPART 106 : IMODEL = KMODEL(IPART) 107 : SEGACT IMODEL 108 : NOMFOR = FORMOD(1) 109 : 110 : * a priori, cette partition convient 111 : LPART(IPART) = .TRUE. 112 : * 113 : * ------------------------------------------------------------------ 114 : * Vérification de la correspondance modèle - support demandé 115 : * 116 : IF ((IPOSI.EQ.2).OR.(IPOSI.EQ.6)) THEN 117 : * On vérifie que le modèle est bien un modèle de méca-flux 118 : * (NAVIER-STOKES ou DARCY), qui ,seuls, permettent l'usage de la 119 : * table domaine, que l'on récupère au passage. 120 : * Sinon, on passe à la partition suivante. 121 : IF ((NOMFOR.NE.'NAVIER_STOKES').AND.(NOMFOR.NE.'DARCY') 122 : & .AND.(NOMFOR.NE.'EULER')) THEN 123 : LPART(IPART) = .FALSE. 124 : GOTO 4 125 : ELSE 125 : CALL LEKMOD(MMODEL,IPTABL,IRET) 126 : ENDIF 127 : ENDIF 128 : 129 : IF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN 130 : * On vérifie que le modèle est bien un modèle MECANIQUE 131 : * sinon, on passe à la partition suivante. 132 : IF ((NOMFOR.EQ.'NAVIER_STOKES').OR.(NOMFOR.EQ.'DARCY') 133 : & .AND.(NOMFOR.NE.'EULER')) THEN 134 : LPART(IPART) = .FALSE. 135 : GOTO 4 136 : ENDIF 137 : ENDIF 138 : 139 : IF ((IPOSI.EQ.7).OR.(IPOSI.EQ.8)) THEN 140 : * On vérifie que le modèle est bien un modèle NAVIER-STOKES 141 : * ou EULER sinon, on passe à la partition suivante. 142 : IF ((NOMFOR.NE.'NAVIER_STOKES').AND.(NOMFOR.NE.'EULER')) THEN 143 : LPART(IPART) = .FALSE. 144 : GOTO 4 145 : ENDIF 146 : ENDIF 147 : * 148 : * N° d'élément fini NELE 149 : * nombre d'éléments NBEL 150 : * nombre de noeuds par éléments NBPT 151 : * nombre de fonctions de forme NBNN (pts de gauss seulement) 152 : NELE = NEFMOD 153 : 154 : MELENT = IMAMOD 155 : SEGACT MELENT 156 : NBPT = MELENT.NUM(/1) 157 : NBEL = MELENT.NUM(/2) 158 : 159 : * ------------------------------------------------------------------ 160 : * 2a. Détermination du support des points de valeur MELSUP s'ils existent, 161 : * du nombre de valeurs à stocker par élément NBPGAU, 162 : * et des pointeurs sur les segments d'intégration MINTE 163 : 164 : IF (IPOSI.EQ.1) THEN 165 : * Noeuds 166 : 167 : * On utilise tous les noeuds du maillage fourni 168 : * En particulier, en MECA-FLUX, tous les noeuds des QUAF 169 : MELSUP = MELENT 170 : SEGACT MELSUP 171 : NBPGAU = MELSUP.NUM(/1) 172 : SEGDES MELSUP 173 : MINTE = 0 174 : 175 : ELSEIF(IPOSI.EQ.2) THEN 176 : * Points centres 177 : 178 : * on les tire du maillage QUAF par la table domaine 178 : CALL LEKTAB(IPTABL,'CENTRE',MELSUP) 179 : IF (IERR.NE.0) THEN 180 : * On ne trouve pas le support qui contient les points 180 : CALL ERREUR(248) 181 : RETURN 182 : ENDIF 183 : if(infmod(/1).lt.iposi+2) then 183 : CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL) 184 : MINTE=infell(11) 185 : SEGSUP INFO 186 : else 187 : MINTE = infmod(4) 188 : endif 189 : NBPGAU = 1 190 : 191 : ELSEIF(IPOSI.EQ.6) THEN 192 : * Points faces 193 : 194 : * !!!!!!!!!! La structure de Castem ne permet pas encore de 195 : * créer un CHAMELEM appuyé aux faces. 196 : * Option noeud indisponible 196 : CALL ERREUR(507) 197 : RETURN 198 : 199 : ** on les tire du maillage QUAF par la table domaine 199 : * CALL LEKTAB(IPTABL,'FACE',MELSUP) 200 : * IF (IERR.NE.0) THEN 201 : ** On ne trouve pas le support qui contient les points 201 : * CALL ERREUR(248) 202 : * RETURN 203 : * ENDIF 203 : * CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL) 204 : * MINTE = INFELE(11) 205 : * SEGSUP INFO 206 : * SEGACT MELSUP 207 : * NBPGAU = MELSUP.NUM(/1) 208 : * SEGDES MELSUP 209 : 210 : ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN 211 : * Points de gauss MECANIQUE 212 : 212 : * CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL) 213 : IF(IPOSI.EQ.3) then 214 : NBPGAU = INFELE(6) 215 : minte=infmod(5) 216 : ELSEIF(IPOSI.EQ.4) then 217 : NBPGAU = INFELE(3) 218 : minte=infmod(6) 219 : ELSEIF(IPOSI.EQ.5) then 220 : NBPGAU = INFELE(4) 221 : MINTE = INFMOD(7) 222 : ENDIF 223 : NBNN = INFELE(8) 224 : * SEGSUP INFO 225 : 226 : ELSEIF ((IPOSI.EQ.7).OR.(IPOSI.EQ.8)) THEN 227 : * !!!!!!! Implémentation possible ci-après pour 228 : * les points de gauss du modèle NAVIER-STOKES (VITESSE et PRESSION) 229 : 230 : * Option noeud indisponible (pour l'instant) 230 : CALL ERREUR(507) 231 : RETURN 232 : 233 : ELSE 234 : * Erreur anormale.contactez votre support 234 : CALL ERREUR(5) 235 : ENDIF 236 : SEGDES MELENT 237 : * 238 : * ------------------------------------------------------------------ 239 : * 2b. Contenu par partition du CHAMELEM : 240 : * 241 : * on ajoute une partition au CHAMELEM 242 : N1 = N1 + 1 243 : SEGADJ MCHELM 244 : 245 : * Le maillage pointé par le CHAMELEM est le maillage du modèle. 246 : 247 : IF (NOMFOR.EQ.'DARCY') THEN 248 : * !!!!!!!! les tests des calculs DARCY exigent encore le 249 : * maillage de sommets, récupéré par la table domaine, alors 250 : * que le modèle pointe sur le maillage de QUAF. 250 : CALL LEKMOD(MMODEL,IPTABL,IRET) 251 : IF (IERR.NE.0) THEN 252 : C Données incompatibles 252 : CALL ERREUR(21) 253 : RETURN 254 : ENDIF 255 : * la table domaine existe. Le CHAMELEM va pointer sur le maillage 256 : * de sommets. On suppose bien entendu que les partitions du 257 : * modèle sont parcourues dans le même ordre que les 258 : * sous-maillages du maillage. 258 : CALL LEKTAB(IPTABL,'MAILLAGE',IPT1) 259 : IF (IERR.NE.0) RETURN 260 : SEGACT IPT1 261 : IF (IPT1.LISOUS(/1).EQ.0) THEN 262 : MELCHM = IPT1 263 : ELSE 264 : MELCHM = IPT1.LISOUS(IPART) 265 : ENDIF 266 : SEGDES IPT1 267 : ELSE 268 : * Le CHAMELEM pointe sur le maillage du modèle 269 : * (NAVIER-STOKES ou MECANIQUE). 270 : MELCHM = MELENT 271 : ENDIF 272 : 273 : CONCHE(N1) = CONMOD 274 : IMACHE(N1) = MELCHM 275 : INFCHE(N1,1) = 0 276 : INFCHE(N1,2) = 0 277 : INFCHE(N1,3) = NIFOUR 278 : INFCHE(N1,4) = MINTE 279 : INFCHE(N1,5) = 0 280 : INFCHE(N1,6) = IPOSI 281 : * 282 : * ------------------------------------------------------------------ 283 : * 2c. Construction des coordonnées des points supports : 284 : * dans le repère réel. 285 : * 286 : * combien de points de valeur à calculer en tout ? 287 : L = NBL 288 : NBL = NBL + (NBEL * NBPGAU) 289 : SEGADJ TABCOR 290 : 291 : * Cas où les points existent, et leurs coordonnées sont dans XCOOR. 292 : IF ((IPOSI.EQ.1).OR.(IPOSI.EQ.2).OR.(IPOSI.EQ.6)) THEN 293 : 294 : * Ce sont les points sommets, centres ou faces. 295 : SEGACT MELSUP 296 : DO 9 I=1,NBEL 297 : DO 10 J=1,NBPGAU 298 : L = L + 1 299 : IREF = (IDIM+1) * (MELSUP.NUM(J,I)-1) 300 : DO 11 K=1,IDIM 301 : COR(L,K) = XCOOR(IREF+K) 302 : 11 CONTINUE 303 : 10 CONTINUE 304 : 9 CONTINUE 305 : SEGDES MELSUP 306 : 307 : ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN 308 : 309 : * Les points n'existent pas. Ce sont les points de gauss MECANIQUE. 310 : * On calcule leurs coordonnées à partir des fonctions de forme 311 : * et des coordonnées des points sommets associés. 312 : SEGACT MELENT 313 : SEGACT MINTE 314 : NBBB = NBEL 315 : SEGINI WRK2 316 : 317 : DO 5 IEL=1,NBEL 318 : 319 : * Coordonnées XE des points supports des fonctions de forme 320 : * de l'élément IEL 321 : * Le maillage MELENT contient nécessairement au moins ces points 321 : CALL DOXE(XCOOR,IDIM,NBNN,MELENT.NUM,IEL,XE) 322 : * 323 : * Boucle sur les points de gauss 324 : DO 1100 IGAU=1,NBPGAU 325 : * 326 : L = L + 1 327 : 328 : XX=0.D0 329 : YY=0.D0 330 : ZZ=0.D0 331 : DO 1101 ISH=1,NBNN 332 : * valeur de la fonction de forme n°ISH au point de gauss IGAU 333 : CGAUSS = SHPTOT(1,ISH,IGAU) 334 : * coordonnée réelle du point support des fonctions de forme 335 : * coordonnées réelles du point de gauss 336 : XX = XX + XE(1,ISH)*CGAUSS 337 : YY = YY + XE(2,ISH)*CGAUSS 338 : IF (IDIM.EQ.3) THEN 339 : ZZ = ZZ + XE(3,ISH)*CGAUSS 340 : ENDIF 341 : 1101 CONTINUE 342 : 343 : COR(L,1) = XX 344 : COR(L,2) = YY 345 : COR(L,3) = ZZ 346 : 347 : 1100 CONTINUE 348 : 349 : 5 CONTINUE 350 : SEGDES MELENT,MINTE,WRK2 351 : 352 : ELSE 353 : 354 : * Erreur anormale.contactez votre support 354 : CALL ERREUR(5) 355 : RETURN 356 : ENDIF 357 : * 358 : * ------------------------------------------------------------------ 359 : * 2d. Transformation des coordonnées dans le repère adimensionné 360 : * par la matrice lambda de dimension LADIM. 361 : * 362 : IF (IDIM.EQ.2) THEN 363 : IF (LADIM.EQ.1) THEN 364 : * 2D stat 1D 365 : DO 20 L=1,NBL 366 : XX = COR(L,1) 367 : YY = COR(L,2) 368 : COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1 369 : 20 CONTINUE 370 : ELSE 371 : * 2D stat 2D 372 : DO 21 L=1,NBL 373 : XX = COR(L,1) 374 : YY = COR(L,2) 375 : COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1 376 : COR(L,2) = ((XX * XDIR2(1)) + (YY * XDIR2(2))) / CLAMD2 377 : 21 CONTINUE 378 : ENDIF 379 : ELSE 380 : IF (LADIM.EQ.1) THEN 381 : * 3D stat 1D 382 : DO 22 L=1,NBL 383 : XX = COR(L,1) 384 : YY = COR(L,2) 385 : ZZ = COR(L,3) 386 : COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2)) 387 : & + (ZZ * XDIR1(3)) ) / CLAMD1 388 : 22 CONTINUE 389 : ELSEIF (LADIM.EQ.2) THEN 390 : * 3D stat 2D 391 : DO 23 L=1,NBL 392 : XX = COR(L,1) 393 : YY = COR(L,2) 394 : ZZ = COR(L,3) 395 : COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2)) 396 : & + (ZZ * XDIR1(3)) ) / CLAMD1 397 : COR(L,2) = ( (XX * XDIR2(1)) + (YY * XDIR2(2)) 398 : & + (ZZ * XDIR2(3)) ) / CLAMD2 399 : 23 CONTINUE 400 : ELSE 401 : * 3D stat 3D 402 : DO 24 L=1,NBL 403 : XX = COR(L,1) 404 : YY = COR(L,2) 405 : ZZ = COR(L,3) 406 : COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2)) 407 : & + (ZZ * XDIR1(3)) ) / CLAMD1 408 : COR(L,2) = ( (XX * XDIR2(1)) + (YY * XDIR2(2)) 409 : & + (ZZ * XDIR2(3)) ) / CLAMD2 410 : COR(L,3) = ( (XX * XDIR3(1)) + (YY * XDIR3(2)) 411 : & + (ZZ * XDIR3(3)) ) / CLAMD3 412 : 24 CONTINUE 413 : ENDIF 414 : ENDIF 415 : 416 : * 417 : * ------------------------------------------------------------------ 418 : * 2e. Contenu par composante (1 seule, 'SCAL') du CHAMELEM : 419 : * 420 : N2 = 1 421 : SEGINI MCHAML 422 : ICHAML(N1) = MCHAML 423 : NOMCHE(1) = 'SCAL ' 424 : TYPCHE(1) = 'REAL*8 ' 425 : N1PTEL = NBPGAU 426 : N1EL = NBEL 427 : N2PTEL = 0 428 : N2EL = 0 429 : SEGINI MELVAL 430 : IELVAL(1) = MELVAL 431 : 432 : SEGDES IMODEL 433 : 434 : 4 CONTINUE 435 : SEGDES MMODEL 436 : 437 : IF (NBL.EQ.0) THEN 438 : * on n'a pas trouvé de partition avec le modèle correspondant 439 : * aux points de gauss demandés : 440 : C Données incompatibles 440 : CALL ERREUR(21) 441 : RETURN 442 : ENDIF 443 : * 444 : * ------------------------------------------------------------------ 445 : * 3. Génération aléatoire : 446 : * ------------------------------------------------------------------ 447 : * 448 : * Si les points supports appartiennent à différentes mailles en 449 : * même temps, il faut savoir que les valeurs en des points 450 : * confondus sont égales avec la méthode des bandes tournantes. 451 : 451 : CALL BANTOU(TABCOR,TABVAL,LADIM,ZSIG,VALMOY,DELZET,OMMAX) 452 : SEGDES TABCOR 453 : * 454 : * ------------------------------------------------------------------ 455 : * 4. Construction champ résultat : 456 : * ------------------------------------------------------------------ 457 : * 458 : * BANTOU renvoie le résultat sous la forme d'une table de valeur, 459 : * TABVAL active, que l'on doit maintenant convertir en MELVAL, 460 : * partition par partition (on parcourt les partitions valides). 461 : 462 : L = 1 463 : ICHELM = 0 464 : DO 8 IPART=1,NBPART 465 : IF (LPART(IPART)) THEN 466 : ICHELM = ICHELM + 1 467 : MCHAML = ICHAML(ICHELM) 468 : MELVAL = IELVAL(1) 469 : 470 : NBEL = VELCHE(/2) 471 : NBPGAU = VELCHE(/1) 472 : DO 1 I=1,NBEL 473 : DO 19 J=1,NBPGAU 474 : VELCHE(J,I) = VAL(L) 475 : L = L + 1 476 : 19 CONTINUE 477 : 1 CONTINUE 478 : ENDIF 479 : 8 CONTINUE 480 : 481 : * ------------------------------------------------------------------ 482 : * 5. Ecriture et sortie 483 : * ------------------------------------------------------------------ 484 : * 485 : * Fermeture des segments 486 : * 487 : DO 12 IPART=1,N1 488 : MCHAML = ICHAML(IPART) 489 : MELVAL = IELVAL(1) 490 : SEGDES MELVAL,MCHAML 491 : 12 CONTINUE 492 : SEGDES MCHELM 493 : SEGSUP TABVAL,TPART 494 : * 495 : * Résultat 496 : * 496 : CALL ECROBJ('MCHAML ',MCHELM) 497 : 498 : RETURN 499 : END 500 : 501 : 502 : 503 : 504 : 505 : 506 : 507 : 508 : 509 : 510 : 511 : 512 : 513 :
© Cast3M 2003 - All rights reserved.
Disclaimer