Télécharger la plaquette Cast3M


fr en

Télécharger matp.eso

Retour à la liste

Numérotation des lignes : 

matp

  1. C MATP      SOURCE    CB215821  24/04/12    21:16:41     11897          
  2.       SUBROUTINE MATP
  3. C-----------------------------------------------------------------------
  4. C Cette subroutine permet de creer les matrices elementaires du système
  5. C matriciel en trace de charge dans le cadre de la résolution des
  6. C équations de DARCY par une méthode d'éléments finis mixtes hybride.
  7. C On obtient ce systeme en exprimant les inconnues de charge et de
  8. C flux en fonction des traces de charge.
  9. C        -1          t              -1
  10. C M1 = M2  * ( COEF * DIV * DIV * M2  - Id )
  11. C
  12. C avec M2    : Matrice massse hybride
  13. C      DIV   : Matrice uniligne representant la divergence
  14. C      COEF  : Scalaire valant 1/d dans le cas permanant, THETA*BETA/d
  15. C              dans le cas transitoire.
  16. C                                        -1  t
  17. C      d     : Scalaire valant ( DIV * M2  *  DIV )
  18. C
  19. C      THETA : Scalaire de la THETA-méthode            |    Données
  20. C      BETA  : Scalaire d*dt/(Ck|K| + THETA*d*dt)      |     utiles
  21. C      dt    : Scalaire pas de temps                   |   uniquement
  22. C      Ck    : Scalaire coefficient d'emmagasinement   |      en
  23. C      |K|   : Scalaire aire de l'élément K            |  transitoire
  24. C
  25. C---------------------------
  26. C Phrase d'appel (GIBIANE) :
  27. C---------------------------
  28. C
  29. C  MRIGI1 = MATP MMODEL  MRIGI2 (TABLE2) ;
  30. C
  31. C------------------------
  32. C Operandes et resultat :
  33. C------------------------
  34. C
  35. C  MRIGI1  : Matrices elementaires en trace de charge.
  36. C  MMODEL  : Objet modele specifiant la formulation.
  37. C  MRIGI2  : Matrices masses hybrides elementaires creees par MHYB.
  38. C  TABLE2  : TABLE DARCY_TRANSITOIRE contenant les infos transitoires.
  39. C
  40. C-----------------------------------------------------------------------
  41. C
  42. C Langage : ESOPE + FORTRAN77
  43. C
  44. C Auteurs : 09/93 F.DABBENE - Cas permanent
  45. C           09/94 X.NOUVELLON - Extension au cas transitoire
  46. C
  47. C-----------------------------------------------------------------------
  48.       IMPLICIT INTEGER(I-N)
  49.       IMPLICIT REAL*8 (A-H,O-Z)
  50. *
  51.  
  52. -INC PPARAM
  53. -INC CCOPTIO
  54. -INC SMELEME
  55. -INC SMRIGID
  56. -INC SMMODEL
  57. -INC SMTABLE
  58. -INC SMCHAML
  59. *
  60.       SEGMENT IPMAHY
  61.          INTEGER MAHYBR(NSOUS)
  62.       ENDSEGMENT
  63. *
  64.       LOGICAL LOGRE,LOGIN
  65.       CHARACTER*6 CHAR6
  66.       CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE
  67. *
  68. * Initialisations
  69. *
  70.       IVALIN = 0
  71.       XVALIN = 0.D0
  72.       LOGIN  = .TRUE.
  73.       IOBIN  = 0
  74.       TAPIND = 'MOT     '
  75. *
  76. * Lecture du MMODEL
  77. *
  78.       CALL LIROBJ('MMODEL',IPMODE,1,IRET1)
  79.       IF (IERR.NE.0) RETURN
  80.       MMODEL = IPMODE
  81. *
  82. * Lecture de la TABLE domaine et récupération des infos
  83. *
  84.       IPTABL = 0
  85.       CALL LEKMOD(IPMODE,IPTABL,IRET)
  86.       IF (IERR.NE.0) RETURN
  87.       CHARIN = 'MAILLAGE'
  88.       TYPOBJ = 'MAILLAGE'
  89.       CALL LEKTAB(IPTABL,CHARIN,IOBRE)
  90.       IF (IERR.NE.0) RETURN
  91.       IPGEOM = IOBRE
  92.       CALL LEKTAB(IPTABL,'ELTFA',IOBRE)
  93.       IF (IERR.NE.0) RETURN
  94.       IELTFA = IOBRE
  95. *
  96. * Lecture de RIGIDITE
  97. *
  98.       CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET1)
  99.       IF (IERR.NE.0) RETURN
  100.       MRIGID = IPRIGI
  101. *
  102. * Lecture éventuelle de la TABLE DARCY_TRANSITOIRE
  103. *
  104.       IPTABL = 0
  105.       CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET)
  106.       IF (IERR.NE.0) RETURN
  107.       IF (IRET.EQ.1) THEN
  108.          TYPOBJ = 'FLOTTANT'
  109.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN,
  110.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  111.          IF (IERR.NE.0) RETURN
  112.          THETA  = XVALRE
  113.          THEMIN = -1.D-12
  114.          THEMAX = 1.D0 - THEMIN
  115.          IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN
  116.             REAERR(1)   = REAL(THETA)
  117.             REAERR(2)   = REAL(0.D0)
  118.             REAERR(3)   = REAL(1.D0)
  119.             MOTERR(1:8) = '  THETA '
  120.             CALL ERREUR(42)
  121.             RETURN
  122.          ENDIF
  123.          IF (THETA.LT.0.D0) THETA=0.D0
  124.          IF (THETA.GT.1.D0) THETA=1.D0
  125.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN,
  126.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  127.          IF (IERR.NE.0) RETURN
  128.          DELTAT = XVALRE
  129.          IF (DELTAT.LT.0.D0) THEN
  130.             REAERR(1)   = REAL(DELTAT)
  131.             REAERR(2)   = REAL(0.D0)
  132.             MOTERR(1:8) = ' DELTAT '
  133.             CALL ERREUR(41)
  134.             RETURN
  135.          ENDIF
  136.          TYPOBJ = 'MCHAML  '
  137.          CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN,
  138.      S                      TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
  139.          CALL ACTOBJ(TYPOBJ,IOBRE,1)
  140.          IF (IERR.NE.0) RETURN
  141.          IPCK = IOBRE
  142.       ELSE
  143.          THETA  = 0.D0
  144.          DELTAT = 0.D0
  145.          IPCK   = 0
  146.       ENDIF
  147. *
  148. * Test de l'ojet MODELE
  149. * Récupération des zones élémentaires ou DARCY est défini
  150. *
  151.       SEGACT MMODEL
  152.       NSOUS  = KMODEL(/1)
  153.       SEGINI IPMAHY
  154.       IDARCY = 0
  155.       IPT2=IPGEOM
  156.       IPT1=IPGEOM
  157.       DO 10 ISOUS=1,NSOUS
  158.          IF(NSOUS.GT.1)THEN
  159.               SEGACT IPT2
  160.               IPT1=IPT2.LISOUS(ISOUS)
  161.          ENDIF
  162.          IMODEL = KMODEL(ISOUS)
  163.          SEGACT IMODEL
  164.          LETYPE = FORMOD(1)
  165.          IF (LETYPE.EQ.'DARCY') THEN
  166.             IDARCY = IDARCY + 1
  167.             MAHYBR(ISOUS) = IPT1
  168.          ENDIF
  169.  10   CONTINUE
  170.       IF (IDARCY.EQ.0) THEN
  171.          MOTERR = LETYPE
  172.          CALL ERREUR(193)
  173.          GOTO 100
  174.       ENDIF
  175. *
  176. * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini
  177. *
  178.       MELEME = IELTFA
  179.       SEGACT MELEME
  180.       LZONES = LISOUS(/1)
  181.       IF (LZONES.EQ.0) LZONES = 1
  182.       IPT1   = IPGEOM
  183.       SEGACT IPT1
  184.       DO 30 ISOUS=1,NSOUS
  185.          IMAMEL = MAHYBR(ISOUS)
  186.          IF (IMAMEL.NE.0) THEN
  187.             DO 20 ISZ=1,LZONES
  188.                IF (LZONES.EQ.1) THEN
  189.                   IPT2 = IPT1
  190.                   IPT3 = MELEME
  191.                ELSE
  192.                   IPT2 = IPT1.LISOUS(ISZ)
  193.                   IPT3 = LISOUS(ISZ)
  194.                ENDIF
  195.                IF (IPT2.EQ.IMAMEL) THEN
  196.                   MAHYBR(ISOUS) = IPT3
  197.                   GOTO 30
  198.                ENDIF
  199.  20         CONTINUE
  200.             IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN
  201.                MOTERR(1:8)  = ' MODELE '
  202.                MOTERR(9:16) = ' ELTFA  '
  203.                INTERR(1)    = ISOUS
  204.                CALL ERREUR(698)
  205.                GOTO 100
  206.             ENDIF
  207.          ENDIF
  208.  30   CONTINUE
  209. *
  210. * Test du sous-type de la matrice de rigidité récupérée
  211. * Controle des MELEME de l'objet RIGIDITE % au ELTFA liée au MMODEL
  212. *
  213.       SEGACT MRIGID
  214.       LETYPE = MTYMAT
  215.       IF (LETYPE.NE.'DARCY') THEN
  216.          MOTERR(1:8)  = 'RIGIDITE'
  217.          MOTERR(9:16) = 'DARCY'
  218.          CALL ERREUR(79)
  219.          SEGDES MRIGID
  220.          GOTO 100
  221.       ENDIF
  222.       DO 40 ISOUS=1,NSOUS
  223.          IF (MAHYBR(ISOUS).NE.0) THEN
  224.             IPTEST = IRIGEL(1,ISOUS)
  225.             IF (IPTEST.NE.MAHYBR(ISOUS)) THEN
  226.                MOTERR(1:8)  = 'DARCY'
  227.                MOTERR(9:16) = ' ELTFA  '
  228.                INTERR(1)    = ISOUS
  229.                CALL ERREUR(698)
  230.                SEGDES MRIGID
  231.                GOTO 100
  232.             ENDIF
  233.          ENDIF
  234.  40   CONTINUE
  235.       SEGDES MRIGID
  236. *
  237. * Vérification du support du MCHAML % au MMODEL
  238. * Controle du nombre de composantes reelles
  239. *
  240.       IF (IPCK.NE.0) THEN
  241.          ISUP  = 2
  242.          ICOND = 1
  243.          CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2)
  244.          IF (IRET.GT.1) GOTO 100
  245.          MCHELM = IPCK
  246.          SEGACT MCHELM
  247.          DO 50 ISOUS=1,NSOUS
  248.             IF (MAHYBR(ISOUS).NE.0) THEN
  249.                MCHAML = ICHAML(ISOUS)
  250.                SEGACT MCHAML
  251.                N2     = TYPCHE(/2)
  252.                IF (N2.GT.1) THEN
  253.                   CALL ERREUR(320)
  254.                   GOTO 100
  255.                ENDIF
  256.                CHAR6  = TYPCHE(1)(1:6)
  257.                IF (CHAR6.NE.'REAL*8') THEN
  258.                   MOTERR(1:8) = NOMCHE(1)
  259.                   CALL ERREUR(679)
  260.                   GOTO 100
  261.                ENDIF
  262.             ENDIF
  263.  50      CONTINUE
  264.       ENDIF
  265.       SEGDES IPMAHY
  266. *
  267. * Construction de IPRIG1
  268. *
  269.       CALL MATP1(NSOUS,IPGEOM,IPMAHY,IPRIGI,THETA,DELTAT,IPCK,IPRIG1)
  270.       CALL ECROBJ('RIGIDITE',IPRIG1)
  271. *
  272. * Ménage
  273. *
  274.  100  CONTINUE
  275.       SEGSUP IPMAHY
  276.       END
  277.  
  278.  
  279.  
  280.  
  281.  

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