Télécharger la plaquette Cast3M


fr en

Télécharger hybp.eso

Retour à la liste

Numérotation des lignes : 

hybp

  1. C HYBP      SOURCE    CB215821  24/04/12    21:16:18     11897          
  2.       SUBROUTINE HYBP
  3. C-----------------------------------------------------------------------
  4. C Cet operateur permet de calculer la charge en fonction des traces de
  5. C charge dans le cadre d'une formulation variationnelle mixte hybride
  6. C pour les equations de DARCY.
  7. C
  8. C        (n)               (n-1)         -1 t    -1
  9. C       H  = (1 - BETA) * H   + (DIV * M1  * DIV)  * BETA *
  10. C                      -1             (n)              (n-1)
  11. C            ( DIV * M1  * (THETA * TH + (1-THETA) * TH  - M * F) + SOUR  )
  12. C
  13. C---------------------------
  14. C Phrase d'appel (GIBIANE) :
  15. C---------------------------
  16. C
  17. C  CHPO1 = HYBP MMODEL MRIGI1 CHPO2 (TABLE2) (CHPO3) (MRIGI2 CHPO4);
  18. C
  19. C------------------------
  20. C Operandes et resultat :
  21. C------------------------
  22. C
  23. C  CHPO1   : CHPO centre de nom de composante H.
  24. C  MRIGI1  : RIGIDITE masses hybrides élémentaires de sous-type DARCY.
  25. C  MRIGI2    : Matrices masses hybrides elementaires de sous-type MASSE.
  26. C  MMODEL  : MODELE spécifiant la formulation.
  27. C  CHPO2   : CHPO face de nom de composante TH.
  28. C  TABLE2  : TABLE DARCY_TRANSITOIRE a inclure en transitoire; contient
  29. C            Theta, le pas de temps, Ck|K|, TH(n-1) et H(n-1).
  30. C  CHPO3   : CHPO centre de nom de composante SOUR.
  31. C  CHPO4    : CHPO centre de composante FX,FY(,FZ),
  32. C                                  densite de force moyenne par élément.
  33. C
  34. C-----------------------------------------------------------------------
  35. C
  36. C Langage : ESOPE + FORTRAN77
  37. C
  38. C Auteurs : 09/93 F.DABBENE - Cas permanent
  39. C           09/94 X.NOUVELLON - Extension au cas transitoire
  40. C           02/96 L.V.BENET - Prise en compte de forces de volume
  41. C
  42. C-----------------------------------------------------------------------
  43.       IMPLICIT INTEGER(I-N)
  44.       IMPLICIT REAL*8 (A-H,O-Z)
  45. *
  46.  
  47. -INC PPARAM
  48. -INC CCOPTIO
  49. -INC SMCHAML
  50. -INC SMCHPOI
  51. -INC SMELEME
  52. -INC SMMODEL
  53. -INC SMRIGID
  54. -INC SMTABLE
  55. *
  56.       SEGMENT IPMAHY
  57.          INTEGER MAHYBR(NSOUS)
  58.       ENDSEGMENT
  59. *
  60.       LOGICAL LOGRE,LOGIN
  61.       CHARACTER*4 NOMTOT(1)
  62.       CHARACTER*4 NOMTO3(3)
  63.       CHARACTER*6 CHAR6
  64.       CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  65. *
  66. * Initialisations
  67. *
  68.       IVALIN = 0
  69.       XVALIN = 0.D0
  70.       LOGIN  = .TRUE.
  71.       IOBIN  = 0
  72.       TAPIND = 'MOT     '
  73. *
  74. * Lecture du MMODEL et test de la formulation
  75. *
  76.       CALL LIROBJ('MMODEL ',IPMODE,1,IRET)
  77.       CALL ACTOBJ('MMODEL ',IPMODE,1)
  78.       IF (IERR.NE.0) RETURN
  79.       MMODEL = IPMODE
  80. *
  81. * Lecture de la TABLE domaine
  82. *
  83.       IPTABL = 0
  84. C     CALL LIRTAB('DOMAINE',IPTABL,1,IRET)
  85.       CALL LEKMOD(MMODEL,IPTABL,IRET)
  86.       IF (IERR.NE.0) RETURN
  87.       CHARIN = 'MAILLAGE'
  88.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  89.       IF (IERR.NE.0) RETURN
  90.       IPGEOM = IOBRE
  91.       CHARIN = 'ELTFA'
  92.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  93.       IF (IERR.NE.0) RETURN
  94.       IELTFA = IOBRE
  95.       CHARIN = 'CENTRE'
  96.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  97.       IF (IERR.NE.0) RETURN
  98.       ICENTR = IOBRE
  99.       CHARIN = 'FACE'
  100.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  101.       IF (IERR.NE.0) RETURN
  102.       IPFACE = IOBRE
  103. *
  104. * récup. MCHAML orientation normale
  105. *
  106.       CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  107.       IF (IERR.NE.0) RETURN
  108. *
  109. * récup. CHPO orientation normale
  110. *
  111.       CALL LEKTAB(IPTABL,'XXNORMAF',INORM)
  112.       IF (IERR.NE.0) RETURN
  113. *
  114. * récup. CHPO surface des faces
  115. *
  116.       CALL LEKTAB(IPTABL,'XXSURFAC',ISURF)
  117.       IF (IERR.NE.0) RETURN
  118. *
  119. * Lecture de RIGIDITE obligatoire
  120. *
  121.       CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET)
  122.       IF (IERR.NE.0) RETURN
  123.       RI1 = IPRIGI
  124. *
  125. * Lecture du CHAMPOINT obligatoire
  126. *
  127.       CALL LIROBJ('CHPOINT ',ITPN,1,IRET)
  128.       CALL ACTOBJ('CHPOINT ',ITPN,1)
  129.       IF (IERR.NE.0) RETURN
  130. *
  131. * Lecture eventuelle du 1er CHAMPOINT optionnel
  132. *
  133.       ICHP1 = 0
  134.       CALL LIROBJ('CHPOINT',ICHP1,0,IRET1)
  135.       IF(IRET1 .EQ.1) CALL ACTOBJ('CHPOINT',ICHP1,1)
  136.       IF (IERR.NE.0) RETURN
  137. *
  138. * Lecture eventuelle du 2me CHAMPOINT optionnel
  139. *
  140.       ICHP2 = 0
  141.       CALL LIROBJ('CHPOINT',ICHP2,0,IRET2)
  142.       IF(IRET2 .EQ.1) CALL ACTOBJ('CHPOINT',ICHP2,1)
  143.       IF (IERR.NE.0) RETURN
  144. *
  145. * Lecture de la RIGIDITE complementaire
  146. *
  147.        IPRIGI = 0
  148.        CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR)
  149.       IF (IERR.NE.0) RETURN
  150.        RI2 = IPRIGI
  151. *
  152. * Lecture éventuelle de la TABLE Darcy_transitoire
  153. *
  154.       IPTABL = 0
  155.       CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET)
  156.       IF (IERR.NE.0) RETURN
  157.       IF (IRET.EQ.1) THEN
  158.          TYPOBJ = 'FLOTTANT'
  159.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  160.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  161.          IF (IERR.NE.0) RETURN
  162.          THETA  = XVALRE
  163.          THEMIN = -1.D-12
  164.          THEMAX = 1.D0 - THEMIN
  165.          IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  166.             REAERR(1)   = REAL(THETA)
  167.             REAERR(2)   = REAL(0.D0)
  168.             REAERR(3)   = REAL(1.D0)
  169.             MOTERR(1:8) = '  THETA '
  170.             CALL ERREUR(42)
  171.             RETURN
  172.          ENDIF
  173.          IF (THETA.LT.0.D0) THETA=0.D0
  174.          IF (THETA.GT.1.D0) THETA=1.D0
  175.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  176.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  177.          IF (IERR.NE.0) RETURN
  178.          DELTAT = XVALRE
  179.          IF (DELTAT.LT.0.D0) THEN
  180.             REAERR(1)   = REAL(DELTAT)
  181.             REAERR(2)   = REAL(0.D0)
  182.             MOTERR(1:8) = ' DELTAT '
  183.             CALL ERREUR(41)
  184.             RETURN
  185.          ENDIF
  186.          TYPOBJ = 'MCHAML  '
  187.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  188.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  189.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  190.          IF (IERR.NE.0) RETURN
  191.          IPCK = IOBRE
  192.          TYPOBJ = 'CHPOINT '
  193.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN,
  194.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  195.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  196.          IF (IERR.NE.0) RETURN
  197.          ITP = IOBRE
  198.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN,
  199.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  200.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  201.          IF (IERR.NE.0) RETURN
  202.          IH = IOBRE
  203.       ELSE
  204.          THETA  = 0.D0
  205.          DELTAT = 0.D0
  206.          IPCK   = 0
  207.          ITP    = 0
  208.          IH     = 0
  209.       ENDIF
  210. *
  211. * Test CHAMPOINT centre H
  212. *
  213.       IF (IH.NE.0) THEN
  214.          NBCOMP    = 1
  215.          INDIC     = 1
  216.          NOMTOT(1) = 'H'
  217.          CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT)
  218.          IF (IERR.NE.0) RETURN
  219.       ENDIF
  220. *                          (n-1)
  221. * Test du CHAMPOINT face TH
  222. *
  223.       IF (ITP.NE.0) THEN
  224.          NBCOMP    = 1
  225.          INDIC     = 1
  226.          NOMTOT(1) = 'TH'
  227.          CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT)
  228.          IF (IERR.NE.0) RETURN
  229.       ENDIF
  230. *
  231. * Test de la formulation contenu dans MMODEL
  232. *
  233.       SEGACT MMODEL
  234.       NSOUS  = KMODEL(/1)
  235.       SEGINI IPMAHY
  236.       IDARCY = 0
  237.       IPT1=IPGEOM
  238.       IPT2=IPGEOM
  239.       DO 10 ISOUS=1,NSOUS
  240.          IF(NSOUS.GT.1)THEN
  241.               SEGACT IPT2
  242.               IPT1= IPT2.LISOUS(ISOUS)
  243.          ENDIF
  244.          IMODEL = KMODEL(ISOUS)
  245.          SEGACT IMODEL
  246.          LETYPE = FORMOD(1)
  247.          IF (LETYPE.EQ.'DARCY') THEN
  248.             IDARCY = IDARCY + 1
  249.             MAHYBR(ISOUS) = IPT1
  250.          ENDIF
  251.  10   CONTINUE
  252.       IF (IDARCY.EQ.0) THEN
  253.          MOTERR = LETYPE
  254.          CALL ERREUR(193)
  255.          GOTO 100
  256.       ENDIF
  257. *
  258. * Tri et test des CHAMPOINTs ARGUMENTS
  259. *
  260.       ISOUR=0
  261.       IFORC=0
  262.       IF(IRET1.EQ.1)THEN
  263.           NBCOMP    = 0
  264.           INDIC     = 1
  265.           NOMTO3(1) = ' '
  266.           NOMTO3(2) = ' '
  267.           NOMTO3(3) = ' '
  268.           CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTO3)
  269.           IF (IERR.NE.0) RETURN
  270.           IF(NOMTO3(1).EQ.'SOUR')THEN
  271.              ISOUR=ICHP1
  272.           ELSE IF(NOMTO3(1).EQ.'FX')THEN
  273.              IFORC=ICHP1
  274.           ENDIF
  275.           IF(IRET2.EQ.1)THEN
  276.              NBCOMP    = 0
  277.              INDIC     = 1
  278.              NOMTO3(1) = ' '
  279.              NOMTO3(2) = ' '
  280.              NOMTO3(3) = ' '
  281.              CALL QUEPOI(ICHP2,ICENTR,INDIC,NBCOMP,NOMTO3)
  282.              IF (IERR.NE.0) RETURN
  283.              IF(NOMTO3(1).EQ.'SOUR')THEN
  284.                 ISOUR=ICHP2
  285.              ELSE IF(NOMTO3(1).EQ.'FX')THEN
  286.                 IFORC=ICHP2
  287.              ENDIF
  288.           ENDIF
  289.       ENDIF
  290. *
  291. * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  292. *
  293.       MELEME = IELTFA
  294.       SEGACT MELEME
  295.       LZONES = LISOUS(/1)
  296.       IF (LZONES.EQ.0) LZONES = 1
  297.       IPT4   = IPGEOM
  298.       SEGACT IPT4
  299.       DO 30 ISOUS=1,NSOUS
  300.          IMAMEL = MAHYBR(ISOUS)
  301.          IF (IMAMEL.NE.0) THEN
  302.             DO 20 ISZ=1,LZONES
  303.                IF (LZONES.EQ.1) THEN
  304.                   IPT2 = IPT4
  305.                   IPT3 = MELEME
  306.                ELSE
  307.                   IPT2 = IPT4.LISOUS(ISZ)
  308.                   IPT3 = LISOUS(ISZ)
  309.                ENDIF
  310.                IF (IPT2.EQ.IMAMEL) THEN
  311.                   MAHYBR(ISOUS) = IPT3
  312.                   GOTO 30
  313.                ENDIF
  314.  20         CONTINUE
  315.             IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  316.                MOTERR(1:8)  = ' MODELE '
  317.                MOTERR(9:16) = ' ELTFA  '
  318.                INTERR(1)    = ISOUS
  319.                CALL ERREUR(698)
  320.                GOTO 100
  321.             ENDIF
  322.          ENDIF
  323.  30   CONTINUE
  324. *
  325. * Test et tri du sous-type des matrices de rigiditée récupérées
  326. *
  327.       IPFORC=0
  328.       IPDARC=0
  329.       SEGACT RI1
  330.       LETYPE =  RI1.MTYMAT
  331.       IF(LETYPE.EQ.'DARCY')THEN
  332.         IPDARC=RI1
  333.       ELSE IF(LETYPE.EQ.'MASSE')THEN
  334.         IPFORC=RI1
  335.       ENDIF
  336.       IF(RI2.NE.0)THEN
  337.         SEGACT RI2
  338.         LETYPE =  RI2.MTYMAT
  339.         IF(LETYPE.EQ.'MASSE')THEN
  340.            IPFORC=RI2
  341.         ELSE IF(LETYPE.EQ.'DARCY')THEN
  342.            IPDARC=RI2
  343.         ENDIF
  344.         IF(IPFORC.EQ.0)THEN
  345.           MOTERR(1:8)  = 'RIGIDITE'
  346.           MOTERR(9:16) = 'MASSE'
  347.           CALL ERREUR(79)
  348.           SEGDES RI1
  349.           SEGDES RI2
  350.           GOTO 100
  351.         ELSE IF(IFORC.EQ.0)THEN
  352.               MOTERR(1:8)  = 'CHPOINT'
  353.               MOTERR(9:16) = 'FORCE'
  354.               CALL ERREUR(79)
  355.               SEGDES RI1
  356.               SEGDES RI2
  357.               GOTO 100
  358.         ENDIF
  359.       ENDIF
  360.       IF(IPDARC.EQ.0)THEN
  361.         MOTERR(1:8)  = 'RIGIDITE'
  362.         MOTERR(9:16) = 'DARCY'
  363.         CALL ERREUR(79)
  364.         SEGDES RI1
  365.         IF(IPFORC.NE.0)SEGDES RI2
  366.         GOTO 100
  367.       ENDIF
  368. *
  369. * Controle des pointeurs de MELEME de la rigidité
  370. *
  371.       DO 40 ISOUS=1,NSOUS
  372.          IMAMEL = MAHYBR(ISOUS)
  373.          IF (IMAMEL.NE.0) THEN
  374.             RI1=IPDARC
  375.             IPTEST = RI1.IRIGEL(1,ISOUS)
  376.             CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  377.             IF (.NOT. LOGRE) THEN
  378.                MOTERR(1:8)  = 'DARCY'
  379.                MOTERR(9:16) = ' ELTFA  '
  380.                INTERR(1)    = ISOUS
  381.                CALL ERREUR(698)
  382.                SEGDES RI1
  383.                GOTO 100
  384.             ENDIF
  385.             IF(IPFORC.NE.0)THEN
  386.               RI2=IPFORC
  387.               IPTEST = RI2.IRIGEL(1,ISOUS)
  388.                CALL CMPMAI(IMAMEL,IPTEST,LOGRE)
  389.                IF (.NOT. LOGRE) THEN
  390.                  MOTERR(1:8)  = 'MASSE'
  391.                  MOTERR(9:16) = ' ELTFA  '
  392.                  INTERR(1)    = ISOUS
  393.                  CALL ERREUR(698)
  394.                  SEGDES RI1
  395.                  SEGDES RI2
  396.                  GOTO 100
  397.                ENDIF
  398.             ENDIF
  399.          ENDIF
  400.  40   CONTINUE
  401.       SEGDES RI1
  402.       IF(IPFORC.NE.0)SEGDES RI2
  403. *
  404. * Vérification du support du MCHAML Ck|K| % au MMODEL
  405. * Controle du nombre de composantes reelles
  406. *
  407.       IF (IPCK.NE.0) THEN
  408.          IF(IPFORC.NE.0)CALL ERREUR(34)
  409.          ISUP  = 2
  410.          ICOND = 1
  411.          CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  412.          IF (IRET.GT.1) GOTO 100
  413.          MCHELM = IPCK
  414.          SEGACT MCHELM
  415.          DO 50 ISOUS=1,NSOUS
  416.             IF (MAHYBR(ISOUS).NE.0) THEN
  417.                MCHAML = ICHAML(ISOUS)
  418.                SEGACT MCHAML
  419.                N2     = TYPCHE(/2)
  420.                IF (N2.GT.1) THEN
  421.                   CALL ERREUR(320)
  422.                   GOTO 100
  423.                ENDIF
  424.                CHAR6  = TYPCHE(1)(1:6)
  425.                IF (CHAR6.NE.'REAL*8') THEN
  426.                   MOTERR(1:8) = NOMCHE(1)
  427.                   CALL ERREUR(679)
  428.                   GOTO 100
  429.                ENDIF
  430.             ENDIF
  431.  50      CONTINUE
  432.       ENDIF
  433. *                          (n)
  434. * Test du CHAMPOINT face TH
  435. * On ne teste que la présence d'une composante de nom TH
  436. *
  437.       MCHPO1 = ITPN
  438.       SEGACT MCHPO1
  439.       NSOUP1 = MCHPO1.IPCHP(/1)
  440.       DO 70 ISOUS=1,NSOUP1
  441.          MSOUPO = MCHPO1.IPCHP(ISOUS)
  442.          SEGACT MSOUPO
  443.          NC = NOCOMP(/2)
  444.          DO 60 J=1,NC
  445.             IF (NOCOMP(J).EQ.'TH  ') GOTO 80
  446.  60      CONTINUE
  447.  70   CONTINUE
  448.       MOTERR(1:4) = ' TH '
  449.       CALL ERREUR(181)
  450.       GOTO 100
  451.  80   CONTINUE
  452. *
  453. * Construction du CHPO des charges
  454. *
  455.       SEGDES IPMAHY
  456.       CALL HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE,
  457.      S  ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,ICHP1)
  458.  
  459.       IF (ICHP1.NE.0) THEN
  460.         CALL ACTOBJ('CHPOINT ',ICHP1,1)
  461.         CALL ECROBJ('CHPOINT ',ICHP1)
  462.       ENDIF
  463. *
  464. * Ménage
  465. *
  466.  100  CONTINUE
  467.       SEGSUP IPMAHY
  468.       END
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  

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