Télécharger alea.eso

Retour à la liste

Numérotation des lignes :

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

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