Télécharger la plaquette Cast3M


fr en

Télécharger dmmu.eso

Retour à la liste

Numérotation des lignes : 

dmmu

  1. C DMMU      SOURCE    CB215821  24/04/12    21:15:37     11897          
  2.       SUBROUTINE DMMU
  3. C-----------------------------------------------------------------------
  4. C                        -1   t
  5. C  Calcul du produit D (M  -U) B  CHPface
  6. C                       -1 t
  7. C    ou   du produit D M    B  CHPface
  8. C         le résultat est un CHPcentre
  9. C-----------------------------------------------------------------------
  10. C
  11. C---------------------------
  12. C Phrase d'appel (GIBIANE) :
  13. C---------------------------
  14. C
  15. C  CHP3 = 'DMMU' MMODEL RIG1  (CHP1) CHP2 ;
  16. C
  17. C------------------------
  18. C Operandes et resultat :
  19. C------------------------
  20. C
  21. C  MMODEL  : modèle formulation DARCY.
  22. C  RIG1    : Matrices hybrides elementaires de DARCY crees par MHYB.
  23. C  CHP1    : flux aux faces
  24. C  CHP2    : CHPO face à plusieur composantes.
  25. C  CHP3    : CHPO centre à plusieur composantes dont le support géometrique
  26. C            est le maillage centre du domaine les noms des composantes
  27. C            sont ceux de CHP2
  28. C
  29. C-----------------------------------------------------------------------
  30. C
  31. C Langage : ESOPE + FORTRAN77
  32. C
  33. C
  34. C-----------------------------------------------------------------------
  35.       IMPLICIT INTEGER(I-N)
  36.       IMPLICIT REAL*8 (A-H,O-Z)
  37. *
  38.  
  39. -INC PPARAM
  40. -INC CCOPTIO
  41. -INC SMCOORD
  42. -INC SMCHPOI
  43. -INC SMCHAML
  44. -INC SMELEME
  45. -INC SMRIGID
  46. -INC SMTABLE
  47. -INC SMMODEL
  48. *
  49.       SEGMENT ICCPR
  50.          INTEGER ICPR(NNGOT)
  51.       ENDSEGMENT
  52.       SEGMENT IORGA
  53.          INTEGER  IAA(ITES)
  54.       ENDSEGMENT
  55. *
  56.       LOGICAL LOGRE,LOGIN
  57.       CHARACTER*8 TAPIND,TYPOBJ,CHARRE,LETYPE
  58.       CHARACTER*4 NOMTOT(1)
  59. *
  60. * Initialisations
  61. *
  62.       IVALIN = 0
  63.       XVALIN = 0.D0
  64.       LOGIN  = .TRUE.
  65.       IOBIN  = 0
  66.       TAPIND = 'MOT     '
  67. *
  68. *
  69. * Lecture du MMODEL
  70. *
  71.       CALL LIROBJ('MMODEL',IPMODE,1,IRET)
  72.       IF(IERR.NE.0)RETURN
  73.       MMODEL=IPMODE
  74.       SEGACT MMODEL
  75.       N1=KMODEL(/1)
  76.       DO 7 I=1,N1
  77.          IMODEL=KMODEL(I)
  78.          SEGACT IMODEL
  79.          IF(FORMOD(1).NE.'DARCY')THEN
  80.              MOTERR(1:16)  = 'DARCY           '
  81.              CALL ERREUR(719)
  82.              RETURN
  83.          ENDIF
  84.     7 CONTINUE
  85. C  on récupère la table DOMAINE
  86.       IPTABL = 0
  87.       CALL LEKMOD(MMODEL,IPTABL,IRET)
  88.       IF (IERR.NE.0) RETURN
  89.       TYPOBJ = 'MAILLAGE'
  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.       IFACE = 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. *
  107. *
  108. * Test du sous-type de la matrice de rigiditée récupérée
  109. *
  110.       SEGACT MRIGID
  111.       LETYPE = MTYMAT
  112.       IF (LETYPE.NE.'DARCY') THEN
  113.          MOTERR(1:8)  = 'RIGIDITE'
  114.          MOTERR(9:16) = 'DARCY   '
  115.          CALL ERREUR(79)
  116.          SEGDES MRIGID
  117.          GOTO 100
  118.       ENDIF
  119. *
  120. * Controle des pointeurs de MELEME de la rigidité
  121. *
  122.       NRIGEL=IRIGEL(/2)
  123. C     write(6,*)' NRIGEL----- ' ,nrigel
  124.       MELEME=IELTFA
  125.       SEGACT MELEME
  126.       NBSOUS=LISOUS(/1)
  127.       IF(NBSOUS.EQ.0)THEN
  128.           IF((NRIGEL.NE.1).OR.(IRIGEL(1,1).NE.MELEME))THEN
  129.                MOTERR(1:8)  = 'DARCY   '
  130.                MOTERR(9:16) = 'ELTFA   '
  131.                INTERR(1)    = 1
  132.                CALL ERREUR(698)
  133.                SEGDES MRIGID
  134.                GOTO 100
  135.           ENDIF
  136.       ELSE
  137.           IF(NRIGEL.NE.NBSOUS)THEN
  138.                MOTERR(1:8)  = 'DARCY   '
  139.                MOTERR(9:16) = 'ELTFA   '
  140.                INTERR(1)    = 1
  141.                CALL ERREUR(698)
  142.                SEGDES MRIGID
  143.                GOTO 100
  144.           ENDIF
  145.           DO 10 ISOUS=1,NBSOUS
  146.             IF (LISOUS(ISOUS).NE.IRIGEL(1,ISOUS)) THEN
  147.                MOTERR(1:8)  = 'DARCY   '
  148.                MOTERR(9:16) = 'ELTFA   '
  149.                INTERR(1)    = ISOUS
  150.                CALL ERREUR(698)
  151.                SEGDES MRIGID
  152.                GOTO 100
  153.             ENDIF
  154.  10       CONTINUE
  155.       ENDIF
  156. *
  157. *              LECTURE DES FLUX  OU
  158. *               DU CHPOIN FACE A MULTIPLIER
  159. *
  160.       CALL LIROBJ('CHPOINT',IPCHF,1,IRET)
  161.       IF(IRET.NE.1) RETURN
  162.       MCHPO1=IPCHF
  163.       NOMTOT(1)='    '
  165.       NBCOMP=-1
  166.       CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT)
  167.       IF(IERR.NE.0)RETURN
  168.       SEGACT MCHPO1
  169.       MSOUP1=MCHPO1.IPCHP(1)
  170.       SEGACT MSOUP1
  171.       IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN
  172.             MPOVA1=MSOUP1.IPOVAL
  173.             SEGACT MPOVA1
  174. *
  175. *
  176.             CALL LIROBJ('CHPOINT',IPCHF,1,IRET)
  177.             IF(IRET.NE.1) RETURN
  178.             MCHPO2=IPCHF
  179.             NOMTOT(1)='    '
  181.             NBCOMP=-1
  182.             CALL QUEPOI(MCHPO2,IFACE,INDIC,NBCOMP,NOMTOT)
  183.             IF(IERR.NE.0)RETURN
  184.             SEGACT MCHPO2
  185.             MSOUP2=MCHPO2.IPCHP(1)
  186.             SEGACT MSOUP2
  187.             MPOVA2=MSOUP2.IPOVAL
  188.             SEGACT MPOVA2
  189.       ELSE
  190.             MCHPO2=MCHPO1
  191.             MCHPO1=0
  192.             MSOUP2=MSOUP1
  193.             MPOVA2=MSOUP2.IPOVAL
  194.             SEGACT MPOVA2
  195.             CALL LIROBJ('CHPOINT',IPCHF,0,IRET)
  196.             IF(IRET.EQ.1) THEN
  197.                   MCHPO1=IPCHF
  198.                   NOMTOT(1)='    '
  200.                   NBCOMP=-1
  201.                   CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT)
  202.                   IF(IERR.NE.0)RETURN
  203.                   SEGACT MCHPO1
  204.                   MSOUP1=MCHPO1.IPCHP(1)
  205.                   SEGACT MSOUP1
  206.                   IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN
  207.                         MPOVA1=MSOUP1.IPOVAL
  208.                         SEGACT MPOVA1
  209.                   ELSE
  210.                         CALL ERREUR(21)
  211.                         RETURN
  212.                   ENDIF
  213.             ENDIF
  214.       ENDIF
  215. *
  216. * Construction de MCHPOI
  217. *
  218. *
  219.       IPT1=ICENTR
  220.       SEGACT IPT1
  221.       NPN=IPT1.NUM(/2)
  222.       NSOUPO=1
  223.       NAT=1
  224.       SEGINI MCHPOI
  225.       MTYPOI='        '
  226.       MOCHDE='      CHPOIN CREE PAR DMMU                      '
  227.       IFOPOI=IFOUR
  228.       JATTRI(1)=2
  229.       NC=NBCOMP
  230.       SEGINI MSOUPO
  231.       IPCHP(1)=MSOUPO
  232.       DO 5 L=1,NBCOMP
  233.             NOCOMP(L)=MSOUP2.NOCOMP(L)
  234.             NOHARM(L)=MSOUP2.NOHARM(L)
  235.     5 CONTINUE
  236.       IGEOC=ICENTR
  237.       N=NPN
  238.       SEGINI MPOVAL
  239.       IPOVAL=MPOVAL
  240.       NB=N*NC
  241.       CALL INITD(VPOCHA,NB,0.D0)
  242. *
  243. *       Creation du tableau ICPR
  244. *
  245.       IK     = 0
  246.       NNGOT  = nbpts
  247.       SEGINI ICCPR
  248.       MELEME = IFACE
  249.       SEGACT MELEME
  250.       N2 = NUM(/2)
  251.       DO 15 I2=1,N2
  252.          K = NUM(1,I2)
  253.          IF (ICPR(K).EQ.0) THEN
  254.             IK      = IK + 1
  255.             ICPR(K) = IK
  256.          ENDIF
  257.   15  CONTINUE
  258.       ITES=IK
  259.       SEGDES MELEME
  260.       SEGINI IORGA
  261.             MELEME = MSOUP2.IGEOC
  262.             SEGACT MELEME
  263.             N2 = NUM(/2)
  264.             DO 25 I2=1,N2
  265.                K = NUM(1,I2)
  266.                IF (ICPR(K).EQ.0) THEN
  267.                   INTERR(1)   = K
  268.                   MOTERR(1:8) = 'FACE    '
  269.                   CALL ERREUR(64)
  270.                   SEGSUP ICCPR, IORGA
  271.                   RETURN
  272.                ELSE
  273.                   IAA(ICPR(K)) = I2
  274.                ENDIF
  275.  25         CONTINUE
  276.       SEGDES MELEME
  277.       CALL LEKTAB(IPTABL,'XXNORMAE',IORIE)
  278.       MCHELM=IORIE
  279.       SEGACT  MCHELM
  280. C
  281. C     Calcul du produit
  282. C
  283.       MELEME=IELTFA
  284.       SEGACT MELEME
  285.       IPT1=MELEME
  286.       ITELEM=0
  287.       IF(MCHPO1.NE.0)THEN
  288. C
  289. C                     Cas avec vitesses
  290. C
  291.       DO 50 ISOUS=1,NRIGEL
  292.       MCHAML= ICHAML(ISOUS)
  293.       SEGACT MCHAML
  294.       MELVAL=IELVAL(1)
  295.       SEGACT MELVAL
  296.          IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS)
  297.          SEGACT IPT1
  298.          xMATRI=IRIGEL(4,ISOUS)
  299.          SEGACT xMATRI
  300.          NELRIG=re(/3)
  301. C        write(6,*)' ISOUS NELRIG ',ISOUS,NELRIG
  302.          DO 60 IEL=1,NELRIG
  303.          ITELEM=ITELEM+1
  304. *         XMATRI=IMATTT(IEL)
  305. *         SEGACT XMATRI
  306.          NLIGRD=RE(/1)
  307.          NLIGRP=RE(/2)
  308. C        write(6,*) ' NLIGRD NLIGRP ',NLIGRD,NLIGRP
  309.          DO 40 I=1,NLIGRD
  310.          RLIGN=0.D0
  311.          DO 30 J=1,NLIGRP
  312.          RLIGN=RLIGN+RE(I,J,iel)
  313.    30    CONTINUE
  314.          IPOPTS=ICPR(IPT1.NUM(I,IEL))
  315.          III=IAA(IPOPTS)
  316. C        write(6,*)' dmmu ',itelem,i,rlign,MPOVA1.VPOCHA(IPOPTS,1)
  317. C    *  , VELCHE(I,IEL)
  318.          RLIGN=RLIGN-MPOVA1.VPOCHA(IPOPTS,1)*VELCHE(I,IEL)
  319.          DO 20 K=1,NBCOMP
  320.               VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN*
  321.      *                   MPOVA2.VPOCHA(III,K)
  322.    20    CONTINUE
  323.  
  324.    40    CONTINUE
  325. *         SEGDES XMATRI
  326.    60    CONTINUE
  327.          SEGDES xMATRI,MCHAML,MELVAL
  328.    50 CONTINUE
  329.       ELSE
  330. C
  331. C                      Cas ou il n'y a pas de vitesse
  332. C         On a dupliqué toute la boucle 50 pour ne pas faire
  333. C         NRIGEL*NELRIG*NLIGRD fois le meme test
  334. C
  335.       DO 51 ISOUS=1,NRIGEL
  336.       MCHAML= ICHAML(ISOUS)
  337.       SEGACT MCHAML
  338.       MELVAL=IELVAL(1)
  339.       SEGACT MELVAL
  340.          IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS)
  341.          SEGACT IPT1
  342.          xMATRI=IRIGEL(4,ISOUS)
  343.          SEGACT xMATRI
  344.          NELRIG=re(/3)
  345.          DO 61 IEL=1,NELRIG
  346.          ITELEM=ITELEM+1
  347. *         XMATRI=IMATTT(IEL)
  348. *         SEGACT XMATRI
  349.          NLIGRD=RE(/1)
  350.          NLIGRP=RE(/2)
  351.          DO 41 I=1,NLIGRD
  352.          RLIGN=0.D0
  353.          DO 31 J=1,NLIGRP
  354.          RLIGN=RLIGN+RE(I,J,iel)
  355.    31    CONTINUE
  356.          IPOPTS=ICPR(IPT1.NUM(I,IEL))
  357.          III=IAA(IPOPTS)
  358.          DO 21 K=1,NBCOMP
  359.               VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN*
  360.      *                   MPOVA2.VPOCHA(III,K)
  361.    21    CONTINUE
  362.  
  363.    41    CONTINUE
  364. *         SEGDES XMATRI
  365.    61    CONTINUE
  366.          SEGDES xMATRI,MCHAML,MELVAL
  367.    51 CONTINUE
  368.       ENDIF
  369.       SEGDES MRIGID
  370.        CALL ECROBJ('CHPOINT',MCHPOI)
  371. *
  372. * Ménage
  373. *
  374.       SEGDES MELEME
  375.       SEGSUP ICCPR ,IORGA
  376.       SEGDES MPOVAL,MSOUPO,MCHPOI
  377.       SEGDES MPOVA2,MSOUP2,MCHPO2
  378.       IF(MCHPO1.NE.0) THEN
  379.             SEGDES MPOVA1,MSOUP1,MCHPO1
  380.       ENDIF
  381.  100  CONTINUE
  382.       RETURN
  383.       END
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  

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