Télécharger alea1.eso

Retour à la liste

Numérotation des lignes :

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

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