Télécharger la plaquette Cast3M


fr en

Télécharger alea1.eso

Retour à la liste

Numérotation des lignes : 

alea1

  1. C ALEA1     SOURCE    MB234859  25/09/08    21:15:03     12358          
  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     Coordonnées reduites et fonctions de forme pour trouver
  62. C     les coordonnées réelles des points de gauss
  63.       SEGMENT WRK2
  64.         REAL*8 XE(3,NBBB),SHP(6,NBBB)
  65.       ENDSEGMENT
  66.  
  67.       DIMENSION    XDIR1(3),XDIR2(3),XDIR3(3)
  68.       CHARACTER*16 NOMFOR
  69.       POINTEUR     MELCHP.MELEME, MELENT.MELEME, MELSUP.MELEME
  70.  
  71.       INTEGER      IPOSI
  72. C
  73. C     epsilon servant à différents tests
  74.       EPS = 1.D-12
  75. *
  76. *     ------------------------------------------------------------------
  77. *     2. On construit les coordonnées des points supports
  78. *        on initialise le champ résultat en parallèle (CHAMELEM)
  79. *     ------------------------------------------------------------------
  80.  
  81.       SEGACT MMODEL
  82. *     Nombre de partitions :
  83.       NBPART = KMODEL(/1)
  84.  
  85. *     En-tête du CHAMELEM
  86. *     Le nombre de partitions est fixé à zéro pour l'instant, et
  87. *     il va augmenter à chaque partition invalide du modèle.
  88.       L1 = 40
  89.       N1 = 0
  90.       N3 = 6
  91.       SEGINI MCHELM
  92.       TITCHE='   CHAMELEM cree par l operateur ALEA   '
  93.       IFOCHE = IFOUR
  94. *
  95. *     Boucle sur les partitions du modèle.
  96. *     Le chamelem résultat aura autant de partitions qu'il y a de
  97. *     sous-maillages supports correspondant.
  98. *
  99.       NBL = 1
  100.       SEGINI TABCOR
  101.       SEGINI TPART
  102.  
  103.       NBL = 0
  104.       NBL1 = 1
  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.           MINTE  = INFMOD(4)
  187.           NBPGAU = 1
  188.  
  189.         ELSEIF(IPOSI.EQ.6) THEN
  190. *         Points faces
  191.  
  192. *         !!!!!!!!!! La structure de Castem ne permet pas encore de
  193. *         créer un CHAMELEM appuyé aux faces.
  194. *         Option noeud indisponible
  195.           CALL ERREUR(507)
  196.           RETURN
  197.  
  198. **         on les tire du maillage QUAF par la table domaine
  199. *          CALL LEKTAB(IPTABL,'FACE',MELSUP)
  200. *          IF (IERR.NE.0) THEN
  201. **           On ne trouve pas le support qui contient les points
  202. *            CALL ERREUR(248)
  203. *            RETURN
  204. *          ENDIF
  205. *          CALL ELQUOI(NELE,0,IPOSI,INFO,IMODEL)
  206. *          MINTE  = INFELE(11)
  207. *          SEGSUP INFO
  208. *          SEGACT MELSUP
  209. *          NBPGAU = MELSUP.NUM(/1)
  210. *          SEGDES MELSUP
  211.  
  212.         ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN
  213. *         Points de gauss MECANIQUE
  214.           IF(IPOSI.EQ.3) then
  215.             NBPGAU = INFELE(6)
  216.           ELSEIF(IPOSI.EQ.4) then
  217.             NBPGAU = INFELE(3)
  218.           ELSEIF(IPOSI.EQ.5) then
  219.             NBPGAU = INFELE(4)
  220.           ENDIF
  221.           MINTE = INFMOD(2+IPOSI)
  222.           NBNN  = INFELE(8)
  223.  
  224.         ELSEIF ((IPOSI.EQ.7).OR.(IPOSI.EQ.8)) THEN
  225. *         !!!!!!! Implémentation possible ci-après pour
  226. *         les points de gauss du modèle NAVIER-STOKES (VITESSE et PRESSION)
  227.  
  228. *         Option noeud indisponible (pour l'instant)
  229.           CALL ERREUR(507)
  230.           RETURN
  231.  
  232.         ELSE
  233. *         Erreur anormale.contactez votre support
  234.           CALL ERREUR(5)
  235.         ENDIF
  236.         SEGDES MELENT
  237. *
  238. *       ------------------------------------------------------------------
  239. *       2b. Contenu par partition du CHAMELEM :
  240. *
  241. *       on ajoute une partition au CHAMELEM
  242.         N1 = N1 + 1
  243.         SEGADJ MCHELM
  244.  
  245. *       Le maillage pointé par le CHAMELEM est le maillage du modèle.
  246.  
  247.         IF (NOMFOR.EQ.'DARCY') THEN
  248. *         !!!!!!!! les tests des calculs DARCY exigent encore le
  249. *         maillage de sommets, récupéré par la table domaine, alors
  250. *         que le modèle pointe sur le maillage de QUAF.
  251.           CALL LEKMOD(MMODEL,IPTABL,IRET)
  252.           IF (IERR.NE.0) THEN
  253. C           Données incompatibles
  254.             CALL ERREUR(21)
  255.             RETURN
  256.           ENDIF
  257. *         la table domaine existe. Le CHAMELEM va pointer sur le maillage
  258. *         de sommets. On suppose bien entendu que les partitions du
  259. *         modèle sont parcourues dans le même ordre que les
  260. *         sous-maillages du maillage.
  261.           CALL LEKTAB(IPTABL,'MAILLAGE',IPT1)
  262.           IF (IERR.NE.0) RETURN
  263.           SEGACT IPT1
  264.           IF (IPT1.LISOUS(/1).EQ.0) THEN
  265.             MELCHM = IPT1
  266.           ELSE
  267.             MELCHM = IPT1.LISOUS(IPART)
  268.           ENDIF
  269.           SEGDES IPT1
  270.         ELSE
  271. *         Le CHAMELEM pointe sur le maillage du modèle
  272. *         (NAVIER-STOKES ou MECANIQUE).
  273.           MELCHM = MELENT
  274.         ENDIF
  275.  
  276.         CONCHE(N1)   = CONMOD
  277.         IMACHE(N1)   = MELCHM
  278.         INFCHE(N1,1) = 0
  279.         INFCHE(N1,2) = 0
  280.         INFCHE(N1,3) = NIFOUR
  281.         INFCHE(N1,4) = MINTE
  282.         INFCHE(N1,5) = 0
  283.         INFCHE(N1,6) = IPOSI
  284. *
  285. *       ------------------------------------------------------------------
  286. *       2c. Construction des coordonnées des points supports :
  287. *           dans le repère réel.
  288. *
  289. *       combien de points de valeur à calculer en tout ?
  290.         L    = NBL
  291.         NBL  = NBL + (NBEL * NBPGAU)
  292.         SEGADJ TABCOR
  293.  
  294. *       Cas où les points existent, et leurs coordonnées sont dans XCOOR.
  295.         IF ((IPOSI.EQ.1).OR.(IPOSI.EQ.2).OR.(IPOSI.EQ.6)) THEN
  296.  
  297. *         Ce sont les points sommets, centres ou faces.
  298.           SEGACT MELSUP
  299.           DO 9 I=1,NBEL
  300.             DO 10 J=1,NBPGAU
  301.               L = L + 1
  302.               IREF = (IDIM+1) * (MELSUP.NUM(J,I)-1)
  303.               DO 11 K=1,IDIM
  304.                 COR(L,K) = XCOOR(IREF+K)
  305.  11           CONTINUE
  306.  10         CONTINUE
  307.  9        CONTINUE
  308.           SEGDES MELSUP
  309.  
  310.         ELSEIF ((IPOSI.EQ.3).OR.(IPOSI.EQ.4).OR.(IPOSI.EQ.5)) THEN
  311.  
  312. *         Les points n'existent pas. Ce sont les points de gauss MECANIQUE.
  313. *         On calcule leurs coordonnées à partir des fonctions de forme
  314. *         et des coordonnées des points sommets associés.
  315.           SEGACT MELENT
  316.           SEGACT MINTE
  317.           NBBB = NBEL
  318.           SEGINI WRK2
  319.  
  320.           DO 5 IEL=1,NBEL
  321.  
  322. *           Coordonnées XE des points supports des fonctions de forme
  323. *           de l'élément IEL
  324. *           Le maillage MELENT contient nécessairement au moins ces points
  325.             CALL DOXE(XCOOR,IDIM,NBNN,MELENT.NUM,IEL,XE)
  326. *
  327. *           Boucle sur les points de gauss
  328.             DO 1100 IGAU=1,NBPGAU
  329. *
  330.               L = L + 1
  331.  
  332.               XX=0.D0
  333.               YY=0.D0
  334.               ZZ=0.D0
  335.               DO 1101 ISH=1,NBNN
  336. *               valeur de la fonction de forme n°ISH au point de gauss IGAU
  337.                 CGAUSS = SHPTOT(1,ISH,IGAU)
  338. *               coordonnée réelle du point support des fonctions de forme
  339. *               coordonnées réelles du point de gauss
  340.                 XX = XX + XE(1,ISH)*CGAUSS
  341.                 YY = YY + XE(2,ISH)*CGAUSS
  342.                 IF (IDIM.EQ.3) THEN
  343.                   ZZ = ZZ + XE(3,ISH)*CGAUSS
  344.                 ENDIF
  345.  1101         CONTINUE
  346.  
  347.               COR(L,1) = XX
  348.               COR(L,2) = YY
  349.               COR(L,3) = ZZ
  350.  
  351.  1100       CONTINUE
  352.  
  353.  5        CONTINUE
  354.           SEGDES MELENT,MINTE,WRK2
  355.  
  356.         ELSE
  357.  
  358. *         Erreur anormale.contactez votre support
  359.           CALL ERREUR(5)
  360.           RETURN
  361.         ENDIF
  362. *
  363. *       ------------------------------------------------------------------
  364. *       2d. Transformation des coordonnées dans le repère adimensionné
  365. *           par la matrice lambda de dimension LADIM.
  366. *
  367.         IF (IDIM.EQ.2) THEN
  368.           IF (LADIM.EQ.1) THEN
  369. *           2D stat 1D
  370.             DO 20 L=NBL1,NBL
  371.               XX = COR(L,1)
  372.               YY = COR(L,2)
  373.               COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1
  374.  20         CONTINUE
  375.           ELSE
  376. *           2D stat 2D
  377.             DO 21 L=NBL1,NBL
  378.               XX = COR(L,1)
  379.               YY = COR(L,2)
  380.               COR(L,1) = ((XX * XDIR1(1)) + (YY * XDIR1(2))) / CLAMD1
  381.               COR(L,2) = ((XX * XDIR2(1)) + (YY * XDIR2(2))) / CLAMD2
  382.  21         CONTINUE
  383.           ENDIF
  384.         ELSE
  385.           IF (LADIM.EQ.1) THEN
  386. *           3D stat 1D
  387.             DO 22 L=NBL1,NBL
  388.               XX = COR(L,1)
  389.               YY = COR(L,2)
  390.               ZZ = COR(L,3)
  391.               COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2))
  392.      &                   + (ZZ * XDIR1(3)) ) / CLAMD1
  393.  22         CONTINUE
  394.           ELSEIF (LADIM.EQ.2) THEN
  395. *           3D stat 2D
  396.             DO 23 L=NBL1,NBL
  397.               XX = COR(L,1)
  398.               YY = COR(L,2)
  399.               ZZ = COR(L,3)
  400.               COR(L,1) = ( (XX * XDIR1(1)) + (YY * XDIR1(2))
  401.      &                   + (ZZ * XDIR1(3)) ) / CLAMD1
  402.               COR(L,2) = ( (XX * XDIR2(1)) + (YY * XDIR2(2))
  403.      &                   + (ZZ * XDIR2(3)) ) / CLAMD2
  404.  23         CONTINUE
  405.           ELSE
  406. *           3D stat 3D
  407.             DO 24 L=NBL1,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.               COR(L,3) = ( (XX * XDIR3(1)) + (YY * XDIR3(2))
  416.      &                   + (ZZ * XDIR3(3)) ) / CLAMD3
  417.  24         CONTINUE
  418.           ENDIF
  419.         ENDIF
  420.         NBL1 = NBL
  421.  
  422. *
  423. *       ------------------------------------------------------------------
  424. *       2e. Contenu par composante (1 seule, 'SCAL') du CHAMELEM :
  425. *
  426.         N2 = 1
  427.         SEGINI MCHAML
  428.         ICHAML(N1) = MCHAML
  429.         NOMCHE(1)  = 'SCAL    '
  430.         TYPCHE(1)  = 'REAL*8          '
  431.         N1PTEL     = NBPGAU
  432.         N1EL       = NBEL
  433.         N2PTEL     = 0
  434.         N2EL       = 0
  435.         SEGINI MELVAL
  436.         IELVAL(1)  = MELVAL
  437.  
  438.         SEGDES IMODEL
  439.  
  440.  4    CONTINUE
  441.       SEGDES MMODEL
  442.  
  443.       IF (NBL.EQ.0) THEN
  444. *       on n'a pas trouvé de partition avec le modèle correspondant
  445. *       aux points de gauss demandés :
  446. C       Données incompatibles
  447.         CALL ERREUR(21)
  448.         RETURN
  449.       ENDIF
  450. *
  451. *     ------------------------------------------------------------------
  452. *     3. Génération aléatoire :
  453. *     ------------------------------------------------------------------
  454. *
  455. *     Si les points supports appartiennent à différentes mailles en
  456. *     même temps, il faut savoir que les valeurs en des points
  457. *     confondus sont égales avec la méthode des bandes tournantes.
  458.  
  459.       CALL BANTOU(TABCOR,TABVAL,LADIM,ZSIG,VALMOY,DELZET,OMMAX)
  460.       SEGDES TABCOR
  461. *
  462. *     ------------------------------------------------------------------
  463. *     4. Construction champ résultat :
  464. *     ------------------------------------------------------------------
  465. *
  466. *     BANTOU renvoie le résultat sous la forme d'une table de valeur,
  467. *     TABVAL active, que l'on doit maintenant convertir en MELVAL,
  468. *     partition par partition (on parcourt les partitions valides).
  469.  
  470.       L = 1
  471.       ICHELM = 0
  472.       DO 8 IPART=1,NBPART
  473.         IF (LPART(IPART)) THEN
  474.           ICHELM = ICHELM + 1
  475.           MCHAML = ICHAML(ICHELM)
  476.           MELVAL = IELVAL(1)
  477.  
  478.           NBEL   = VELCHE(/2)
  479.           NBPGAU = VELCHE(/1)
  480.           DO 1 I=1,NBEL
  481.             DO 19 J=1,NBPGAU
  482.               VELCHE(J,I) = VAL(L)
  483.               L = L + 1
  484.  19         CONTINUE
  485.  1        CONTINUE
  486.         ENDIF
  487.  8    CONTINUE
  488.  
  489. *     ------------------------------------------------------------------
  490. *     5. Ecriture et sortie
  491. *     ------------------------------------------------------------------
  492. *
  493. *     Fermeture des segments
  494. *
  495.       DO 12 IPART=1,N1
  496.         MCHAML = ICHAML(IPART)
  497.         MELVAL = IELVAL(1)
  498.         SEGDES MELVAL,MCHAML
  499.  12   CONTINUE
  500.       SEGDES MCHELM
  501.       SEGSUP TABVAL,TPART
  502. *
  503. *     Résultat
  504. *
  505.       CALL ACTOBJ('MCHAML  ',MCHELM,1)
  506.       CALL ECROBJ('MCHAML  ',MCHELM)
  507.  
  508.       END
  509.  
  510.  
  511.  
  512.  
  513.  
  514.  
  515.  

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