Télécharger alea1.eso

Retour à la liste

Numérotation des lignes :

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

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