Télécharger la plaquette Cast3M


fr en

Télécharger hybp1.eso

Retour à la liste

Numérotation des lignes : 

hybp1

  1. C HYBP1     SOURCE    FANDEUR   22/01/03    21:15:20     11136          
  2.       SUBROUTINE HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE,
  3.      S     ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,IRET)
  4. C-----------------------------------------------------------------------
  5. C Calcul de la charge au centre de chaque maille
  6. C-----------------------------------------------------------------------
  7. C
  8. C---------------------------
  9. C Parametres Entree/Sortie :
  10. C---------------------------
  11. C
  12. C E/  IPMAHY  : Segment contenant le pointeur vers le meleme des
  13. C               connectivites elements/faces pour les zones du MMODEL
  14. C               ou on a defini DARCY.
  15. C E/  ICENTR  : Pointeur vers l'objet maillage CENTRE
  16. C E/  IPFACE  : Pointeur vers l'objet maillage FACE       -1
  17. C E/  IPDARC  : Matrice de sous type DARCY (contient RE  ).
  18. C E/  ITPN    : CHPO face des Traces de charge TH au temps n
  19. C E/  THETA   : Parametre de discretisation temporelle (Theta-methode)
  20. C E/  DELTAT  : Pas de temps
  21. C E/  IPCK    : MCHAML donnant pour chaque element Ck|K|
  22. C E/  ITP     : CHPO face des traces de charge TH au temps n-1
  23. C E/  IH      : CHPO centre des charges H au temps n-1
  24. C E/  ISOUR   : CHPO centre SOUR contenant un eventuel terme source
  25. C E/  IPFORC  : Matrice de sous type MASSE (contient RE  ).
  26. C E/  IFORC   : CHPO centre force contenant la force volumique
  27. C  /S IRET    : CHPO centre H resultat contenant la charge au temps n
  28. C
  29. C----------------------
  30. C Variables en COMMON :
  31. C----------------------
  32. C
  33. C E/  IDIM    : cf CCOPTIO
  34. C E/  IFOUR   : cf CCOPTIO
  35. C E/  NIFOUR  : cf CCOPTIO
  36. C
  37. C----------------------
  38. C Tableaux de travail :
  39. C----------------------
  40. C
  41. C  ICPR(I)=J   : Le noeud I a le numero local J
  42. C                Correspondance numerotation globale/locale
  43. C  NNGOT       : Nombre de noeuds total du domaine
  44. C
  45. C-----------------------------------------------------------------------
  46. C
  47. C Langage : ESOPE + FORTRAN77
  48. C
  49. C Auteurs : 09/93 F.DABBENE - Cas permanent
  50. C           09/94 X.NOUVELLON - Extension au cas transitoire
  51. C           02/96 L.V.BENET - Prise en compte de forces de volume
  52. C
  53. C-----------------------------------------------------------------------
  54.       IMPLICIT INTEGER(I-N)
  55.       IMPLICIT REAL*8 (A-H,O-Z)
  56. *
  57.  
  58. -INC PPARAM
  59. -INC CCOPTIO
  60. -INC SMCHAML
  61. -INC SMCHPOI
  62.       POINTEUR MCHPO5.MCHPOI,MCHPO6.MCHPOI,MSOUP6.MSOUPO
  63. -INC SMCOORD
  64. -INC SMELEME
  65. -INC SMRIGID
  66. *
  67.       CHARACTER*(LOCOMP) MOREFD, MOREFP
  68.       CHARACTER*1 BLAN1
  69.       SEGMENT ICCPR
  70.       INTEGER ICPR(NNGOT)
  71.       ENDSEGMENT
  72.       SEGMENT IORGA
  73.       REAL*8  VAA(ITES)
  74.       ENDSEGMENT
  75.       SEGMENT IPMAHY
  76.       INTEGER MAHYBR(NSOUS)
  77.       ENDSEGMENT
  78.       SEGMENT MTRAV
  79.       REAL*8  VECT(NBDDL),FLFOR(NBDDL),RFOR(NBDDL)
  80.       ENDSEGMENT
  81. *
  82.       BLAN1=' '
  83. *
  84. * Initialisations
  85. *
  86.       IRET   = 0
  87.       IPT1   = IPFACE
  88.       MRIGID  = IPDARC
  89.       RI1    = IPFORC
  90.       MCHPOI = ITPN
  91.       MCHELM = IPCK
  92.       MCHEL1 = IORIE
  93.       MCHPO1 = ISOUR
  94.       MCHPO2 = ITP
  95.       MCHPO3 = IH
  96.       MCHPO4 = IFORC
  97.       MCHPO5 = INORM
  98.       MCHPO6 = ISURF
  99.       VSOUR  = 0.D0
  100. *
  101. * Initialisations se servant de IPMAHY et MRIGID
  102. *
  103.       SEGACT IPMAHY
  104.       NBMAIL = MAHYBR(/1)
  105.       SEGACT MRIGID
  106.       DO 5 IMAIL=1,NBMAIL
  107.          IF (MAHYBR(IMAIL).NE.0) THEN
  108.             DESCR  = IRIGEL(3,IMAIL)
  109.             SEGACT DESCR
  110.             MOREFP = LISINC(1)
  111.             SEGDES DESCR
  112.             GOTO 10
  113.          ENDIF
  114.  5    CONTINUE
  115.  10   CONTINUE
  116.       NNGOT  = nbpts
  117.       SEGINI ICCPR
  118. *
  119. *= Creation du tableau ICPR pour le maillage IPT1
  120. *
  121.       SEGACT IPT1
  122.       N2 = IPT1.NUM(/2)
  123.       IK = 0
  124.       DO 15 I2=1,N2
  125.          K = IPT1.NUM(1,I2)
  126.          IF (ICPR(K).EQ.0) THEN
  127.             IK      = IK + 1
  128.             ICPR(K) = IK
  129.          ENDIF
  130.  15   CONTINUE
  131.       ITES = IK
  132. *
  133. *= Extraction des TH et réduction sur le support IPT1 du CHPO ITPN
  134. *= Les valeurs du CHPOINT en TH sont dans VAA avec l'ordre IPT1
  135. *
  136.       SEGINI IORGA
  137.       SEGACT MCHPOI
  138.       NSOUPO = IPCHP(/1)
  139.       DO 30 ISOU=1,NSOUPO
  140.          MSOUPO = IPCHP(ISOU)
  141.          SEGACT MSOUPO
  142.          NC     = NOCOMP(/2)
  143.          NUMHAR = NOHARM(1)
  144.          ITTP = 0
  145.          DO 20 I=1,NC
  146.             IF (NOCOMP(I).EQ.MOREFP) THEN
  147.                ITTP = I
  148.             ENDIF
  149.  20      CONTINUE
  150.          IF (ITTP.NE.0) THEN
  151.             MELEME = IGEOC
  152.             SEGACT MELEME
  153.             N2 = NUM(/2)
  154.             MPOVAL = IPOVAL
  155.             SEGACT MPOVAL
  156.             DO 25 I2=1,N2
  157.                K = NUM(1,I2)
  158.                IF (ICPR(K).EQ.0) THEN
  159.                   INTERR(1)   = K
  160.                   MOTERR(1:8) = 'FACE    '
  161.                   CALL ERREUR(64)
  162.                   SEGDES MRIGID, IPMAHY
  163.                   SEGSUP ICCPR, IORGA
  164.                   RETURN
  165.                ELSE
  166.                   VAA(ICPR(K)) = VPOCHA(I2,ITTP)
  167.                ENDIF
  168.  25         CONTINUE
  169.          ENDIF
  170.  30   CONTINUE
  171. *
  172. *= Création des chapeaux du CHPO centre P résultat
  173. *
  174.       NSOUPO = 1
  175.       NAT    = 1
  176.       SEGINI MCHPOI
  177.       IRET      = MCHPOI
  178.       DO 31 ITIT=1,72
  179.          MOCHDE(ITIT:ITIT)=BLAN1
  180.  31   CONTINUE
  181.       MTYPOI    = 'CENTRE  '
  182.       IFOPOI    = IFOUR
  183.       JATTRI(1) = 1
  184.       NC        = 1
  185.       SEGINI MSOUPO
  186.       IPCHP(1)  = MSOUPO
  187.       NOHARM(1) = NUMHAR
  188.       IGEOC     = ICENTR
  189.       NOCOMP(1) = 'H   '
  190.       IPT2 = ICENTR
  191.       SEGACT IPT2
  192.       N    = IPT2.NUM(/2)
  193.       SEGINI MPOVAL
  194.       IPOVAL = MPOVAL
  195. *
  196. * Activation du MPOVAL du CHPO centre SOUR
  197. *
  198.       IF (ISOUR.NE.0) THEN
  199.          SEGACT MCHPO1
  200.          MSOUP1 = MCHPO1.IPCHP(1)
  201.          SEGACT MSOUP1
  202.          MPOVA1 = MSOUP1.IPOVAL
  203.          SEGACT MPOVA1
  204.       ENDIF
  205. *
  206. * Activation du MPOVAL du CHPO face trace de charge au temps (n-1)
  207. *
  208.       IF (ITP.NE.0) THEN
  209.          SEGACT MCHPO2
  210.          MSOUP2 = MCHPO2.IPCHP(1)
  211.          SEGACT MSOUP2
  212.          MPOVA2 = MSOUP2.IPOVAL
  213.          SEGACT MPOVA2
  214.       ENDIF
  215. *
  216. * Activation du MPOVAL du CHPO centre charge au temps (n-1)
  217. *
  218.       IF (IH.NE.0) THEN
  219.          SEGACT MCHPO3
  220.          MSOUP3 = MCHPO3.IPCHP(1)
  221.          SEGACT MSOUP3
  222.          MPOVA3 = MSOUP3.IPOVAL
  223.          SEGACT MPOVA3
  224.       ENDIF
  225. *
  226. * activation des objets liés à la présence d'une force volumique
  227. *
  228.       IF (IFORC.NE.0) THEN
  229. *
  230. * Activation du MPOVAL du CHPO force appuyé au centre des éléments volumiques
  231. *
  232.          SEGACT MCHPO4
  233.          MSOUP4 = MCHPO4.IPCHP(1)
  234.          SEGACT MSOUP4
  235.          MPOVA4 = MSOUP4.IPOVAL
  236.          SEGACT MPOVA4
  237. *
  238. * Activation du MPOVAL du CHPO des vecteurs normales appuyé au centre des faces
  239. *
  240.          SEGACT MCHPO5
  241.          MSOUP5 = MCHPO5.IPCHP(1)
  242.          SEGACT MSOUP5
  243.          MPOVA5 = MSOUP5.IPOVAL
  244.          SEGACT MPOVA5
  245. *
  246. * Activation du MPOVAL du CHPO des surfaces appuyé au centre des faces
  247. *
  248.          SEGACT MCHPO6
  249.          MSOUP6 = MCHPO6.IPCHP(1)
  250.          SEGACT MSOUP6
  251.          MPOVA6 = MSOUP6.IPOVAL
  252.          SEGACT MPOVA6
  253. *
  254. * Activation du MCHEL des orientations des faces
  255. *
  256.          SEGACT MCHEL1
  257. *
  258. * Activation du MRIGI de la matrice masse hybride
  259. *
  260.          SEGACT RI1
  261.       ENDIF
  262. *
  263. * Activation de Ck|K|
  264. *
  265.       IF (IPCK.NE.0) SEGACT MCHELM
  266. *
  267. *---------------------------------------
  268. *= Boucle sur les maillages elementaires
  269. *---------------------------------------
  270. *
  271.       ITELEM = 0
  272.       DO 110 IMAIL=1,NBMAIL
  273.          IF (MAHYBR(IMAIL).EQ.0) GOTO 110
  274. *
  275. * Récupération des infos pour la zone IMAIL ou Darcy est définie
  276. *
  277.          MELEME = IRIGEL(1,IMAIL)
  278.          xMATRI = IRIGEL(4,IMAIL)
  279.          SEGACT MELEME
  280.          NBDDL  = NUM(/1)
  281.          NBELEM = NUM(/2)
  282.          SEGACT xMATRI
  283.          SEGINI MTRAV
  284. *
  285. * Activation des objets necessaires à la prise en compte des forces de volumes
  286. *
  287.          IF (IFORC.NE.0) THEN
  288.             MCHAM1 = MCHEL1.ICHAML(IMAIL)
  289.             SEGACT MCHAM1
  290.             MELVA1 = MCHAM1.IELVAL(1)
  291.             SEGACT MELVA1
  292.             xMATR1 = RI1.IRIGEL(4,IMAIL)
  293.             SEGACT xMATR1
  294.          ELSE
  295.             DO 35 I=1,NBDDL
  296.                RFOR(I)=0.D0
  297.  35         CONTINUE
  298.          ENDIF
  299. *
  300. C
  301. C Activation du MELVAL du MCHAML Ck|K|
  302. C
  303.          IF (IPCK.NE.0) THEN
  304.             MCHAML = ICHAML(IMAIL)
  305.             SEGACT MCHAML
  306.             MELVAL = IELVAL(1)
  307.             SEGACT MELVAL
  308.          ENDIF
  309. *
  310. *-------------------------------------------------------
  311. * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL
  312. *-------------------------------------------------------
  313. *
  314.          DO 100 IEL=1,NBELEM
  315.             ITELEM = ITELEM + 1
  316. *
  317. *- Calcul de VECT = DIV * RE  , CONSD = VECT * DIV
  318. *
  319. *            XMATRI = IMATTT(IEL)
  320. *            SEGACT XMATRI
  321.             CONSD = 0.D0
  322.             DO 45 J=1,NBDDL
  323.                VECT(J) = 0.D0
  324.                DO 40 I=1,NBDDL
  325.                   VECT(J) = VECT(J) + RE(I,J,iel)
  326.  40            CONTINUE
  327.                CONSD = CONSD + VECT(J)
  328.  45         CONTINUE
  329. *            SEGDES XMATRI
  330. *
  331. *- Recuperation de SOUR et de H(n-1)
  332. *
  333.             IF (ISOUR.NE.0) VSOUR=MPOVA1.VPOCHA(ITELEM,1)
  334.             IF (IFORC.NE.0) THEN
  335. *
  336. *- calcul des flux de forces aux faces de l'element
  337. *
  338.                DO 55 IDDL=1,NBDDL
  339.                   FLFOR(IDDL)= 0.D0
  340.                   IPOPTS = ICPR(NUM(IDDL,IEL))
  341.                   DO 50 I=1,IDIM
  342.                      FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I)
  343.      $                    *MELVA1.VELCHE(IDDL,IEL) * MPOVA4
  344.      $                    .VPOCHA(ITELEM,I) *MPOVA6.VPOCHA(IPOPTS,1)
  345.  50               CONTINUE
  346.  55            CONTINUE
  347. *
  348. *- Construction du tableau aux faces M.FORCE
  349. *
  350. *               XMATR1 = IMATR1.IMATTT(IEL)
  351. *               SEGACT XMATR1
  352.                DO 65 I=1,NBDDL
  353.                   RFOR(I)=0.D0
  354.                   DO 60 J=1,NBDDL
  355.                      RFOR(I)   = RFOR(I) + XMATR1.RE(I,J,iel)*FLFOR(J)
  356.  60               CONTINUE
  357.  65            CONTINUE
  358. *               SEGDES XMATR1
  359.             ENDIF
  360. *
  361. *- Calcul de p pour l'element IEL
  362. *
  363.             VALP = 0.D0
  364.             IF (IPCK.EQ.0) THEN
  365.                DO 80 I=1,NBDDL
  366.                   NUMFA = ICPR(NUM(I,IEL))
  367.                   VALP  = VALP + VECT(I) * (VAA(NUMFA) - RFOR(I))
  368.  80            CONTINUE
  369.                VPOCHA(ITELEM,1) = (VALP + VSOUR) / CONSD
  370.             ELSE
  371.                BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT)
  372.                DO 90 I=1,NBDDL
  373.                   NUMFA = ICPR(NUM(I,IEL))
  374.                   VALP  = VALP + VECT(I) * ( THETA*VAA(NUMFA)
  375.      S                 + (1.D0-THETA)*MPOVA2.VPOCHA(NUMFA,1))
  376.  90            CONTINUE
  377.                VALP = (VALP + VSOUR) * BETA / CONSD
  378.                VH   = MPOVA3.VPOCHA(ITELEM,1)
  379.                VPOCHA(ITELEM,1) = VALP + (1.D0-BETA)*VH
  380.             ENDIF
  381.  100     CONTINUE
  382.          IF (IFORC.NE.0)  SEGDES xMATR1
  383.          SEGDES xMATRI
  384.          SEGSUP MTRAV
  385.  110  CONTINUE
  386.       SEGDES MRIGID, IPMAHY
  387.       IF (IFORC.NE.0) SEGDES RI1
  388.       SEGSUP ICCPR , IORGA
  389. *
  390.       END
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  

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