Télécharger la plaquette Cast3M


fr en

Télécharger smtp.eso

Retour à la liste

Numérotation des lignes : 

smtp

  1. C SMTP      SOURCE    CB215821  24/04/12    21:17:15     11897          
  2.       SUBROUTINE SMTP
  3. C-----------------------------------------------------------------------
  4. C 1) Calcul du second membre du systeme en trace de charge dans le cas
  5. C de la résolution des équations de Darcy par EFMH avec le modèle DARCY.
  6. C 2) Calcul du second membre du systeme en trace de charge et/ou d'une
  7. C contribution matricielle provenant de la convection dans le cas de la
  8. C résolution d'une équation de diffusion-convection par EFMH avec le
  9. C modèle DARCY.
  10. C-----------------------------------------------------------------------
  11. C
  12. C---------------------------
  13. C Phrase d'appel (GIBIANE) :
  14. C---------------------------
  15. C
  16. C  (RIG2) (CHP3) = 'SMTP' MMODEL TABLE1 RIG1 (TABLE2) (CHP1) (CHP2) ;
  17. C
  18. C------------------------
  19. C Operandes et resultat :
  20. C------------------------
  21. C
  22. C  MMODEL  : Objet modele décrivant la formulation.
  23. C  TABLE1  : TABLE DOMAINE contenant les maillages et les connectivités.
  24. C  RIG1    : Matrices hybrides elementaires de DARCY crees par MHYB.
  25. C  TABLE2  : TABLE DARCY_TRANSITOIRE contenant les infos transitoires.
  26. C  CHP1    : CHPO centre de composante SOUR, source moyenne par élément.
  27. C  CHP2    : CHPO face de composante FLUX, flux de vitesse convective.
  28. C  CHP3    : CHPO face de composante FLUX, second membre du pb en TH.
  29. C  RIG2    : RIGIDITE, contribution du terme convectif au LHS en TH.
  30. C
  31. C-----------------------------------------------------------------------
  32. C
  33. C Langage : ESOPE + FORTRAN77
  34. C
  35. C Auteurs : 09/93 F.DABBENE   - Cas permanent
  36. C           09/94 X.NOUVELLON - Extension au transitoire
  37. C           12/95 F.DABBENE   - Extension à la convection
  38. C
  39. C-----------------------------------------------------------------------
  40.       IMPLICIT INTEGER(I-N)
  41.       IMPLICIT REAL*8 (A-H,O-Z)
  42. *
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC SMCHAML
  47. -INC SMCHPOI
  48. -INC SMELEME
  49. -INC SMMODEL
  50. -INC SMRIGID
  51. -INC SMTABLE
  52. *
  53.       SEGMENT IPMAHY
  54.          INTEGER MAHYBR(NSOUS)
  55.       ENDSEGMENT
  56. *
  57.       LOGICAL LOGRE,LOGIN
  58.       CHARACTER*4 NOMTOT(1)
  59.       CHARACTER*6 CHAR6
  60.       CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  61. *
  62. * Initialisations
  63. *
  64.       ISOU1  = 0
  65.       IHN1   = 0
  66.       ITPN1  = 0
  67.       IVALIN = 0
  68.       XVALIN = 0.D0
  69.       LOGIN  = .TRUE.
  70.       IOBIN  = 0
  71.       TAPIND = 'MOT     '
  72. *
  73. * Lecture du MMODEL
  74. *
  75.       CALL LIROBJ('MMODEL  ',IPMODE,1,IRET)
  76.       CALL ACTOBJ('MMODEL  ',IPMODE,1)
  77.       IF (IERR.NE.0) RETURN
  78.       MMODEL = IPMODE
  79. *
  80. * Lecture de la TABLE domaine
  81. *
  82.       IPTABL = 0
  83.       CALL LEKMOD(IPMODE,IPTABL,IRET)
  84.       IF (IERR.NE.0) RETURN
  85.       CHARIN = 'MAILLAGE'
  86.       TYPOBJ = 'MAILLAGE'
  87.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  88.       IF (IERR.NE.0) RETURN
  89.       IPGEOM = IOBRE
  90.       CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  91.       IF (IERR.NE.0) RETURN
  92.       IELTFA = IOBRE
  93.       CALL LEKTAB(IPTABL,'CENTRE',IOBRE)
  94.       IF (IERR.NE.0) RETURN
  95.       ICENTR = IOBRE
  96.       CALL LEKTAB(IPTABL,'FACE',IOBRE)
  97.       IF (IERR.NE.0) RETURN
  98.       IPFACE = IOBRE
  99. *
  100. * Lecture de RIGIDITE
  101. *
  102.       CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  103.       IF (IERR.NE.0) RETURN
  104.       MRIGID = IPRIGI
  105. *
  106. * Lecture éventuelle de la TABLE Darcy_transitoire
  107. *
  108.       IPTAB1 = 0
  109.       CALL LIRTAB('DARCY_TRANSITOIRE',IPTAB1,0,IRET)
  110.       IF (IERR.NE.0) RETURN
  111.       IF (IRET.EQ.1) THEN
  112.          THEMIN = -1.D-12
  113.          THEMAX = 1.D0 - THEMIN
  114.          TYPOBJ = 'FLOTTANT'
  115.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  116.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  117.          IF (IERR.NE.0) RETURN
  118.          THETA  = XVALRE
  119.          IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  120.             REAERR(1)   = REAL(THETA)
  121. C            REAERR(2)   = REAL(0.D0)
  122. C            REAERR(3)   = REAL(1.D0)
  123.             MOTERR(1:8) = '  THETA '
  124.             CALL ERREUR(42)
  125.             RETURN
  126.          ENDIF
  127.          IF (THETA.LT.0.D0) THETA=0.D0
  128.          IF (THETA.GT.1.D0) THETA=1.D0
  129.          TYPOBJ = '        '
  130.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'THETA_CONVECTION',
  131.      S                                              LOGIN,IOBIN,
  132.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  133.          IF (IERR.NE.0) RETURN
  134.          IF (TYPOBJ.EQ.'FLOTTANT') THEN
  135.            THETAC  = XVALRE
  136.            IF ((THETAC.LT.THEMIN).OR.(THETAC.GT.THEMAX)) THEN
  137.               REAERR(1)   = REAL(THETAC)
  138.               REAERR(2)   = REAL(0.D0)
  139.               REAERR(3)   = REAL(1.D0)
  140.               MOTERR(1:8) = 'TETACONV'
  141.               CALL ERREUR(42)
  142.               RETURN
  143.            ENDIF
  144.            IF (THETAC.LT.0.D0) THETA=0.D0
  145.            IF (THETAC.GT.1.D0) THETA=1.D0
  146.          ELSE
  147.            THETAC  = THETA
  148.          ENDIF
  149.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  150.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  151.          IF (IERR.NE.0) RETURN
  152.          DELTAT = XVALRE
  153.          IF (DELTAT.LE.0.D0) THEN
  154.             REAERR(1)   = REAL(DELTAT)
  155.             REAERR(2)   = REAL(0.D0)
  156.             MOTERR(1:8) = ' DELTAT '
  157.             CALL ERREUR(41)
  158.             RETURN
  159.          ENDIF
  160.          TYPOBJ = 'MCHAML  '
  161.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  162.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  163.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  164.          IF (IERR.NE.0) RETURN
  165.          IPCK = IOBRE
  166.          TYPOBJ = 'CHPOINT '
  167.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN,
  168.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  169.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  170.          IF (IERR.NE.0) RETURN
  171.          ITP = IOBRE
  172.          CALL ACCTAB(IPTAB1,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN,
  173.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  174.          IF (IERR.NE.0) RETURN
  175.          IH = IOBRE
  176.       ELSE
  177.          THETA  = 0.D0
  178.          DELTAT = 0.D0
  179.          IPCK   = 0
  180.          ITP    = 0
  181.          IH     = 0
  182.       ENDIF
  183. *
  184. * Lecture des champoints facultatifs SOUR et FLUX
  185. *     prise en compte des termes sources
  186. *     prise en compte de la convection
  187. *
  188.       ISOUR = 0
  189.       CALL LIROBJ('CHPOINT ',ISOUR,0,IRET)
  190.       IF(IRET .EQ. 1) CALL ACTOBJ('CHPOINT ',ISOUR,1)
  191.       IF (IERR.NE.0) RETURN
  192. *
  193.       IFLUX = 0
  194.       CALL LIROBJ('CHPOINT ',IFLUX,0,IRET)
  195.       IF(IRET .EQ. 1) CALL ACTOBJ('CHPOINT ',IFLUX,1)
  196.       IF (IERR.NE.0) RETURN
  197. *
  198. * Tri des deux champoints facultatifs
  199. *
  200.       IFAC = ISOUR + IFLUX
  201.       IF (IFAC.EQ.0.AND.IPTAB1.EQ.0) THEN
  202.          MOTERR(1:8) = 'CHAMPOIN'
  203.          CALL ERREUR(37)
  204.          RETURN
  205.       ENDIF
  206.       IF (IFAC.NE.0) THEN
  207.          NBCOMP = 0
  208.          INDIC  = 0
  209.          CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
  210.          IF (NOMTOT(1).NE.'SOUR'.AND.NOMTOT(1).NE.'FLUX') THEN
  211.             MOTERR(1:4) = NOMTOT(1)
  212.             CALL ERREUR(197)
  213.             RETURN
  214.          ELSE IF (NOMTOT(1).EQ.'FLUX') THEN
  215.             IFAC  = IFLUX
  216.             IFLUX = ISOUR
  217.             ISOUR = IFAC
  218.          ENDIF
  219.       ENDIF
  220. *
  221. * Dans le cas de la convection, récup. MCHAML orientation normale
  222. *
  223.       IF (IFLUX.NE.0) THEN
  224.          CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  225.          IF (IERR.NE.0) RETURN
  226.       ENDIF
  227. *
  228. * Test du CHAMPOINT de composante SOUR de spg CENTRE
  229. *
  230.       IF (ISOUR.NE.0) THEN
  231.          NBCOMP    = 1
  232.          INDIC     = 1
  233.          NOMTOT(1) = 'SOUR'
  234.          CALL QUEPOI(ISOUR,ICENTR,INDIC,NBCOMP,NOMTOT)
  235.          IF (IERR.NE.0) RETURN
  236.       ENDIF
  237. *
  238. * Test du CHAMPOINT de composante FLUX de spg FACE
  239. *
  240.       IF (IFLUX.NE.0) THEN
  241.          NBCOMP    = 1
  242.          INDIC     = 1
  243.          NOMTOT(1) = 'FLUX'
  244.          CALL QUEPOI(IFLUX,IPFACE,INDIC,NBCOMP,NOMTOT)
  245.          IF (IERR.NE.0) RETURN
  246.       ENDIF
  247. *
  248. * Test du CHAMPOINT de composante TH de spg FACES
  249. *
  250.       IF (ITP.NE.0) THEN
  251.          NBCOMP    = 1
  252.          INDIC     = 1
  253.          NOMTOT(1) = 'TH'
  254.          CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT)
  255.          IF (IERR.NE.0) RETURN
  256.       ENDIF
  257. *
  258. * Test du CHAMPOINT de composante H de spg CENTRE
  259. *
  260.       IF (IH.NE.0) THEN
  261.          NBCOMP = 1
  262.          INDIC  = 1
  263.          NOMTOT(1) = 'H'
  264.          CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT)
  265.          IF (IERR.NE.0) RETURN
  266.       ENDIF
  267. *
  268. * Test de la formulation
  269. *
  270.       SEGACT MMODEL
  271.       NSOUS  = KMODEL(/1)
  272.       SEGINI IPMAHY
  273.       IDARCY = 0
  274.       IPT1=IPGEOM
  275.       IPT2=IPGEOM
  276.       DO 10 ISOUS=1,NSOUS
  277.          IF(NSOUS.GT.1)THEN
  278.               SEGACT IPT2
  279.               IPT1=IPT2.LISOUS(ISOUS)
  280.          ENDIF
  281.          IMODEL = KMODEL(ISOUS)
  282.          SEGACT IMODEL
  283.          LETYPE = FORMOD(1)
  284.          IF (LETYPE.EQ.'DARCY') THEN
  285.             IDARCY = IDARCY + 1
  286.             MAHYBR(ISOUS) = IPT1
  287.          ENDIF
  288.  10   CONTINUE
  289.       IF (IDARCY.EQ.0) THEN
  290.          MOTERR = LETYPE
  291.          CALL ERREUR(193)
  292.          GOTO 100
  293.       ENDIF
  294. *
  295. * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  296. *
  297.       MELEME = IELTFA
  298.       SEGACT MELEME
  299.       LZONES = LISOUS(/1)
  300.       IF (LZONES.EQ.0) LZONES = 1
  301.       IPT1   = IPGEOM
  302.       SEGACT IPT1
  303.       DO 30 ISOUS=1,NSOUS
  304.          IMAMEL = MAHYBR(ISOUS)
  305.          IF (IMAMEL.NE.0) THEN
  306.             DO 20 ISZ=1,LZONES
  307.                IF (LZONES.EQ.1) THEN
  308.                   IPT2 = IPT1
  309.                   IPT3 = MELEME
  310.                ELSE
  311.                   IPT2 = IPT1.LISOUS(ISZ)
  312.                   IPT3 = LISOUS(ISZ)
  313.                ENDIF
  314.                IF (IPT2.EQ.IMAMEL) THEN
  315.                   MAHYBR(ISOUS) = IPT3
  316.                   GOTO 30
  317.                ENDIF
  318.  20         CONTINUE
  319.             IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  320.                INTERR(1) = ISOUS
  321.                CALL ERREUR(698)
  322.                GOTO 100
  323.             ENDIF
  324.          ENDIF
  325.  30   CONTINUE
  326. *
  327. * Test du sous-type de la matrice de rigiditée récupérée
  328. *
  329.       SEGACT MRIGID
  330.       LETYPE = MTYMAT
  331.       IF (LETYPE.NE.'DARCY') THEN
  332.          MOTERR(1:8)  = 'RIGIDITE'
  333.          MOTERR(9:16) = 'DARCY   '
  334.          CALL ERREUR(79)
  335.          SEGDES MRIGID
  336.          GOTO 100
  337.       ENDIF
  338. *
  339. * Controle des pointeurs de MELEME de la rigidité
  340. *
  341.       DO 40 ISOUS=1,NSOUS
  342.          IMAMEL = MAHYBR(ISOUS)
  343.          IF (IMAMEL.NE.0) THEN
  344.             IPTEST = IRIGEL(1,ISOUS)
  345.             IF (MAHYBR(ISOUS).NE.IRIGEL(1,ISOUS)) THEN
  346.                MOTERR(1:8)  = 'DARCY   '
  347.                MOTERR(9:16) = 'ELTFA   '
  348.                INTERR(1)    = ISOUS
  349.                CALL ERREUR(698)
  350.                SEGDES MRIGID
  351.                GOTO 100
  352.             ENDIF
  353.          ENDIF
  354.  40   CONTINUE
  355.       SEGDES MRIGID
  356. *
  357. * Vérification du support du MCHAML % au MMODEL
  358. * Controle du nombre de composantes reelles
  359. *
  360.       IF (IPCK.NE.0) THEN
  361.          ISUP  = 2
  362.          ICOND = 1
  363.          CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  364.          IF (IRET.GT.1) GOTO 100
  365.          MCHELM = IPCK
  366.          SEGACT MCHELM
  367.          DO 50 ISOUS=1,NSOUS
  368.             IF (MAHYBR(ISOUS).NE.0) THEN
  369.                MCHAML = ICHAML(ISOUS)
  370.                SEGACT MCHAML
  371.                N2     = TYPCHE(/2)
  372.                IF (N2.GT.1) THEN
  373.                   CALL ERREUR(320)
  374.                   GOTO 100
  375.                ENDIF
  376.                CHAR6  = TYPCHE(1)(1:6)
  377.                IF (CHAR6.NE.'REAL*8') THEN
  378.                   MOTERR(1:8) = NOMCHE(1)
  379.                   CALL ERREUR(679)
  380.                   GOTO 100
  381.                ENDIF
  382.             ENDIF
  383.  50      CONTINUE
  384.       ENDIF
  385. *
  386. * Construction de MCHPOI
  387. *
  388.       SEGDES IPMAHY
  389.       CALL SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK,
  390.      S           ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1)
  391. *
  392.       IF (MCHPOI.NE.0) THEN
  393.         CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  394.         CALL ECROBJ('CHPOINT ',MCHPOI)
  395.       ENDIF
  396.       IF (IRI1.NE.0)   CALL ECROBJ('RIGIDITE',IRI1)
  397. *
  398. * Ménage
  399. *
  400.  100  CONTINUE
  401.       SEGSUP IPMAHY
  402.       END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  

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