Télécharger alea1.eso

Retour à la liste

Numérotation des lignes :

  1. C ALEA1 SOURCE FD218221 17/01/11 21:15:01 9277
  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. IFOCHE = 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
  126. CALL LEKMOD(MMODEL,IPTABL,IRET)
  127. ENDIF
  128. ENDIF
  129.  
  130. IF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN
  131. * On vérifie que le modèle est bien un modèle MECANIQUE
  132. * sinon, on passe à la partition suivante.
  133. IF ((NOMFOR.EQ.'NAVIER_STOKES').OR.(NOMFOR.EQ.'DARCY')
  134. & .AND.(NOMFOR.NE.'EULER')) THEN
  135. LPART(IPART) = .FALSE.
  136. GOTO 4
  137. ENDIF
  138. ENDIF
  139.  
  140. IF ((IPOSI.EQ.7).OR.(IPOSI.EQ.8)) THEN
  141. * On vérifie que le modèle est bien un modèle NAVIER-STOKES
  142. * ou EULER sinon, on passe à la partition suivante.
  143. IF ((NOMFOR.NE.'NAVIER_STOKES').AND.(NOMFOR.NE.'EULER')) THEN
  144. LPART(IPART) = .FALSE.
  145. GOTO 4
  146. ENDIF
  147. ENDIF
  148. *
  149. * N° d'élément fini NELE
  150. * nombre d'éléments NBEL
  151. * nombre de noeuds par éléments NBPT
  152. * nombre de fonctions de forme NBNN (pts de gauss seulement)
  153. NELE = NEFMOD
  154.  
  155. MELENT = IMAMOD
  156. SEGACT MELENT
  157. NBPT = MELENT.NUM(/1)
  158. NBEL = MELENT.NUM(/2)
  159.  
  160. * ------------------------------------------------------------------
  161. * 2a. Détermination du support des points de valeur MELSUP s'ils existent,
  162. * du nombre de valeurs à stocker par élément NBPGAU,
  163. * et des pointeurs sur les segments d'intégration MINTE
  164.  
  165. IF (IPOSI.EQ.1) THEN
  166. * Noeuds
  167.  
  168. * On utilise tous les noeuds du maillage fourni
  169. * En particulier, en MECA-FLUX, tous les noeuds des QUAF
  170. MELSUP = MELENT
  171. SEGACT MELSUP
  172. NBPGAU = MELSUP.NUM(/1)
  173. SEGDES MELSUP
  174. MINTE = 0
  175.  
  176. ELSEIF(IPOSI.EQ.2) THEN
  177. * Points centres
  178.  
  179. * on les tire du maillage QUAF par la table domaine
  180. CALL LEKTAB(IPTABL,'CENTRE',MELSUP)
  181. IF (IERR.NE.0) THEN
  182. * On ne trouve pas le support qui contient les points
  183. CALL ERREUR(248)
  184. RETURN
  185. ENDIF
  186. if(infmod(/1).lt.iposi+2) then
  187. CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL)
  188. MINTE=infell(11)
  189. SEGSUP INFO
  190. else
  191. MINTE = infmod(4)
  192. endif
  193. NBPGAU = 1
  194.  
  195. ELSEIF(IPOSI.EQ.6) THEN
  196. * Points faces
  197.  
  198. * !!!!!!!!!! La structure de Castem ne permet pas encore de
  199. * créer un CHAMELEM appuyé aux faces.
  200. * Option noeud indisponible
  201. CALL ERREUR(507)
  202. RETURN
  203.  
  204. ** on les tire du maillage QUAF par la table domaine
  205. * CALL LEKTAB(IPTABL,'FACE',MELSUP)
  206. * IF (IERR.NE.0) THEN
  207. ** On ne trouve pas le support qui contient les points
  208. * CALL ERREUR(248)
  209. * RETURN
  210. * ENDIF
  211. * CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL)
  212. * MINTE = INFELE(11)
  213. * SEGSUP INFO
  214. * SEGACT MELSUP
  215. * NBPGAU = MELSUP.NUM(/1)
  216. * SEGDES MELSUP
  217.  
  218. ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN
  219. * Points de gauss MECANIQUE
  220.  
  221. * CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL)
  222. IF(IPOSI.EQ.3) then
  223. NBPGAU = INFELE(6)
  224. minte=infmod(5)
  225. ELSEIF(IPOSI.EQ.4) then
  226. NBPGAU = INFELE(3)
  227. minte=infmod(6)
  228. ELSEIF(IPOSI.EQ.5) then
  229. NBPGAU = INFELE(4)
  230. MINTE = INFMOD(7)
  231. ENDIF
  232. NBNN = INFELE(8)
  233. * SEGSUP INFO
  234.  
  235. ELSEIF ((IPOSI.EQ.7).OR.(IPOSI.EQ.8)) THEN
  236. * !!!!!!! Implémentation possible ci-après pour
  237. * les points de gauss du modèle NAVIER-STOKES (VITESSE et PRESSION)
  238.  
  239. * Option noeud indisponible (pour l'instant)
  240. CALL ERREUR(507)
  241. RETURN
  242.  
  243. ELSE
  244. * Erreur anormale.contactez votre support
  245. CALL ERREUR(5)
  246. ENDIF
  247. SEGDES MELENT
  248. *
  249. * ------------------------------------------------------------------
  250. * 2b. Contenu par partition du CHAMELEM :
  251. *
  252. * on ajoute une partition au CHAMELEM
  253. N1 = N1 + 1
  254. SEGADJ MCHELM
  255.  
  256. * Le maillage pointé par le CHAMELEM est le maillage du modèle.
  257.  
  258. IF (NOMFOR.EQ.'DARCY') THEN
  259. * !!!!!!!! les tests des calculs DARCY exigent encore le
  260. * maillage de sommets, récupéré par la table domaine, alors
  261. * que le modèle pointe sur le maillage de QUAF.
  262. CALL LEKMOD(MMODEL,IPTABL,IRET)
  263. IF (IERR.NE.0) THEN
  264. C Données incompatibles
  265. CALL ERREUR(21)
  266. RETURN
  267. ENDIF
  268. * la table domaine existe. Le CHAMELEM va pointer sur le maillage
  269. * de sommets. On suppose bien entendu que les partitions du
  270. * modèle sont parcourues dans le même ordre que les
  271. * sous-maillages du maillage.
  272. CALL LEKTAB(IPTABL,'MAILLAGE',IPT1)
  273. IF (IERR.NE.0) RETURN
  274. SEGACT IPT1
  275. IF (IPT1.LISOUS(/1).EQ.0) THEN
  276. MELCHM = IPT1
  277. ELSE
  278. MELCHM = IPT1.LISOUS(IPART)
  279. ENDIF
  280. SEGDES IPT1
  281. ELSE
  282. * Le CHAMELEM pointe sur le maillage du modèle
  283. * (NAVIER-STOKES ou MECANIQUE).
  284. MELCHM = MELENT
  285. ENDIF
  286.  
  287. CONCHE(N1) = CONMOD
  288. IMACHE(N1) = MELCHM
  289. INFCHE(N1,1) = 0
  290. INFCHE(N1,2) = 0
  291. INFCHE(N1,3) = NIFOUR
  292. INFCHE(N1,4) = MINTE
  293. INFCHE(N1,5) = 0
  294. INFCHE(N1,6) = IPOSI
  295. *
  296. * ------------------------------------------------------------------
  297. * 2c. Construction des coordonnées des points supports :
  298. * dans le repère réel.
  299. *
  300. * combien de points de valeur à calculer en tout ?
  301. L = NBL
  302. NBL = NBL + (NBEL * NBPGAU)
  303. SEGADJ TABCOR
  304.  
  305. * Cas où les points existent, et leurs coordonnées sont dans XCOOR.
  306. IF ((IPOSI.EQ.1).OR.(IPOSI.EQ.2).OR.(IPOSI.EQ.6)) THEN
  307.  
  308. * Ce sont les points sommets, centres ou faces.
  309. SEGACT MELSUP
  310. DO 9 I=1,NBEL
  311. DO 10 J=1,NBPGAU
  312. L = L + 1
  313. IREF = (IDIM+1) * (MELSUP.NUM(J,I)-1)
  314. DO 11 K=1,IDIM
  315. COR(L,K) = XCOOR(IREF+K)
  316. 11 CONTINUE
  317. 10 CONTINUE
  318. 9 CONTINUE
  319. SEGDES MELSUP
  320.  
  321. ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN
  322.  
  323. * Les points n'existent pas. Ce sont les points de gauss MECANIQUE.
  324. * On calcule leurs coordonnées à partir des fonctions de forme
  325. * et des coordonnées des points sommets associés.
  326. SEGACT MELENT
  327. SEGACT MINTE
  328. NBBB = NBEL
  329. SEGINI WRK2
  330.  
  331. DO 5 IEL=1,NBEL
  332.  
  333. * Coordonnées XE des points supports des fonctions de forme
  334. * de l'élément IEL
  335. * Le maillage MELENT contient nécessairement au moins ces points
  336. CALL DOXE(XCOOR,IDIM,NBNN,MELENT.NUM,IEL,XE)
  337. *
  338. * Boucle sur les points de gauss
  339. DO 1100 IGAU=1,NBPGAU
  340. *
  341. L = L + 1
  342.  
  343. XX=0.D0
  344. YY=0.D0
  345. ZZ=0.D0
  346. DO 1101 ISH=1,NBNN
  347. * valeur de la fonction de forme n°ISH au point de gauss IGAU
  348. CGAUSS = SHPTOT(1,ISH,IGAU)
  349. * coordonnée réelle du point support des fonctions de forme
  350. * coordonnées réelles du point de gauss
  351. XX = XX + XE(1,ISH)*CGAUSS
  352. YY = YY + XE(2,ISH)*CGAUSS
  353. IF (IDIM.EQ.3) THEN
  354. ZZ = ZZ + XE(3,ISH)*CGAUSS
  355. ENDIF
  356. 1101 CONTINUE
  357.  
  358. COR(L,1) = XX
  359. COR(L,2) = YY
  360. COR(L,3) = ZZ
  361.  
  362. 1100 CONTINUE
  363.  
  364. 5 CONTINUE
  365. SEGDES MELENT,MINTE,WRK2
  366.  
  367. ELSE
  368.  
  369. * Erreur anormale.contactez votre support
  370. CALL ERREUR(5)
  371. RETURN
  372. ENDIF
  373. *
  374. * ------------------------------------------------------------------
  375. * 2d. Transformation des coordonnées dans le repère adimensionné
  376. * par la matrice lambda de dimension LADIM.
  377. *
  378. IF (IDIM.EQ.2) THEN
  379. IF (LADIM.EQ.1) THEN
  380. * 2D stat 1D
  381. DO 20 L=1,NBL
  382. XX = COR(L,1)
  383. YY = COR(L,2)
  384. COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1
  385. 20 CONTINUE
  386. ELSE
  387. * 2D stat 2D
  388. DO 21 L=1,NBL
  389. XX = COR(L,1)
  390. YY = COR(L,2)
  391. COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1
  392. COR(L,2) = ((XX * XDIR2(1)) + (YY * XDIR2(2))) / CLAMD2
  393. 21 CONTINUE
  394. ENDIF
  395. ELSE
  396. IF (LADIM.EQ.1) THEN
  397. * 3D stat 1D
  398. DO 22 L=1,NBL
  399. XX = COR(L,1)
  400. YY = COR(L,2)
  401. ZZ = COR(L,3)
  402. COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2))
  403. & + (ZZ * XDIR1(3)) ) / CLAMD1
  404. 22 CONTINUE
  405. ELSEIF (LADIM.EQ.2) THEN
  406. * 3D stat 2D
  407. DO 23 L=1,NBL
  408. XX = COR(L,1)
  409. YY = COR(L,2)
  410. ZZ = COR(L,3)
  411. COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2))
  412. & + (ZZ * XDIR1(3)) ) / CLAMD1
  413. COR(L,2) = ( (XX * XDIR2(1)) + (YY * XDIR2(2))
  414. & + (ZZ * XDIR2(3)) ) / CLAMD2
  415. 23 CONTINUE
  416. ELSE
  417. * 3D stat 3D
  418. DO 24 L=1,NBL
  419. XX = COR(L,1)
  420. YY = COR(L,2)
  421. ZZ = COR(L,3)
  422. COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2))
  423. & + (ZZ * XDIR1(3)) ) / CLAMD1
  424. COR(L,2) = ( (XX * XDIR2(1)) + (YY * XDIR2(2))
  425. & + (ZZ * XDIR2(3)) ) / CLAMD2
  426. COR(L,3) = ( (XX * XDIR3(1)) + (YY * XDIR3(2))
  427. & + (ZZ * XDIR3(3)) ) / CLAMD3
  428. 24 CONTINUE
  429. ENDIF
  430. ENDIF
  431.  
  432. *
  433. * ------------------------------------------------------------------
  434. * 2e. Contenu par composante (1 seule, 'SCAL') du CHAMELEM :
  435. *
  436. N2 = 1
  437. SEGINI MCHAML
  438. ICHAML(N1) = MCHAML
  439. NOMCHE(1) = 'SCAL '
  440. TYPCHE(1) = 'REAL*8 '
  441. N1PTEL = NBPGAU
  442. N1EL = NBEL
  443. N2PTEL = 0
  444. N2EL = 0
  445. SEGINI MELVAL
  446. IELVAL(1) = MELVAL
  447.  
  448. SEGDES IMODEL
  449.  
  450. 4 CONTINUE
  451. SEGDES MMODEL
  452.  
  453. IF (NBL.EQ.0) THEN
  454. * on n'a pas trouvé de partition avec le modèle correspondant
  455. * aux points de gauss demandés :
  456. C Données incompatibles
  457. CALL ERREUR(21)
  458. RETURN
  459. ENDIF
  460. *
  461. * ------------------------------------------------------------------
  462. * 3. Génération aléatoire :
  463. * ------------------------------------------------------------------
  464. *
  465. * Si les points supports appartiennent à différentes mailles en
  466. * même temps, il faut savoir que les valeurs en des points
  467. * confondus sont égales avec la méthode des bandes tournantes.
  468.  
  469. CALL BANTOU(TABCOR,TABVAL,LADIM,ZSIG,VALMOY,DELZET,OMMAX)
  470. SEGDES TABCOR
  471. *
  472. * ------------------------------------------------------------------
  473. * 4. Construction champ résultat :
  474. * ------------------------------------------------------------------
  475. *
  476. * BANTOU renvoie le résultat sous la forme d'une table de valeur,
  477. * TABVAL active, que l'on doit maintenant convertir en MELVAL,
  478. * partition par partition (on parcourt les partitions valides).
  479.  
  480. L = 1
  481. ICHELM = 0
  482. DO 8 IPART=1,NBPART
  483. IF (LPART(IPART)) THEN
  484. ICHELM = ICHELM + 1
  485. MCHAML = ICHAML(ICHELM)
  486. MELVAL = IELVAL(1)
  487.  
  488. NBEL = VELCHE(/2)
  489. NBPGAU = VELCHE(/1)
  490. DO 1 I=1,NBEL
  491. DO 19 J=1,NBPGAU
  492. VELCHE(J,I) = VAL(L)
  493. L = L + 1
  494. 19 CONTINUE
  495. 1 CONTINUE
  496. ENDIF
  497. 8 CONTINUE
  498.  
  499. * ------------------------------------------------------------------
  500. * 5. Ecriture et sortie
  501. * ------------------------------------------------------------------
  502. *
  503. * Fermeture des segments
  504. *
  505. DO 12 IPART=1,N1
  506. MCHAML = ICHAML(IPART)
  507. MELVAL = IELVAL(1)
  508. SEGDES MELVAL,MCHAML
  509. 12 CONTINUE
  510. SEGDES MCHELM
  511. SEGSUP TABVAL,TPART
  512. *
  513. * Résultat
  514. *
  515. CALL ECROBJ('MCHAML ',MCHELM)
  516.  
  517. RETURN
  518. END
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  

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