Télécharger alea.eso

Retour à la liste

Numérotation des lignes :

alea
  1. C ALEA SOURCE CB215821 19/07/31 21:15:20 10277
  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.  
  102. -INC PPARAM
  103. -INC CCOPTIO
  104. -INC CCREEL
  105. -INC CCGEOME
  106. -INC SMCOORD
  107. -INC SMELEME
  108. -INC SMMODEL
  109.  
  110. DIMENSION XDIR1(3),XDIR2(3),XDIR3(3)
  111. CHARACTER*20 MOTMET(2)
  112. CHARACTER*8 MOTL(1), MOTS(1), MOTLA(4), MOTM(1)
  113. CHARACTER*8 MOTOBJ(2), MOTSUP(8)
  114.  
  115. INTEGER IOBJ, IPOSI
  116.  
  117. LOGICAL AXES
  118. C
  119. DATA MOTMET /'BANDES_TOURNANTES','MATRICIELLE'/
  120. DATA MOTL /'EXPO '/
  121. DATA MOTS /'SIGMA '/
  122. DATA MOTLA /'LAMBDA ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/
  123. DATA MOTM /'MOYENNE '/
  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 -----------------------------
  135. CALL LIRMOT(MOTMET,2,IMET,1)
  136. IF (IERR.NE.0) RETURN
  137. *
  138. * Lecture du modèle ou du maillage de points
  139. * La présence du modèle indique que l'on désire un objet CHAMELEM
  140. * en sortie. Sinon, un maillage est fourni et on rendra un CHAMPOINT
  141. *
  142. CALL LIROBJ('MMODEL ',MMODEL,0,IRET)
  143. IF (IERR.NE.0) RETURN
  144. IF (IRET.NE.0) THEN
  145. * Il y a un modèle
  146. CALL ACTOBJ('MMODEL ',MMODEL,1)
  147. IOBJ = 1
  148. ELSE
  149. * Il y a un maillage
  150. IOBJ = 2
  151. CALL LIROBJ('MAILLAGE',MELENT,1,IRET)
  152. IF (IERR.NE.0) RETURN
  153. ENDIF
  154. C
  155. C Lecture du support (points centre ou sommet) - facultatif
  156. C ------------------
  157. CALL LIRMOT(MOTSUP,8,IPOSI,0)
  158. IF (IERR.NE.0) RETURN
  159. IF (IPOSI.EQ.0) THEN
  160. C valeur par défaut = NOEUD
  161. IPOSI = 1
  162. ENDIF
  163. C
  164. C Lecture du mot-clé 'EXPO'
  165. C -------------------------
  166. CALL LIRMOT(MOTL,1,IRET,1)
  167. IF (IERR.NE.0) RETURN
  168. IF (IRET.EQ.0) THEN
  169. MOTERR(1:4) = 'EXPO'
  170. c Il manque le mot-clé %m1:4
  171. CALL ERREUR(396)
  172. RETURN
  173. ENDIF
  174. C
  175. C Lecture du mot-clé 'SIGMA'
  176. C --------------------------
  177. CALL LIRMOT(MOTS,1,IRET,1)
  178. IF (IERR.NE.0) RETURN
  179. IF (IRET.EQ.0) THEN
  180. MOTERR(1:4) = 'SIGM'
  181. c Il manque le mot-clé %m1:4
  182. CALL ERREUR(396)
  183. RETURN
  184. ENDIF
  185. C Lecture de la valeur de sigma (strictement supérieure à 0.D0)
  186. CALL LIRREE(ZSIG,1,IRET)
  187. IF (IERR.NE.0) RETURN
  188. IF (IRET.EQ.0) THEN
  189. MOTERR(1:4) = 'SIGM'
  190. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  191. CALL ERREUR(166)
  192. RETURN
  193. ENDIF
  194. IF (ZSIG.LE.0.D0) THEN
  195. REAERR(1) = REAL(ZSIG)
  196. REAERR(2) = REAL(0.D0)
  197. MOTERR(1:8) = ' SIGMA '
  198. c %m1:8 = %r1 inférieur à %r2
  199. CALL ERREUR(41)
  200. RETURN
  201. ENDIF
  202. C
  203. C Lecture du mot-clé 'MOYENNE'
  204. C ----------------------------
  205. CALL LIRMOT(MOTM,1,IRET,0)
  206. IF (IERR.NE.0) RETURN
  207. IF (IRET.EQ.0) THEN
  208. * valeur par défaut = 0
  209. VALMOY = 0.D0
  210. ELSE
  211. C Lecture de la valeur de la moyenne
  212. CALL LIRREE(VALMOY,1,IRET)
  213. IF (IERR.NE.0) RETURN
  214. IF (IRET.EQ.0) THEN
  215. MOTERR(1:4) = 'MOYE'
  216. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  217. CALL ERREUR(166)
  218. RETURN
  219. ENDIF
  220. ENDIF
  221. C
  222. C Lecture du ou des mot(s)-clé(s) 'LAMBDA' (isotrope) ou 'LAMBDAx' (anisotrope)
  223. C Initialisation des CLAMDx et XDIRx.
  224. C -----------------------------------------------------------------------------
  225. CALL LIRMOT(MOTLA,4,ILAM,1)
  226. IF (IERR.NE.0) RETURN
  227. IF (ILAM.EQ.0) THEN
  228. MOTERR(1:4) = 'LAMB'
  229. c Il manque le mot-clé %m1:4
  230. CALL ERREUR(396)
  231. RETURN
  232. ENDIF
  233.  
  234. IF (ILAM.EQ.1) THEN
  235. C Cas ISOTROPE
  236.  
  237. * la statistique a la même dimension que l'espace réel :
  238. LADIM = IDIM
  239.  
  240. C Lecture de la valeur de lambda (strictement supérieure à 0.D0)
  241. CALL LIRREE(CLAMDA,1,IRET)
  242. IF (IERR.NE.0) RETURN
  243. IF (IRET.EQ.0) THEN
  244. MOTERR(1:4) = 'LAMB'
  245. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  246. CALL ERREUR(166)
  247. RETURN
  248. ENDIF
  249. IF (CLAMDA.LE.0.D0) THEN
  250. REAERR(1) = REAL(CLAMDA)
  251. REAERR(2) = REAL(0.D0)
  252. MOTERR(1:8) = 'LAMBDA '
  253. c %m1:8 = %r1 inférieur à %r2
  254. CALL ERREUR(41)
  255. RETURN
  256. ENDIF
  257.  
  258. * mêmes longueurs de corrélation dans toutes les directions
  259. CLAMD1 = CLAMDA
  260. CLAMD2 = CLAMDA
  261. CLAMD3 = CLAMDA
  262. * les vecteurs principaux sont les axes x, y et z
  263. XDIR1(1) = 1.D0
  264. XDIR1(2) = 0.D0
  265. XDIR1(3) = 0.D0
  266. XDIR2(1) = 0.D0
  267. XDIR2(2) = 1.D0
  268. XDIR2(3) = 0.D0
  269. XDIR3(1) = 0.D0
  270. XDIR3(2) = 0.D0
  271. XDIR3(3) = 1.D0
  272.  
  273. ELSE
  274. C Cas ANISOTROPE
  275.  
  276. C axe 1 :
  277. C on doit avoir lu le mot-clef 'LAMBDA1'
  278. IF (ILAM.NE.2) THEN
  279. * Syntaxe incorrecte : on attend %m1:30
  280. MOTERR(1:30) = 'le mot-clef LAMBDA1 '
  281. CALL ERREUR(881)
  282. RETURN
  283. ENDIF
  284. C Lecture de la valeur de lambda1 (strictement supérieure à 0.D0)
  285. CALL LIRREE(CLAMD1,1,IRET)
  286. IF (IERR.NE.0) RETURN
  287. IF (IRET.EQ.0) THEN
  288. MOTERR(1:4) = 'LAMB'
  289. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  290. CALL ERREUR(166)
  291. RETURN
  292. ENDIF
  293. IF (CLAMD1.LE.0.D0) THEN
  294. REAERR(1) = REAL(CLAMD1)
  295. REAERR(2) = REAL(0.D0)
  296. MOTERR(1:8) = 'LAMBDA1 '
  297. c %m1:8 = %r1 inférieur à %r2
  298. CALL ERREUR(41)
  299. RETURN
  300. ENDIF
  301.  
  302. C Lecture optionnelle de la direction VEC1
  303. CALL LIROBJ('POINT',IPTV,0,IRET)
  304. IF (IERR.NE.0) RETURN
  305.  
  306. * si on ne donne pas d'axe, ce sont les axes par défaut, et la
  307. * statistique a la même dimensinonalité que l'espace réel.
  308. AXES = (IRET.NE.0)
  309.  
  310. IF (.NOT.AXES) THEN
  311. * pas de vecteurs => axes par défaut.
  312. XDIR1(1) = 1.D0
  313. XDIR1(2) = 0.D0
  314. XDIR1(3) = 0.D0
  315. XDIR2(1) = 0.D0
  316. XDIR2(2) = 1.D0
  317. XDIR2(3) = 0.D0
  318. XDIR3(1) = 0.D0
  319. XDIR3(2) = 0.D0
  320. XDIR3(3) = 1.D0
  321. * la statistique a la même dimension que l'espace réel :
  322. LADIM = IDIM
  323. ELSE
  324. C on charge le vecteur et on vérifie qu'il n'est pas de longueur nulle
  325. C on le normalise aussi.
  326. IREF = (IPTV-1) * (IDIM+1)
  327. SDIR = 0.D0
  328. DO 13 I=1,IDIM
  329. XDIR1(I) = XCOOR(IREF+I)
  330. SDIR = SDIR + (XDIR1(I) * XDIR1(I))
  331. 13 CONTINUE
  332. SDIR = SQRT(SDIR)
  333. C sinon, erreur :
  334. IF (SDIR.LT.EPS) THEN
  335. c Une direction ne peut pas être définie par un vecteur nul
  336. CALL ERREUR(239)
  337. RETURN
  338. ENDIF
  339. * normalisation
  340. DO I=1,IDIM
  341. XDIR1(I) = XDIR1(I) / SDIR
  342. ENDDO
  343. * la statistique est au moins de dimension 1
  344. * (à surcharger si nécessaire) :
  345. LADIM = 1
  346. ENDIF
  347.  
  348. C axe 2 :
  349. C Si le mot-clef 'LAMBDA2' existe, stat 2D au moins
  350. CALL LIRMOT(MOTLA,4,ILAM,0)
  351. IF (IERR.NE.0) RETURN
  352.  
  353. IF ((ILAM.NE.0).AND.(ILAM.NE.3)) THEN
  354. * Syntaxe incorrecte : on attend %m1:30
  355. MOTERR(1:30) = 'le mot-clef LAMBDA2 '
  356. CALL ERREUR(881)
  357. RETURN
  358. ENDIF
  359. IF (ILAM.EQ.3) THEN
  360. * la statistique est au moins de dimension 2
  361. LADIM = 2
  362.  
  363. C Lecture de la valeur de lambda2 (strictement supérieure à 0.D0)
  364. CALL LIRREE(CLAMD2,1,IRET)
  365. IF (IERR.NE.0) RETURN
  366. IF (IRET.EQ.0) THEN
  367. MOTERR(1:4) = 'LAMB'
  368. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  369. CALL ERREUR(166)
  370. RETURN
  371. ENDIF
  372. IF (CLAMD2.LE.0.D0) THEN
  373. REAERR(1) = REAL(CLAMD2)
  374. REAERR(2) = REAL(0.D0)
  375. MOTERR(1:8) = 'LAMBDA2 '
  376. c %m1:8 = %r1 inférieur à %r2
  377. CALL ERREUR(41)
  378. RETURN
  379. ENDIF
  380.  
  381. C Lecture obligatoire de la direction VEC2 si on a donné VEC1
  382. IF (AXES) THEN
  383. CALL LIROBJ('POINT',IPTV,1,IRET)
  384. IF (IERR.NE.0) RETURN
  385.  
  386. IF (IRET.EQ.0) THEN
  387. * Syntaxe incorrecte : on attend %m1:30
  388. MOTERR(1:30) = 'un vecteur '
  389. CALL ERREUR(881)
  390. RETURN
  391. ENDIF
  392.  
  393. C on charge le vecteur et on vérifie qu'il n'est pas de longueur
  394. C nulle, et qu'il est bien perpendiculaire à VEC1.
  395. C on le normalise aussi.
  396. IREF = (IPTV-1) * (IDIM+1)
  397. SDIR = 0.D0
  398. PSCL = 0.D0
  399. DO 14 I=1,IDIM
  400. XDIR2(I) = XCOOR(IREF+I)
  401. SDIR = SDIR + (XDIR2(I) * XDIR2(I))
  402. PSCL = PSCL + (XDIR1(I) * XDIR2(I))
  403. 14 CONTINUE
  404. SDIR = SQRT(SDIR)
  405. * si VEC2 de longueur nulle
  406. IF (SDIR.LT.EPS) THEN
  407. c Une direction ne peut pas être définie par un vecteur nul
  408. CALL ERREUR(239)
  409. RETURN
  410. ENDIF
  411. * si VEC2 non perpendiculaire à VEC1
  412. IF ((ABS(PSCL)).GE.(EPS*SDIR*SDIR)) THEN
  413. c Les vecteurs définissant le repère local ne sont pas orthogonaux
  414. CALL ERREUR(161)
  415. RETURN
  416. ENDIF
  417. * orthogonalisation à VEC1 :
  418. DO I=1,IDIM
  419. XDIR2(I) = XDIR2(I) - (PSCL * XDIR1(I))
  420. ENDDO
  421. * normalisation
  422. DO I=1,IDIM
  423. XDIR2(I) = XDIR2(I) / SDIR
  424. ENDDO
  425. ENDIF
  426.  
  427. C axe 3 (si 3D) :
  428. IF (IDIM.EQ.3) THEN
  429. C Si le mot-clef 'LAMBDA3' existe, stat 3D
  430. CALL LIRMOT(MOTLA,4,ILAM,0)
  431. IF (IERR.NE.0) RETURN
  432.  
  433. IF ((ILAM.NE.0).AND.(ILAM.NE.4)) THEN
  434. * Syntaxe incorrecte : on attend %m1:30
  435. MOTERR(1:30) = 'le mot-clef LAMBDA3 '
  436. CALL ERREUR(881)
  437. RETURN
  438. ENDIF
  439. IF (ILAM.EQ.4) THEN
  440. * la statistique est de dimension 3
  441. LADIM = 3
  442.  
  443. C Lecture de la valeur de lambda3 (strictement supérieure à 0.D0)
  444. CALL LIRREE(CLAMD3,1,IRET)
  445. IF (IERR.NE.0) RETURN
  446. IF (IRET.EQ.0) THEN
  447. MOTERR(1:4) = 'LAMB'
  448. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  449. CALL ERREUR(166)
  450. RETURN
  451. ENDIF
  452. IF (CLAMD3.LE.0.D0) THEN
  453. REAERR(1) = REAL(CLAMD3)
  454. REAERR(2) = REAL(0.D0)
  455. MOTERR(1:8) = 'LAMBDA3 '
  456. c %m1:8 = %r1 inférieur à %r2
  457. CALL ERREUR(41)
  458. RETURN
  459. ENDIF
  460.  
  461. C Lecture obligatoire de la direction VEC3 si on a donné
  462. C VEC1 et VEC2
  463. IF (AXES) THEN
  464. CALL LIROBJ('POINT',IPTV,1,IRET)
  465. IF (IERR.NE.0) RETURN
  466.  
  467. IF (IRET.EQ.0) THEN
  468. * Syntaxe incorrecte : on attend %m1:30
  469. MOTERR(1:30) = 'un vecteur '
  470. CALL ERREUR(881)
  471. RETURN
  472. ENDIF
  473.  
  474. C on charge le vecteur et on vérifie qu'il n'est pas de longueur
  475. C nulle, et qu'il est bien perpendiculaire à VEC1 et à VEC2.
  476. C on le normalise aussi.
  477. IREF = (IPTV-1) * (IDIM+1)
  478. SDIR = 0.D0
  479. PSCL = 0.D0
  480. PSCL2 = 0.D0
  481. DO 15 I=1,IDIM
  482. XDIR3(I) = XCOOR(IREF+I)
  483. SDIR = SDIR + (XDIR3(I) * XDIR3(I))
  484. PSCL = PSCL + (XDIR1(I) * XDIR3(I))
  485. PSCL2 = PSCL + (XDIR2(I) * XDIR3(I))
  486. 15 CONTINUE
  487. SDIR = SQRT(SDIR)
  488. * si VEC3 de longueur nulle
  489. IF (SDIR.LT.EPS) THEN
  490. c Une direction ne peut pas être définie par un vecteur nul
  491. CALL ERREUR(239)
  492. RETURN
  493. ENDIF
  494. * si VEC3 non perpendiculaire à VEC1 et VEC2
  495. IF (((ABS(PSCL )).GE.(EPS*SDIR*SDIR)).OR.
  496. & ((ABS(PSCL2)).GE.(EPS*SDIR*SDIR))) THEN
  497. c Les vecteurs définissant le repère local ne sont pas
  498. * orthogonaux
  499. CALL ERREUR(161)
  500. RETURN
  501. ENDIF
  502. * orthogonalisation à VEC1 :
  503. DO I=1,IDIM
  504. XDIR2(I)=XDIR2(I) - (PSCL*XDIR1(I)) - (PSCL2*XDIR2(I))
  505. ENDDO
  506. * normalisation
  507. DO I=1,IDIM
  508. XDIR3(I) = XDIR3(I) / SDIR
  509. ENDDO
  510. ENDIF lit vec3
  511. ENDIF stat 3D
  512. ENDIF dime 3D
  513. ENDIF stat 2D
  514. ENDIF anisotrope
  515. c
  516. c ================ Travail ================
  517. c
  518. C Dans la suite, les coordonnées vont être adimensionnées par
  519. C Lambda. Donc toutes les données numériques seront calculées comme
  520. C si Lambda est isotrope égal à 1.
  521.  
  522. C Longueur minimale de description des hétérogénéités
  523. C pour déterminer la plage de variation des vecteurs d'onde
  524. ZDIST = 0.2D0
  525. * Pas de discrétisation de la bande (avec lambda = 1)
  526. IF (LADIM.EQ.1) THEN
  527. DELZET = ZDIST / 3.D0
  528. ELSE
  529. DELZET = ZDIST / 4.D0
  530. ENDIF
  531.  
  532. * Fréquence de coupure (avec lambda = 1),
  533. * elle doit être supérieure à 2 pi / dx
  534. OMMAX = 1.1D0 * (2.D0 * XPI / ZDIST)
  535.  
  536. * Génération
  537. IF (IOBJ.EQ.1) THEN
  538. * On veut un CHAMELEM en sortie
  539. CALL ALEA1(MMODEL,IPOSI,LADIM,XDIR1,XDIR2,XDIR3,
  540. & ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
  541. ELSEIF (IOBJ.EQ.2) THEN
  542. * On veut un CHAMPOINT en sortie
  543. CALL ALEA2(MELENT,LADIM,XDIR1,XDIR2,XDIR3,
  544. & ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
  545. ENDIF
  546. END
  547.  
  548.  
  549.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales