Télécharger alea.eso

Retour à la liste

Numérotation des lignes :

alea
  1. C ALEA SOURCE CB215821 23/01/25 21:15:03 11573
  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. SEGACT,MCOORD
  235. IF (ILAM.EQ.1) THEN
  236. C Cas ISOTROPE
  237.  
  238. * la statistique a la même dimension que l'espace réel :
  239. LADIM = IDIM
  240.  
  241. C Lecture de la valeur de lambda (strictement supérieure à 0.D0)
  242. CALL LIRREE(CLAMDA,1,IRET)
  243. IF (IERR.NE.0) RETURN
  244. IF (IRET.EQ.0) THEN
  245. MOTERR(1:4) = 'LAMB'
  246. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  247. CALL ERREUR(166)
  248. RETURN
  249. ENDIF
  250. IF (CLAMDA.LE.0.D0) THEN
  251. REAERR(1) = REAL(CLAMDA)
  252. REAERR(2) = REAL(0.D0)
  253. MOTERR(1:8) = 'LAMBDA '
  254. c %m1:8 = %r1 inférieur à %r2
  255. CALL ERREUR(41)
  256. RETURN
  257. ENDIF
  258.  
  259. * mêmes longueurs de corrélation dans toutes les directions
  260. CLAMD1 = CLAMDA
  261. CLAMD2 = CLAMDA
  262. CLAMD3 = CLAMDA
  263. * les vecteurs principaux sont les axes x, y et z
  264. XDIR1(1) = 1.D0
  265. XDIR1(2) = 0.D0
  266. XDIR1(3) = 0.D0
  267. XDIR2(1) = 0.D0
  268. XDIR2(2) = 1.D0
  269. XDIR2(3) = 0.D0
  270. XDIR3(1) = 0.D0
  271. XDIR3(2) = 0.D0
  272. XDIR3(3) = 1.D0
  273.  
  274. ELSE
  275. C Cas ANISOTROPE
  276.  
  277. C axe 1 :
  278. C on doit avoir lu le mot-clef 'LAMBDA1'
  279. IF (ILAM.NE.2) THEN
  280. * Syntaxe incorrecte : on attend %m1:30
  281. MOTERR(1:30) = 'le mot-clef LAMBDA1 '
  282. CALL ERREUR(881)
  283. RETURN
  284. ENDIF
  285. C Lecture de la valeur de lambda1 (strictement supérieure à 0.D0)
  286. CALL LIRREE(CLAMD1,1,IRET)
  287. IF (IERR.NE.0) RETURN
  288. IF (IRET.EQ.0) THEN
  289. MOTERR(1:4) = 'LAMB'
  290. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  291. CALL ERREUR(166)
  292. RETURN
  293. ENDIF
  294. IF (CLAMD1.LE.0.D0) THEN
  295. REAERR(1) = REAL(CLAMD1)
  296. REAERR(2) = REAL(0.D0)
  297. MOTERR(1:8) = 'LAMBDA1 '
  298. c %m1:8 = %r1 inférieur à %r2
  299. CALL ERREUR(41)
  300. RETURN
  301. ENDIF
  302.  
  303. C Lecture optionnelle de la direction VEC1
  304. CALL LIROBJ('POINT',IPTV,0,IRET)
  305. IF (IERR.NE.0) RETURN
  306.  
  307. * si on ne donne pas d'axe, ce sont les axes par défaut, et la
  308. * statistique a la même dimensinonalité que l'espace réel.
  309. AXES = (IRET.NE.0)
  310.  
  311. IF (.NOT.AXES) THEN
  312. * pas de vecteurs => axes par défaut.
  313. XDIR1(1) = 1.D0
  314. XDIR1(2) = 0.D0
  315. XDIR1(3) = 0.D0
  316. XDIR2(1) = 0.D0
  317. XDIR2(2) = 1.D0
  318. XDIR2(3) = 0.D0
  319. XDIR3(1) = 0.D0
  320. XDIR3(2) = 0.D0
  321. XDIR3(3) = 1.D0
  322. * la statistique a la même dimension que l'espace réel :
  323. LADIM = IDIM
  324. ELSE
  325. C on charge le vecteur et on vérifie qu'il n'est pas de longueur nulle
  326. C on le normalise aussi.
  327. IREF = (IPTV-1) * (IDIM+1)
  328. SDIR = 0.D0
  329. DO 13 I=1,IDIM
  330. XDIR1(I) = XCOOR(IREF+I)
  331. SDIR = SDIR + (XDIR1(I) * XDIR1(I))
  332. 13 CONTINUE
  333. SDIR = SQRT(SDIR)
  334. C sinon, erreur :
  335. IF (SDIR.LT.EPS) THEN
  336. c Une direction ne peut pas être définie par un vecteur nul
  337. CALL ERREUR(239)
  338. RETURN
  339. ENDIF
  340. * normalisation
  341. DO I=1,IDIM
  342. XDIR1(I) = XDIR1(I) / SDIR
  343. ENDDO
  344. * la statistique est au moins de dimension 1
  345. * (à surcharger si nécessaire) :
  346. LADIM = 1
  347. ENDIF
  348.  
  349. C axe 2 :
  350. C Si le mot-clef 'LAMBDA2' existe, stat 2D au moins
  351. CALL LIRMOT(MOTLA,4,ILAM,0)
  352. IF (IERR.NE.0) RETURN
  353.  
  354. IF ((ILAM.NE.0).AND.(ILAM.NE.3)) THEN
  355. * Syntaxe incorrecte : on attend %m1:30
  356. MOTERR(1:30) = 'le mot-clef LAMBDA2 '
  357. CALL ERREUR(881)
  358. RETURN
  359. ENDIF
  360. IF (ILAM.EQ.3) THEN
  361. * la statistique est au moins de dimension 2
  362. LADIM = 2
  363.  
  364. C Lecture de la valeur de lambda2 (strictement supérieure à 0.D0)
  365. CALL LIRREE(CLAMD2,1,IRET)
  366. IF (IERR.NE.0) RETURN
  367. IF (IRET.EQ.0) THEN
  368. MOTERR(1:4) = 'LAMB'
  369. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  370. CALL ERREUR(166)
  371. RETURN
  372. ENDIF
  373. IF (CLAMD2.LE.0.D0) THEN
  374. REAERR(1) = REAL(CLAMD2)
  375. REAERR(2) = REAL(0.D0)
  376. MOTERR(1:8) = 'LAMBDA2 '
  377. c %m1:8 = %r1 inférieur à %r2
  378. CALL ERREUR(41)
  379. RETURN
  380. ENDIF
  381.  
  382. C Lecture obligatoire de la direction VEC2 si on a donné VEC1
  383. IF (AXES) THEN
  384. CALL LIROBJ('POINT',IPTV,1,IRET)
  385. IF (IERR.NE.0) RETURN
  386.  
  387. IF (IRET.EQ.0) THEN
  388. * Syntaxe incorrecte : on attend %m1:30
  389. MOTERR(1:30) = 'un vecteur '
  390. CALL ERREUR(881)
  391. RETURN
  392. ENDIF
  393.  
  394. C on charge le vecteur et on vérifie qu'il n'est pas de longueur
  395. C nulle, et qu'il est bien perpendiculaire à VEC1.
  396. C on le normalise aussi.
  397. IREF = (IPTV-1) * (IDIM+1)
  398. SDIR = 0.D0
  399. PSCL = 0.D0
  400. DO 14 I=1,IDIM
  401. XDIR2(I) = XCOOR(IREF+I)
  402. SDIR = SDIR + (XDIR2(I) * XDIR2(I))
  403. PSCL = PSCL + (XDIR1(I) * XDIR2(I))
  404. 14 CONTINUE
  405. SDIR = SQRT(SDIR)
  406. * si VEC2 de longueur nulle
  407. IF (SDIR.LT.EPS) THEN
  408. c Une direction ne peut pas être définie par un vecteur nul
  409. CALL ERREUR(239)
  410. RETURN
  411. ENDIF
  412. * si VEC2 non perpendiculaire à VEC1
  413. IF ((ABS(PSCL)).GE.(EPS*SDIR*SDIR)) THEN
  414. c Les vecteurs définissant le repère local ne sont pas orthogonaux
  415. CALL ERREUR(161)
  416. RETURN
  417. ENDIF
  418. * orthogonalisation à VEC1 :
  419. DO I=1,IDIM
  420. XDIR2(I) = XDIR2(I) - (PSCL * XDIR1(I))
  421. ENDDO
  422. * normalisation
  423. DO I=1,IDIM
  424. XDIR2(I) = XDIR2(I) / SDIR
  425. ENDDO
  426. ENDIF
  427.  
  428. C axe 3 (si 3D) :
  429. IF (IDIM.EQ.3) THEN
  430. C Si le mot-clef 'LAMBDA3' existe, stat 3D
  431. CALL LIRMOT(MOTLA,4,ILAM,0)
  432. IF (IERR.NE.0) RETURN
  433.  
  434. IF ((ILAM.NE.0).AND.(ILAM.NE.4)) THEN
  435. * Syntaxe incorrecte : on attend %m1:30
  436. MOTERR(1:30) = 'le mot-clef LAMBDA3 '
  437. CALL ERREUR(881)
  438. RETURN
  439. ENDIF
  440. IF (ILAM.EQ.4) THEN
  441. * la statistique est de dimension 3
  442. LADIM = 3
  443.  
  444. C Lecture de la valeur de lambda3 (strictement supérieure à 0.D0)
  445. CALL LIRREE(CLAMD3,1,IRET)
  446. IF (IERR.NE.0) RETURN
  447. IF (IRET.EQ.0) THEN
  448. MOTERR(1:4) = 'LAMB'
  449. c Le mot-clé %m1:4 n'est pas suivi de la donnée correspondante
  450. CALL ERREUR(166)
  451. RETURN
  452. ENDIF
  453. IF (CLAMD3.LE.0.D0) THEN
  454. REAERR(1) = REAL(CLAMD3)
  455. REAERR(2) = REAL(0.D0)
  456. MOTERR(1:8) = 'LAMBDA3 '
  457. c %m1:8 = %r1 inférieur à %r2
  458. CALL ERREUR(41)
  459. RETURN
  460. ENDIF
  461.  
  462. C Lecture obligatoire de la direction VEC3 si on a donné
  463. C VEC1 et VEC2
  464. IF (AXES) THEN
  465. CALL LIROBJ('POINT',IPTV,1,IRET)
  466. IF (IERR.NE.0) RETURN
  467.  
  468. IF (IRET.EQ.0) THEN
  469. * Syntaxe incorrecte : on attend %m1:30
  470. MOTERR(1:30) = 'un vecteur '
  471. CALL ERREUR(881)
  472. RETURN
  473. ENDIF
  474.  
  475. C on charge le vecteur et on vérifie qu'il n'est pas de longueur
  476. C nulle, et qu'il est bien perpendiculaire à VEC1 et à VEC2.
  477. C on le normalise aussi.
  478. IREF = (IPTV-1) * (IDIM+1)
  479. SDIR = 0.D0
  480. PSCL = 0.D0
  481. PSCL2 = 0.D0
  482. DO 15 I=1,IDIM
  483. XDIR3(I) = XCOOR(IREF+I)
  484. SDIR = SDIR + (XDIR3(I) * XDIR3(I))
  485. PSCL = PSCL + (XDIR1(I) * XDIR3(I))
  486. PSCL2 = PSCL + (XDIR2(I) * XDIR3(I))
  487. 15 CONTINUE
  488. SDIR = SQRT(SDIR)
  489. * si VEC3 de longueur nulle
  490. IF (SDIR.LT.EPS) THEN
  491. c Une direction ne peut pas être définie par un vecteur nul
  492. CALL ERREUR(239)
  493. RETURN
  494. ENDIF
  495. * si VEC3 non perpendiculaire à VEC1 et VEC2
  496. IF (((ABS(PSCL )).GE.(EPS*SDIR*SDIR)).OR.
  497. & ((ABS(PSCL2)).GE.(EPS*SDIR*SDIR))) THEN
  498. c Les vecteurs définissant le repère local ne sont pas
  499. * orthogonaux
  500. CALL ERREUR(161)
  501. RETURN
  502. ENDIF
  503. * orthogonalisation à VEC1 :
  504. DO I=1,IDIM
  505. XDIR2(I)=XDIR2(I) - (PSCL*XDIR1(I)) - (PSCL2*XDIR2(I))
  506. ENDDO
  507. * normalisation
  508. DO I=1,IDIM
  509. XDIR3(I) = XDIR3(I) / SDIR
  510. ENDDO
  511. ENDIF lit vec3
  512. ENDIF stat 3D
  513. ENDIF dime 3D
  514. ENDIF stat 2D
  515. ENDIF anisotrope
  516. c
  517. c ================ Travail ================
  518. c
  519. C Dans la suite, les coordonnées vont être adimensionnées par
  520. C Lambda. Donc toutes les données numériques seront calculées comme
  521. C si Lambda est isotrope égal à 1.
  522.  
  523. C Longueur minimale de description des hétérogénéités
  524. C pour déterminer la plage de variation des vecteurs d'onde
  525. ZDIST = 0.2D0
  526. * Pas de discrétisation de la bande (avec lambda = 1)
  527. IF (LADIM.EQ.1) THEN
  528. DELZET = ZDIST / 3.D0
  529. ELSE
  530. DELZET = ZDIST / 4.D0
  531. ENDIF
  532.  
  533. * Fréquence de coupure (avec lambda = 1),
  534. * elle doit être supérieure à 2 pi / dx
  535. OMMAX = 1.1D0 * (2.D0 * XPI / ZDIST)
  536.  
  537. * Génération
  538. IF (IOBJ.EQ.1) THEN
  539. * On veut un CHAMELEM en sortie
  540. CALL ALEA1(MMODEL,IPOSI,LADIM,XDIR1,XDIR2,XDIR3,
  541. & ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
  542. ELSEIF (IOBJ.EQ.2) THEN
  543. * On veut un CHAMPOINT en sortie
  544. CALL ALEA2(MELENT,LADIM,XDIR1,XDIR2,XDIR3,
  545. & ZSIG,CLAMD1,CLAMD2,CLAMD3,VALMOY,DELZET,OMMAX)
  546. ENDIF
  547. SEGDES,MCOORD
  548. END
  549.  
  550.  
  551.  
  552.  

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