Download alea.eso

Back to the list

   1 : C ALEA      SOURCE    CHAT      11/03/16    21:15:24     6902           
   2 :       SUBROUTINE ALEA
   3 : C
   4 : C=======================================================================        
   5 : C                                                                               
   6 : C  Opérateur ALEA : génération d'un champ scalaire aléatoire gaussien 
   7 : C  stationnaire isotrope par la méthode des bandes tournantes.
   8 : C
   9 : C---------------------------
  10 : C Phrase d'appel (GIBIANE) :
  11 : C---------------------------
  12 : C
  13 : C     OBJ1 = 'ALEA' 'BANDES_TOURNANTES' | MODE1  (MOT1)
  14 : C                                       | MAIL1
  15 : C            'EXPO' 'SIGMA' FLOT1 ('MOYENNE' FLOT2) 
  16 : C                  |  'LAMBDA'  FLOT3 
  17 : C                  |  'LAMBDA1' FLOT4 (VEC4) 
  18 : C                    ('LAMBDA2' FLOT5 (VEC5))
  19 : C                    ('LAMBDA3' FLOT6 (VEC6))
  20 : C     
  21 : C------------------------
  22 : C Opérandes et résultat :
  23 : C------------------------
  24 : C
  25 : C
  26 : C 'BANDES_TOURNANTES' 
  27 : C             : mot-clé indiquant que l'on utilise la méthode des bandes
  28 : C               tournantes
  29 : C
  30 : C       MODE1 : Modèle sur lequel s'appuie le champ résultat
  31 : C               (type MMODEL), pour obtenir un champ par élément en sortie.
  32 : C
  33 : C        MOT1 : mot-clef facultatif valant 'NOEUD','GRAVITE','RIGIDITE',
  34 : C               'MASSE','STRESSES','FACE','VITESSE' ou 'PRESSION'
  35 : C               indiquant quels points supports prendre en compte dans la
  36 : C               génération du champ (défaut = 'NOEUD'). 
  37 : C
  38 : C       MAIL1 : Maillage sur lequel s'appuie le champ résultat
  39 : C               (type MELEME), pour obtenir un champ par point en sortie.
  40 : C
  41 : C      'EXPO' : mot-clé indiquant que la loi de covariance est
  42 : C               exponentielle.
  43 : C
  44 : C     'SIGMA' : mot-clé suivi de :
  45 : C
  46 : C       FLOT1 : écart-type du champ à générer (type FLOTTANT).
  47 : C
  48 : C   'MOYENNE' : mot-clé optionnel suivi de :
  49 : C
  50 : C       FLOT2 : valeur de la moyenne du champ aléatoire (type FLOTTANT)
  51 : C               par défaut = 0.     
  52 : C
  53 : C    'LAMBDA' : mot-clé (cas de la corrélation isotrope) suivi de :
  54 : C
  55 : C       FLOT3 : longueur de corrélation isotrope (type FLOTTANT).
  56 : C
  57 : C   'LAMBDA1' : mots-clés (cas de la corrélation anisotrope) 
  58 : C  ('LAMBDA2')  autant que la dimension de la structure de corrélation, 
  59 : C  ('LAMBDA3')  suivis respectivement de :
  60 : C
  61 : C       FLOT4 : longueurs de corrélation (type FLOTTANT) dans les 3 
  62 : C      (FLOT5)    
  63 : C      (FLOT6)  directions principales.
  64 : C
  65 : C       VEC4  : directions optionnelles orthogonales des axes principaux 
  66 : C      (VEC5)   de corrélation 1, 2, et 3 respectivement (type POINT). 
  67 : C      (VEC6)   Ils doivent être non nuls et orthogonaux. 
  68 : C               Par défaut, ce sont les axes x, y et z.
  69 : C
  70 : C        OBJ1 : champ résultat (type 'CHPOINT ' ou 'MCHAML' selon que
  71 : C               l'on donne un maillage ou un modèle en entrée), 
  72 : C               diffus (si c'est un champ-point), une composante 'SCAL'.
  73 : C
  74 : C------------
  75 : C Remarques :
  76 : C------------
  77 : C
  78 : C           1.  Les options 'RIGIDITE','MASSE','STRESSES' réclament un
  79 : C               modèle de mécanique.
  80 : C
  81 : C           2.  Les options 'VITESSE','PRESSION' réclament un modèle
  82 : C               NAVIER-STOKES (non implémenté pour l'instant).
  83 : C
  84 : C           3.  Les options 'GRAVITE','FACE' réclament un modèle
  85 : C               NAVIER-STOKES ou DARCY (usage de la table domaine)
  86 : C
  87 : C           4.  L'option FACE n'est pas encore mise en service.
  88 : C
  89 : C--------------------------------------
  90 : C Développements ultérieurs possibles :
  91 : C--------------------------------------
  92 : C
  93 : C Génération aux points de gauss pour la MECA-FLUX (pression et vitesse)
  94 : C
  95 : C=======================================================================        
  96 : C
  97 :       IMPLICIT INTEGER(I-N)
  98 :       IMPLICIT REAL*8(A-H,O-Z)                                                  
  99 :       REAL*8 ZDIST, VALMOY, OMMAX
 100 : C
 101 : -INC CCOPTIO
 102 : -INC CCREEL
 103 : -INC CCGEOME
 104 : -INC SMCOORD
 105 : -INC SMELEME
 106 : -INC SMMODEL
 107 : 
 108 :       DIMENSION    XDIR1(3),XDIR2(3),XDIR3(3)
 109 :       CHARACTER*20 MOTMET(2)
 110 :       CHARACTER*9  MOTD(1)
 111 :       CHARACTER*8  MOTL(1), MOTS(1), MOTLA(4), MOTM(1)
 112 :       CHARACTER*8  MOTOBJ(2), MOTSUP(8)
 113 : 
 114 :       INTEGER      IOBJ, IPOSI
 115 : 
 116 :       LOGICAL      AXES
 117 : C
 118 :       DATA MOTMET /'BANDES_TOURNANTES','MATRICIELLE'/
 119 :       DATA MOTL   /'EXPO    '/
 120 :       DATA MOTS   /'SIGMA   '/
 121 :       DATA MOTLA  /'LAMBDA  ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/
 122 :       DATA MOTM   /'MOYENNE '/
 123 :       DATA MOTD   /'DIRECTION'/
 124 :       DATA MOTSUP /'NOEUD','GRAVITE','RIGIDITE','MASSE','STRESSES',
 125 :      &             'FACE','VITESSE','PRESSION'/
 126 : 
 127 : C     epsilon servant à différents tests (1.D-10)
 128 :       EPS = 1.D5 * (SQRT(XPETIT))
 129 : c
 130 : c     ================  Lecture données  ================
 131 : c
 132 : C
 133 : C     Lecture du mode de génération
 134 : C     -----------------------------
 134 :       CALL LIRMOT(MOTMET,2,IMET,1)
 135 :       IF (IERR.NE.0) RETURN 
 136 : *
 137 : *     Lecture du modèle ou du maillage de points
 138 : *     La présence du modèle indique que l'on désire un objet CHAMELEM
 139 : *     en sortie. Sinon, un maillage est fourni et on rendra un CHAMPOINT
 140 : *
 140 :       CALL LIROBJ('MMODEL',MMODEL,0,IRET)
 141 :       IF (IERR.NE.0) RETURN
 142 :       IF (IRET.NE.0) THEN
 143 : *       Il y a un modèle
 144 :         IOBJ = 1
 145 :       ELSE
 146 : *       Il y a un maillage
 147 :         IOBJ = 2
 147 :         CALL LIROBJ('MAILLAGE',MELENT,1,IRET)
 148 :         IF (IERR.NE.0) RETURN
 149 :       ENDIF
 150 : C
 151 : C     Lecture du support (points centre ou sommet) - facultatif
 152 : C     ------------------
 152 :       CALL LIRMOT(MOTSUP,8,IPOSI,0)
 153 :       IF (IERR.NE.0) RETURN 
 154 :       IF (IPOSI.EQ.0) THEN 
 155 : C       valeur par défaut = NOEUD
 156 :         IPOSI = 1 
 157 :       ENDIF
 158 : C
 159 : C     Lecture du mot-clé 'EXPO'
 160 : C     -------------------------
 160 :       CALL LIRMOT(MOTL,1,IRET,1)
 161 :       IF (IERR.NE.0) RETURN 
 162 :       IF (IRET.EQ.0) THEN 
 163 :         MOTERR(1:4) = 'EXPO'
 164 : c       Il manque le mot-clé %m1:4
 164 :         CALL ERREUR(396)
 165 :         RETURN
 166 :       ENDIF
 167 : C
 168 : C     Lecture du mot-clé 'SIGMA'
 169 : C     --------------------------
 169 :       CALL LIRMOT(MOTS,1,IRET,1)
 170 :       IF (IERR.NE.0) RETURN 
 171 :       IF (IRET.EQ.0) THEN 
 172 :         MOTERR(1:4) = 'SIGM'
 173 : c       Il manque le mot-clé %m1:4
 173 :         CALL ERREUR(396)
 174 :         RETURN
 175 :       ENDIF
 176 : C     Lecture de la valeur de sigma (strictement supérieure à 0.D0)
 176 :       CALL LIRREE(ZSIG,1,IRET)
 177 :       IF (IERR.NE.0) RETURN 
 178 :       IF (IRET.EQ.0) THEN 
 179 :         MOTERR(1:4) = 'SIGM'
 180 : c       Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 180 :         CALL ERREUR(166)
 181 :         RETURN
 182 :       ENDIF
 183 :       IF (ZSIG.LE.0.D0) THEN
 184 :         REAERR(1) = REAL(ZSIG)
 185 :         REAERR(2) = REAL(0.D0)
 186 :         MOTERR(1:8) = ' SIGMA '
 187 : c       %m1:8 = %r1 inférieur à %r2
 187 :         CALL ERREUR(41)
 188 :         RETURN
 189 :       ENDIF
 190 : C
 191 : C     Lecture du mot-clé 'MOYENNE'
 192 : C     ----------------------------
 192 :       CALL LIRMOT(MOTM,1,IRET,0)
 193 :       IF (IERR.NE.0) RETURN 
 194 :       IF (IRET.EQ.0) THEN 
 195 : *       valeur par défaut = 0
 196 :         VALMOY = 0.D0
 197 :       ELSE
 198 : C       Lecture de la valeur de la moyenne
 198 :         CALL LIRREE(VALMOY,1,IRET)
 199 :         IF (IERR.NE.0) RETURN 
 200 :         IF (IRET.EQ.0) THEN 
 201 :           MOTERR(1:4) = 'MOYE'
 202 : c         Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 202 :           CALL ERREUR(166)
 203 :           RETURN
 204 :         ENDIF
 205 :       ENDIF
 206 : C
 207 : C     Lecture du ou des mot(s)-clé(s) 'LAMBDA' (isotrope) ou 'LAMBDAx' (anisotrope)
 208 : C     Initialisation des CLAMDx et XDIRx.
 209 : C     -----------------------------------------------------------------------------
 209 :       CALL LIRMOT(MOTLA,4,ILAM,1)
 210 :       IF (IERR.NE.0) RETURN 
 211 :       IF (ILAM.EQ.0) THEN 
 212 :         MOTERR(1:4) = 'LAMB'
 213 : c       Il manque le mot-clé %m1:4
 213 :         CALL ERREUR(396)
 214 :         RETURN
 215 :       ENDIF
 216 : 
 217 :       IF (ILAM.EQ.1) THEN 
 218 : C       Cas ISOTROPE
 219 : 
 220 : *       la statistique a la même dimension que l'espace réel :
 221 :         LADIM = IDIM
 222 : 
 223 : C       Lecture de la valeur de lambda (strictement supérieure à 0.D0)
 223 :         CALL LIRREE(CLAMDA,1,IRET)
 224 :         IF (IERR.NE.0) RETURN 
 225 :         IF (IRET.EQ.0) THEN 
 226 :           MOTERR(1:4) = 'LAMB'
 227 : c         Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 227 :           CALL ERREUR(166)
 228 :           RETURN
 229 :         ENDIF
 230 :         IF (CLAMDA.LE.0.D0) THEN
 231 :           REAERR(1) = REAL(CLAMDA)
 232 :           REAERR(2) = REAL(0.D0)
 233 :           MOTERR(1:8) = 'LAMBDA  '
 234 : c         %m1:8 = %r1 inférieur à %r2
 234 :           CALL ERREUR(41)
 235 :           RETURN
 236 :         ENDIF
 237 : 
 238 : *       mêmes longueurs de corrélation dans toutes les directions
 239 :         CLAMD1 = CLAMDA
 240 :         CLAMD2 = CLAMDA
 241 :         CLAMD3 = CLAMDA
 242 : *       les vecteurs principaux sont les axes x, y et z
 243 :         XDIR1(1) = 1.D0
 244 :         XDIR1(2) = 0.D0
 245 :         XDIR1(3) = 0.D0
 246 :         XDIR2(1) = 0.D0
 247 :         XDIR2(2) = 1.D0
 248 :         XDIR2(3) = 0.D0
 249 :         XDIR3(1) = 0.D0
 250 :         XDIR3(2) = 0.D0
 251 :         XDIR3(3) = 1.D0
 252 :         
 253 :       ELSE
 254 : C       Cas ANISOTROPE
 255 :    
 256 : C       axe 1 : 
 257 : C       on doit avoir lu le mot-clef 'LAMBDA1'
 258 :         IF (ILAM.NE.2) THEN 
 259 : *         Syntaxe incorrecte : on attend %m1:30
 260 :           MOTERR(1:30) = 'le mot-clef LAMBDA1           '
 260 :           CALL ERREUR(881)
 261 :           RETURN
 262 :         ENDIF
 263 : C       Lecture de la valeur de lambda1 (strictement supérieure à 0.D0)
 263 :         CALL LIRREE(CLAMD1,1,IRET)
 264 :         IF (IERR.NE.0) RETURN 
 265 :         IF (IRET.EQ.0) THEN 
 266 :           MOTERR(1:4) = 'LAMB'
 267 : c         Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 267 :           CALL ERREUR(166)
 268 :           RETURN
 269 :         ENDIF
 270 :         IF (CLAMD1.LE.0.D0) THEN
 271 :           REAERR(1) = REAL(CLAMD1)
 272 :           REAERR(2) = REAL(0.D0)
 273 :           MOTERR(1:8) = 'LAMBDA1 '
 274 : c         %m1:8 = %r1 inférieur à %r2
 274 :           CALL ERREUR(41)
 275 :           RETURN
 276 :         ENDIF
 277 : 
 278 : C       Lecture optionnelle de la direction VEC1
 278 :         CALL LIROBJ('POINT',IPTV,0,IRET)
 279 :         IF (IERR.NE.0) RETURN
 280 : 
 281 : *       si on ne donne pas d'axe, ce sont les axes par défaut, et la
 282 : *       statistique a la même dimensinonalité que l'espace réel.
 283 :         AXES = (IRET.NE.0)
 284 : 
 285 :         IF (.NOT.AXES) THEN 
 286 : *         pas de vecteurs => axes par défaut.
 287 :           XDIR1(1) = 1.D0
 288 :           XDIR1(2) = 0.D0
 289 :           XDIR1(3) = 0.D0
 290 :           XDIR2(1) = 0.D0
 291 :           XDIR2(2) = 1.D0
 292 :           XDIR2(3) = 0.D0
 293 :           XDIR3(1) = 0.D0
 294 :           XDIR3(2) = 0.D0
 295 :           XDIR3(3) = 1.D0
 296 : *         la statistique a la même dimension que l'espace réel :
 297 :           LADIM = IDIM
 298 :         ELSE
 299 : C         on charge le vecteur et on vérifie qu'il n'est pas de longueur nulle
 300 : C         on le normalise aussi.
 301 :           IREF = (IPTV-1) * (IDIM+1)
 302 :           SDIR = 0.D0
 303 :           DO 13 I=1,IDIM
 304 :             XDIR1(I) = XCOOR(IREF+I)
 305 :             SDIR = SDIR + (XDIR1(I) * XDIR1(I))
 306 :  13       CONTINUE
 307 :           SDIR     = SQRT(SDIR)
 308 : C         sinon, erreur :
 309 :           IF (SDIR.LT.EPS) THEN
 310 : c           Une direction ne peut pas être définie par un vecteur nul
 310 :             CALL ERREUR(239)
 311 :             RETURN
 312 :           ENDIF
 313 : *         normalisation
 314 :           XDIR1(1) = XDIR1(1) / SDIR
 315 :           XDIR1(2) = XDIR1(2) / SDIR
 316 :           XDIR1(3) = XDIR1(3) / SDIR
 317 : *         la statistique est au moins de dimension 1 
 318 : *         (à surcharger si nécessaire) :
 319 :           LADIM = 1
 320 :         ENDIF
 321 : 
 322 : C       axe 2 : 
 323 : C       Si le mot-clef 'LAMBDA2' existe, stat 2D au moins
 323 :         CALL LIRMOT(MOTLA,4,ILAM,0)
 324 :         IF (IERR.NE.0) RETURN 
 325 : 
 326 :         IF ((ILAM.NE.0).AND.(ILAM.NE.3)) THEN 
 327 : *         Syntaxe incorrecte : on attend %m1:30
 328 :           MOTERR(1:30) = 'le mot-clef LAMBDA2           '
 328 :           CALL ERREUR(881)
 329 :           RETURN
 330 :         ENDIF
 331 :         IF (ILAM.EQ.3) THEN 
 332 : *         la statistique est au moins de dimension 2
 333 :           LADIM = 2
 334 : 
 335 : C         Lecture de la valeur de lambda2 (strictement supérieure à 0.D0)
 335 :           CALL LIRREE(CLAMD2,1,IRET)
 336 :           IF (IERR.NE.0) RETURN 
 337 :           IF (IRET.EQ.0) THEN 
 338 :             MOTERR(1:4) = 'LAMB'
 339 : c           Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 339 :             CALL ERREUR(166)
 340 :             RETURN
 341 :           ENDIF
 342 :           IF (CLAMD2.LE.0.D0) THEN
 343 :             REAERR(1) = REAL(CLAMD2)
 344 :             REAERR(2) = REAL(0.D0)
 345 :             MOTERR(1:8) = 'LAMBDA2 '
 346 : c           %m1:8 = %r1 inférieur à %r2
 346 :             CALL ERREUR(41)
 347 :             RETURN
 348 :           ENDIF
 349 : 
 350 : C         Lecture obligatoire de la direction VEC2 si on a donné VEC1
 351 :           IF (AXES) THEN
 351 :             CALL LIROBJ('POINT',IPTV,1,IRET)
 352 :             IF (IERR.NE.0) RETURN
 353 : 
 354 :             IF (IRET.EQ.0) THEN 
 355 : *             Syntaxe incorrecte : on attend %m1:30
 356 :               MOTERR(1:30) = 'un vecteur                    '
 356 :               CALL ERREUR(881)
 357 :               RETURN
 358 :             ENDIF
 359 : 
 360 : C           on charge le vecteur et on vérifie qu'il n'est pas de longueur
 361 : C           nulle, et qu'il est bien perpendiculaire à VEC1.
 362 : C           on le normalise aussi.
 363 :             IREF = (IPTV-1) * (IDIM+1)
 364 :             SDIR = 0.D0
 365 :             PSCL = 0.D0
 366 :             DO 14 I=1,IDIM
 367 :               XDIR2(I) = XCOOR(IREF+I)
 368 :               SDIR = SDIR + (XDIR2(I) * XDIR2(I))
 369 :               PSCL = PSCL + (XDIR1(I) * XDIR2(I))
 370 :  14         CONTINUE
 371 :             SDIR   = SQRT(SDIR)
 372 : *           si VEC2 de longueur nulle
 373 :             IF (SDIR.LT.EPS) THEN
 374 : c             Une direction ne peut pas être définie par un vecteur nul
 374 :               CALL ERREUR(239)
 375 :               RETURN
 376 :             ENDIF
 377 : *           si VEC2 non perpendiculaire à VEC1
 378 :             IF ((ABS(PSCL)).GE.(EPS*SDIR*SDIR)) THEN
 379 : c             Les vecteurs définissant le repère local ne sont pas orthogonaux
 379 :               CALL ERREUR(161)
 380 :               RETURN
 381 :             ENDIF
 382 : *           orthogonalisation à VEC1 :
 383 :             XDIR2(1) = XDIR2(1) - (PSCL * XDIR1(1))
 384 :             XDIR2(2) = XDIR2(2) - (PSCL * XDIR1(2))
 385 :             XDIR2(3) = XDIR2(3) - (PSCL * XDIR1(3))
 386 : *           normalisation
 387 :             XDIR2(1) = XDIR2(1) / SDIR
 388 :             XDIR2(2) = XDIR2(2) / SDIR
 389 :             XDIR2(3) = XDIR2(3) / SDIR
 390 :           ENDIF
 391 : 
 392 : C         axe 3 (si 3D) : 
 393 :           IF (IDIM.EQ.3) THEN
 394 : C           Si le mot-clef 'LAMBDA3' existe, stat 3D 
 394 :             CALL LIRMOT(MOTLA,4,ILAM,0)
 395 :             IF (IERR.NE.0) RETURN 
 396 : 
 397 :             IF ((ILAM.NE.0).AND.(ILAM.NE.4)) THEN 
 398 : *             Syntaxe incorrecte : on attend %m1:30
 399 :               MOTERR(1:30) = 'le mot-clef LAMBDA3           '
 399 :               CALL ERREUR(881)
 400 :               RETURN
 401 :             ENDIF
 402 :             IF (ILAM.EQ.4) THEN 
 403 : *             la statistique est de dimension 3
 404 :               LADIM = 3
 405 : 
 406 : C             Lecture de la valeur de lambda3 (strictement supérieure à 0.D0)
 406 :               CALL LIRREE(CLAMD3,1,IRET)
 407 :               IF (IERR.NE.0) RETURN 
 408 :               IF (IRET.EQ.0) THEN 
 409 :                 MOTERR(1:4) = 'LAMB'
 410 : c               Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
 410 :                 CALL ERREUR(166)
 411 :                 RETURN
 412 :               ENDIF
 413 :               IF (CLAMD3.LE.0.D0) THEN
 414 :                 REAERR(1) = REAL(CLAMD3)
 415 :                 REAERR(2) = REAL(0.D0)
 416 :                 MOTERR(1:8) = 'LAMBDA3 '
 417 : c               %m1:8 = %r1 inférieur à %r2
 417 :                 CALL ERREUR(41)
 418 :                 RETURN
 419 :               ENDIF
 420 : 
 421 : C             Lecture obligatoire de la direction VEC3 si on a donné 
 422 : C             VEC1 et VEC2
 423 :               IF (AXES) THEN
 423 :                 CALL LIROBJ('POINT',IPTV,1,IRET)
 424 :                 IF (IERR.NE.0) RETURN
 425 :   
 426 :                 IF (IRET.EQ.0) THEN 
 427 : *                 Syntaxe incorrecte : on attend %m1:30
 428 :                   MOTERR(1:30) = 'un vecteur                    '
 428 :                   CALL ERREUR(881)
 429 :                   RETURN
 430 :                 ENDIF
 431 : 
 432 : C               on charge le vecteur et on vérifie qu'il n'est pas de longueur
 433 : C               nulle, et qu'il est bien perpendiculaire à VEC1 et à VEC2.
 434 : C               on le normalise aussi.
 435 :                 IREF  = (IPTV-1) * (IDIM+1)
 436 :                 SDIR  = 0.D0
 437 :                 PSCL  = 0.D0
 438 :                 PSCL2 = 0.D0
 439 :                 DO 15 I=1,IDIM
 440 :                   XDIR3(I) = XCOOR(IREF+I)
 441 :                   SDIR  = SDIR + (XDIR3(I) * XDIR3(I))
 442 :                   PSCL  = PSCL + (XDIR1(I) * XDIR3(I))
 443 :                   PSCL2 = PSCL + (XDIR2(I) * XDIR3(I))
 444 :  15             CONTINUE
 445 :                 SDIR    = SQRT(SDIR)
 446 : *               si VEC3 de longueur nulle
 447 :                 IF (SDIR.LT.EPS) THEN
 448 : c                 Une direction ne peut pas être définie par un vecteur nul
 448 :                   CALL ERREUR(239)
 449 :                   RETURN
 450 :                 ENDIF
 451 : *               si VEC3 non perpendiculaire à VEC1 et VEC2
 452 :                 IF (((ABS(PSCL )).GE.(EPS*SDIR*SDIR)).OR.
 453 :      &              ((ABS(PSCL2)).GE.(EPS*SDIR*SDIR))) THEN
 454 : c                 Les vecteurs définissant le repère local ne sont pas 
 455 : *                 orthogonaux
 455 :                   CALL ERREUR(161)
 456 :                   RETURN
 457 :                 ENDIF
 458 : *               orthogonalisation à VEC1 :
 459 :                 XDIR2(1) = XDIR2(1) - (PSCL  * XDIR1(1)) 
 460 :      &                              - (PSCL2 * XDIR2(1))
 461 :                 XDIR2(2) = XDIR2(2) - (PSCL  * XDIR1(2)) 
 462 :      &                              - (PSCL2 * XDIR2(2))
 463 :                 XDIR2(3) = XDIR2(3) - (PSCL  * XDIR1(3)) 
 464 :      &                              - (PSCL2 * XDIR2(3))
 465 : *               normalisation
 466 :                 XDIR3(1) = XDIR3(1) / SDIR
 467 :                 XDIR3(2) = XDIR3(2) / SDIR
 468 :                 XDIR3(3) = XDIR3(3) / SDIR
 469 :               ENDIF                                                      lit vec3
 470 :             ENDIF                                                        stat 3D
 471 :           ENDIF                                                          dime 3D
 472 :         ENDIF                                                            stat 2D
 473 :       ENDIF                                                              anisotrope
 474 : c
 475 : c     ================  Travail  ================
 476 : c
 477 : C     Dans la suite, les coordonnées vont être adimensionnées par
 478 : C     Lambda. Donc toutes les données numériques seront calculées comme
 479 : C     si Lambda est isotrope égal à 1.
 480 : 
 481 : C     Longueur minimale de description des hétérogénéités
 482 : C     pour déterminer la plage de variation des vecteurs d'onde
 483 :       ZDIST = 0.2D0 
 484 : *     Pas de discrétisation de la bande (avec lambda = 1)
 485 :       IF (LADIM.EQ.1) THEN
 486 :         DELZET = ZDIST / 3.D0
 487 :       ELSE
 488 :         DELZET = ZDIST / 4.D0
 489 :       ENDIF
 490 : 
 491 : *     Fréquence de coupure (avec lambda = 1), 
 492 : *     elle doit être supérieure à 2 pi / dx
 493 :       OMMAX = 1.1D0 * (2.D0 * XPI / ZDIST)
 494 : 
 495 : *     Génération
 496 :       IF (IOBJ.EQ.1) THEN
 497 : *       On veut un CHAMELEM en sortie
 497 :         CALL ALEA1(MMODEL,IPOSI,LADIM,XDIR1,XDIR2,XDIR3,
 498 :      &             ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
 499 :       ELSEIF (IOBJ.EQ.2) THEN
 500 : *       On veut un CHAMPOINT en sortie
 500 :         CALL ALEA2(MELENT,LADIM,XDIR1,XDIR2,XDIR3,
 501 :      &             ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
 502 :       ENDIF
 503 : c
 504 : c     ================  Sortie  ================
 505 : c
 506 :       RETURN
 507 :       END
 508 : 
 509 : 
 510 : 
 511 :  
 512 :  
 513 :  
 514 :  
 515 :  
 516 :  
 517 :  
 518 :  
 519 :  

© Cast3M 2003 - All rights reserved.
Disclaimer