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    CB215821  24/04/12    21:15:03     11897          
  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.  
  528.  

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