Télécharger la plaquette Cast3M


fr en

Télécharger smtp1.eso

Retour à la liste

Numérotation des lignes : 

smtp1

  1. C SMTP1     SOURCE    FANDEUR   22/01/03    21:15:47     11237          
  2.       SUBROUTINE SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK,
  3.      S              ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1)
  4. C-----------------------------------------------------------------------
  5. C 1) Calcul du second membre du systeme en trace de charge dans le cas
  6. C de la résolution des équations de Darcy par EFMH avec le modèle DARCY.
  7. C 2) Calcul du second membre du systeme en trace de charge et/ou d'une
  8. C contribution matricielle provenant de la convection dans le cas de la
  9. C résolution d'une équation de diffusion-convection par EFMH avec le
  10. C modèle DARCY.
  11. C-----------------------------------------------------------------------
  12. C
  13. C---------------------------
  14. C Parametres Entree/Sortie :
  15. C---------------------------
  16. C
  17. C E/  IPMAHY  : Segment contenant le pointeur vers le meleme des
  18. C               connectivites elements/faces pour les zones du MMODEL
  19. C               ou on a defini DARCY.
  20. C E/  IPFACE  : MELEME de type POI1 des faces
  21. C E/  IPRIGI  : RIGIDITE de sous type 'DARCY'
  22. C E/  THETA   : Parametre de discretisation temporelle (Theta-methode)
  23. C               pour le terme de diffusion
  24. C E/  THETAC : Parametre de discretisation temporelle (Theta-methode)
  25. C               pour le terme de convection
  26. C E/  DELTAT  : Pas de temps
  27. C E/  IPCK    : MCHAML donnant pour chaque element Ck|K|
  28. C E/  ITP     : CHPO faces des traces de charge TH au temps precedent
  29. C E/  IH      : CHPO centre des charges H au temps precedent
  30. C E/  ISOUR   : CHPO centre des sources de composante SOUR
  31. C E/  IFLUX   : CHPO face des flux de composante FLUX
  32. C E/  IORIE   : MCHAML orientation des normales
  33. C  /S MCHPOI  : CHPOINT face de composante FLUX
  34. C  /S IRI1    : RIGIDITE de sous type CONVEFMH
  35. C
  36. C----------------------
  37. C Variables en COMMON :
  38. C----------------------
  39. C
  40. C E/  IFOUR   : cf CCOPTIO
  41. C E/  NIFOUR  : cf CCOPTIO
  42. C
  43. C---------------------
  44. C Variables internes :
  45. C---------------------
  46. C
  47. C LHS1 = 1 si on doit créer une matrice (Left Hand Side); 0 sinon.
  48. C RHS1 = 1 si on doit créer un second membre (Right Hand Side); 0 sinon.
  49. C
  50. C-----------------------------------------------------------------------
  51. C
  52. C Langage : ESOPE + FORTRAN77
  53. C
  54. C Auteurs : 09/93 F.DABBENE   - Cas permanent
  55. C           09/94 X.NOUVELLON - Extension au cas transitoire
  56. C           12/95 F.DABBENE   - Extension à la convection
  57. C
  58. C-----------------------------------------------------------------------
  59.       IMPLICIT INTEGER(I-N)
  60.       IMPLICIT REAL*8 (A-H,O-Z)
  61. *
  62.  
  63. -INC PPARAM
  64. -INC CCOPTIO
  65. -INC CCHAMP
  66. -INC SMCOORD
  67. -INC SMCHPOI
  68. -INC SMCHAML
  69. -INC SMELEME
  70. -INC SMRIGID
  71. *
  72.       SEGMENT IPMAHY
  73.          INTEGER MAHYBR(NSOUS)
  74.       ENDSEGMENT
  75.       SEGMENT ICCPR
  76.          INTEGER ICPR(NNGOT)
  77.       ENDSEGMENT
  78.       SEGMENT MTRAV
  79.          REAL*8  RLIGNE(NBDDL),RCOLON(NBDDL),FLUORI(NBDDL)
  80.       ENDSEGMENT
  81. *
  82. *- Initialisations
  83. *
  84.       MCHPOI = 0
  85.       IRI1   = 0
  86.       IF (IPCK.EQ.0) THEN
  87.          IF (ISOUR.EQ.0) THEN
  88.             RHS1 = 0
  89.          ELSE
  90.             RHS1 = 1
  91.          ENDIF
  92.          IF (IFLUX.EQ.0) THEN
  93.             LHS1 = 0
  94.          ELSE
  95.             LHS1 = 1
  96.          ENDIF
  97.       ELSE
  98.          RHS1 = 1
  99.          LHS1 = 0
  100.          IF (THETAC.GT.1.D-4 .AND. IFLUX.NE.0) LHS1=1
  101.       ENDIF
  102. *
  103. *- Creation du tableau ICPR pour le MELEME POI1 IPFACE
  104. *
  105.       IK     = 0
  106.       NNGOT  = nbpts
  107.       SEGINI ICCPR
  108.       MELEME = IPFACE
  109.       SEGACT MELEME
  110.       N2 = NUM(/2)
  111.       DO 10 I2=1,N2
  112.          K = NUM(1,I2)
  113.          IF (ICPR(K).EQ.0) THEN
  114.             IK      = IK + 1
  115.             ICPR(K) = IK
  116.          ENDIF
  117.   10  CONTINUE
  118. *
  119. *- Creation du CHAMPOINT face - On laisse actif le MPOVAL du CHPOINT
  120. *
  121.       IF (RHS1.NE.0) THEN
  122.          NSOUPO = 1
  123.          NAT    = 1
  124.          SEGINI MCHPOI
  125.          IRET      = MCHPOI
  126.          MTYPOI    = 'FACE    '
  127.          IFOPOI    = IFOUR
  128.          JATTRI(1) = 2
  129.          NC        = 1
  130.          SEGINI MSOUPO
  131.          IPCHP(1)  = MSOUPO
  132.          NOHARM(1) = NIFOUR
  133.          IGEOC     = IPFACE
  134.          NOCOMP(1) = 'FLUX'
  135.          N         = N2
  136.          SEGINI MPOVAL
  137.          IPOVAL    = MPOVAL
  138.       ENDIF
  139. *
  140. *- Recuperation du nombre de zone NBMAIL du MMODEL a partir IPMAHY
  141. *
  142.       SEGACT IPMAHY
  143.       NBMAIL = MAHYBR(/1)
  144. *
  145. *- Activation du segment RIGIDITE
  146. *
  147.       MRIGID = IPRIGI
  148.       SEGACT MRIGID
  149. *
  150. *- Activation du segment MPOVAL du CHAMPOINT centre de composante SOUR
  151. *
  152.       IF (ISOUR.NE.0) THEN
  153.          MCHPO1 = ISOUR
  154.          SEGACT MCHPO1
  155.          MSOUP1 = MCHPO1.IPCHP(1)
  156.          SEGACT MSOUP1
  157.          MPOVA1 = MSOUP1.IPOVAL
  158.          SEGACT MPOVA1
  159.       ENDIF
  160. *
  161. *- Données transitoire
  162. *
  163.       IF (IPCK.NE.0) THEN
  164. *
  165. * Activation du segment MPOVAL du CHAMPOINT face de traces de charge
  166. *
  167.          MCHPO2 = ITP
  168.          SEGACT MCHPO2
  169.          MSOUP2 = MCHPO2.IPCHP(1)
  170.          SEGACT MSOUP2
  171.          MPOVA2 = MSOUP2.IPOVAL
  172.          SEGACT MPOVA2
  173. *
  174. * Activation du segment MPOVAL du CHAMPOINT centre de charges
  175. *
  176.          MCHPO3 = IH
  177.          SEGACT MCHPO3
  178.          MSOUP3 = MCHPO3.IPCHP(1)
  179.          SEGACT MSOUP3
  180.          MPOVA3 = MSOUP3.IPOVAL
  181.          SEGACT MPOVA3
  182. *
  183. * Activation du MCHAML contenant Ck|K|
  184. *
  185.          MCHELM = IPCK
  186.          SEGACT MCHELM
  187.       ENDIF
  188. *
  189. *- Activation des données pour la convection
  190. *
  191.       IF (IFLUX.NE.0) THEN
  192. *
  193. * Activation du MCHAML d'orientation des normales
  194. *
  195.          MCHEL1 = IORIE
  196.          SEGACT MCHEL1
  197. *
  198. * Activation du segment MPOVAL du CHAMPOINT flux aux faces
  199. *
  200.          MCHPO4 = IFLUX
  201.          SEGACT MCHPO4
  202.          MSOUP4 = MCHPO4.IPCHP(1)
  203.          SEGACT MSOUP4
  204.          MPOVA4 = MSOUP4.IPOVAL
  205.          SEGACT MPOVA4
  206. *
  207. * Initialisation du chapeau de rigidité
  208. *
  209.          IF (LHS1.NE.0) THEN
  210.             NRIGE  = 7
  211.             NRIGEL = NBMAIL
  212.             SEGINI RI1
  213.             IRI1 = RI1
  214.             RI1.MTYMAT = 'CONVEFMH'
  215.             RI1.IFORIG = IFORIG
  216.          ENDIF
  217.       ENDIF
  218. *
  219. *--------------------------------------------------
  220. *= BOUCLE SUR LES MAILLAGES ELEMENTAIRES,ZONE IMAIL
  221. *--------------------------------------------------
  222. *
  223.       ITELEM = 0
  224.       DO 110 IMAIL=1,NBMAIL
  225. C
  226. C- Activation de l'objet maillage ELTFA pour la zone IMAIL
  227. C
  228.          MELEME = MAHYBR(IMAIL)
  229.          IF (MELEME.EQ.0) GOTO 110
  230.          SEGACT MELEME
  231. *
  232. *- Recuperation des caracteristiques de RIGIDITE de sous type DARCY
  233. *- pour la zone IMAIL en cours de traitement
  234. *
  235.          xMATRI = IRIGEL(4,IMAIL)
  236.          SEGACT xMATRI
  237.          NBELEM = re(/3)
  238.          NELRIG = NBELEM
  239. *         XMATRI = IMATTT(1)
  240. *         SEGACT XMATRI
  241.          NBDDL  = RE(/1)
  242.          NLIGRP = NBDDL
  243.          NLIGRD = NBDDL
  244. *         SEGDES XMATRI
  245.          SEGINI MTRAV
  246. *
  247. *- Activation du segment contenant les valeurs
  248. *- du MCHAML Ck|K| pour la zone IMAIL
  249. *
  250.          IF (IPCK.NE.0) THEN
  251.             MCHAML = ICHAML(IMAIL)
  252.             SEGACT MCHAML
  253.             MELVAL = IELVAL(1)
  254.             SEGACT MELVAL
  255.          ENDIF
  256. *
  257. *- Convection
  258. *
  259.          IF (IFLUX.NE.0) THEN
  260. *
  261. * Activation du segment contenant les valeurs du MCHAML d'orientation
  262. * des normales par face pour la zone IMAIL
  263. *
  264.             MCHAM1 = MCHEL1.ICHAML(IMAIL)
  265.             SEGACT MCHAM1
  266.             MELVA1 = MCHAM1.IELVAL(1)
  267.             SEGACT MELVA1
  268. *
  269. * Initialisation du chapeau de rigidité pour la zone IMAIL
  270. *
  271.             IF (LHS1.NE.0) THEN
  272.                RI1.COERIG(IMAIL)   = 1.D0
  273.                RI1.IRIGEL(1,IMAIL) = IRIGEL(1,IMAIL)
  274.                RI1.IRIGEL(3,IMAIL) = IRIGEL(3,IMAIL)
  275.                SEGINI xMATR1
  276.                RI1.IRIGEL(4,IMAIL) = xMATR1
  277.                RI1.IRIGEL(7,IMAIL) = 2
  278.                xmatr1.symre=2
  279.             ENDIF
  280.          ENDIF
  281. *
  282. *-------------------------------------------------------
  283. * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  284. *-------------------------------------------------------
  285. *
  286.          DO 100 IEL=1,NBELEM
  287.             ITELEM = ITELEM + 1
  288. *
  289. *- Orientation du flux et calcul de la source = -div(u*Th)
  290. *
  291.             DIVUTH = 0.D0
  292.             IF (IFLUX.NE.0) THEN
  293.                DO 20 I=1,NBDDL
  294.                 IPOPTS   = ICPR(NUM(I,IEL))
  295.                 FLUORI(I)= MPOVA4.VPOCHA(IPOPTS,1)*MELVA1.VELCHE(I,IEL)
  296.                 DIVUTH   = DIVUTH - (FLUORI(I)*MPOVA2.VPOCHA(IPOPTS,1))
  297.   20           CONTINUE
  298.             ENDIF
  299. *
  300. *- Calcul de la somme des coefs pour une ligne ; une colonne
  301. *-               -1 t                      -1
  302. *-     LIGNE = RE  * DIV ; COLON = DIV * RE
  303. *-                           -1  t
  304. *- Calcul de CONSD = DIV * RE  *  DIV
  305. *
  306.             CONSD = 0.D0
  307. *            XMATRI = IMATTT(IEL)
  308. *            SEGACT XMATRI
  309.             DO 40 I=1,NBDDL
  310.                RCOLON(I) = 0.D0
  311.                RLIGNE(I) = 0.D0
  312.                DO 30 J=1,NBDDL
  313.                   RCOLON(I)   = RCOLON(I) + RE(J,I,iel)
  314.                   RLIGNE(I)   = RLIGNE(I) + RE(I,J,iel)
  315.   30           CONTINUE
  316.                CONSD = CONSD + RLIGNE(I)
  317.   40        CONTINUE
  318. *            SEGDES XMATRI
  319. *
  320. *- Recuperation de la valeur SOUR au centre de l'element
  321. *- ET ajout du terme de convection éventuelle (1-thetac)*(-div(u*Th))
  322. *
  323.             IF (IPCK.EQ.0) THEN
  324.                VSOU = MPOVA1.VPOCHA(ITELEM,1)
  325.             ELSE
  326.                IF (ISOUR.NE.0) THEN
  327.                   VSOU = MPOVA1.VPOCHA(ITELEM,1)+((1.D0-THETAC)*DIVUTH)
  328.                ELSE
  329.                   VSOU = (1.D0-THETAC) * DIVUTH
  330.                ENDIF
  331.             ENDIF
  332. *
  333. *
  334.             IF (IPCK.EQ.0) THEN
  335.                IF (RHS1.NE.0) THEN
  336.                   DO 45 IDDL=1,NBDDL
  337.                      VALFLU = RLIGNE(IDDL) * VSOU / CONSD
  338.                      IPOPTS = ICPR(NUM(IDDL,IEL))
  339.                      VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - VALFLU
  340.   45              CONTINUE
  341.                ENDIF
  342.                IF (LHS1.NE.0) THEN
  343. *                  SEGINI XMATR1
  344. *                  IMATR1.IMATTT(IEL) = XMATR1
  345.                   DO 55 J=1,NBDDL
  346.                      REMU1 = -1.D0 / CONSD
  347.                      DO 50 I=1,NBDDL
  348.                         XMATR1.RE(I,J,iel) = REMU1*RLIGNE(I)*FLUORI(J)
  349.   50                 CONTINUE
  350.   55              CONTINUE
  351. *                  SEGDES XMATR1
  352.                ENDIF
  353.             ELSE
  354.                BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT)
  355.                VH   = MPOVA3.VPOCHA(ITELEM,1)
  356.                DO 70 IDDL=1,NBDDL
  357.                   SMFLU  = RLIGNE(IDDL) * VSOU * BETA/CONSD
  358.                   SMP    = RLIGNE(IDDL) * VH   * (1.D0-BETA)
  359.                   STP    = 0.D0
  360.                   DO 60 J=1,NBDDL
  361.                      VTH   = MPOVA2.VPOCHA(ICPR(NUM(J,IEL)),1)
  362.                      STP   = STP + RCOLON(J) * VTH
  363.   60              CONTINUE
  364.                   SMTP   = RLIGNE(IDDL) * STP * (1.D0-THETA)*BETA/CONSD
  365.                   SMDDL  = SMFLU + SMP + SMTP
  366.                   IPOPTS = ICPR(NUM(IDDL,IEL))
  367.                   VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - SMDDL
  368.   70           CONTINUE
  369.                IF (LHS1.NE.0) THEN
  370. *                  SEGINI XMATR1
  371. *                  IMATR1.IMATTT(IEL) = XMATR1
  372.                   REMU1 = (-1.D0) * THETAC * BETA / CONSD
  373.                   DO 90 J=1,NBDDL
  374.                      DO 80 I=1,NBDDL
  375.                         XMATR1.RE(I,J,iel)=REMU1*RLIGNE(I) * FLUORI(J)
  376.   80                 CONTINUE
  377.   90              CONTINUE
  378. *                  SEGDES XMATR1
  379.                ENDIF
  380.             ENDIF
  381.  100     CONTINUE
  382.          SEGDES xMATRI
  383.          SEGSUP MTRAV
  384.  110  CONTINUE
  385. *
  386. *- Desactivation des segments
  387. *
  388.       IF (IFLUX.NE.0) THEN
  389.          IF (LHS1.NE.0) SEGDES RI1
  390.       ENDIF
  391.       SEGDES IPMAHY
  392.       SEGDES MRIGID
  393.       SEGSUP ICCPR
  394. *
  395.       END
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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