Télécharger la plaquette Cast3M


fr en

Télécharger tadve1.eso

Retour à la liste

Numérotation des lignes : 

tadve1

  1. C TADVE1    SOURCE    PV090527  26/04/30    21:16:33     12529          
  2.  
  3. ************************************************************************
  4. *
  5. *                             T A D V E 1
  6. *                             -----------
  7. *
  8. * FONCTION:
  9. * ---------
  10. *     CREATION DE LA MATRICE DE ADVECTION
  11. *     GESTION DES SEGMENTS ET TESTS DE COMPATIBILITE
  12. *
  13. * PARAMETRES:  (E)=ENTREE    (S)=SORTIE
  14. * -----------
  15. *     MMODEL   (E)  POINTEUR SUR LE SEGMENT MMODEL
  16. *     IPCHEL   (E)  POINTEUR SUR LE SEGMENT MCHELM
  17. *     IPRIGI   (S)  POINTEUR SUR LE SEGMENT MRIGID
  18. *
  19. * AUTEUR, DATE DE CREATION:
  20. * -------------------------
  21. *      MARINO ARROYO, 18 MAI 1999
  22. *
  23. * LANGAGE:
  24. * --------
  25. *     ESOPE + FORTRAN77
  26. *
  27. ************************************************************************
  28.  
  29.       SUBROUTINE TADVE1 (MMODEL,IPCHEL,IPRIGI,ISYMM)
  30.  
  31.       IMPLICIT INTEGER(I-N)
  32.       IMPLICIT REAL*8(A-H,O-Z)
  33.  
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. -INC CCHAMP
  37.  
  38. -INC SMCOORD
  39. -INC SMCHAML
  40. -INC SMELEME
  41. -INC SMINTE
  42. -INC SMMODEL
  43.       POINTEUR nomid1.NOMID
  44. -INC SMRIGID
  45.  
  46. -INC TMPTVAL
  47.  
  48.       INTEGER OOOVAL
  49.  
  50.       SEGMENT NOTYPE
  51.         CHARACTER*16 TYPE(NBTYPE)
  52.       ENDSEGMENT
  53.  
  54.       PARAMETER ( NINF=3 )
  55.       INTEGER INFOS(NINF)
  56.       LOGICAL LOG1
  57.  
  58.       CHARACTER*8 CMATE
  59.       CHARACTER*(LCONMO) CONM
  60.       CHARACTER*10 PEAU
  61.       CHARACTER*4 MOTADV
  62.  
  63.       PARAMETER ( NFO1=4,INTTYP=3 )
  64.       CHARACTER*16 MOTFOR,MOTFO1(NFO1)
  65.       DATA MOTFO1 /'THERMIQUE' ,'DIFFUSION','NAVIER_STOKES','MECANIQUE'/
  66.       MACRO,(THERMIQUE,DIFFUSION,NAVIER_STOKES,MECANIQUE)
  67.       DATA MOTADV /'ADVE'/
  68.  
  69.       PARAMETER ( LNUCOQ=5 , LINUM=14 , LINUC=12 )
  70.       INTEGER INUCOQ(LNUCOQ), INUMA(LINUM), INUCO(LINUC)
  71. *
  72. *                TRI3 TRI6 QUA4 QUA8 CUB8 CU20 PRI6 PR15
  73.       DATA INUMA/   4,   6,   8,  10,  14,  15,  16,  17,
  74. *                TET4 TET10 PYR5 PY13 TRI7 QUA9
  75.      &             23,   24,  25,  26, 144, 145/
  76. *                 SEG2 SEG3 TRI3 TRI6 QUA4 QUA8
  77.       DATA INUCO /   2,   3,   4,   6,   8,  10,
  78. *                 RAC2 RAC3 LIA3 LIA6 LIA4 LIA8
  79.      &             12, 13,  18,  19,   20, 21 /
  80. *                  COQ2  COQ3  COQ6 COQ4 COQ8
  81.       DATA INUCOQ / 44 , 27  , 56 , 49  ,41 /
  82.  
  83.       IPRIGI = 0
  84. C---
  85. C Verification du lieu support du MCHAML de caracteristiques
  86. C---
  87.       CALL ecrcha('THERMIQUE')
  88.       CALL ecrcha('FORM')
  89.       CALL ecrobj('MMODEL',MMODEL)
  90.       CALL EXIS
  91.       call lirlog(log1,0,iret1)
  92.  
  93.       if (LOG1) then
  94.       CALL ecrcha('THERMIQUE')
  95.       CALL ecrcha('FORM')
  96.       CALL ecrobj('MMODEL',MMODEL)
  97.       CALL EXTRAI
  98.       CALL lirobj('MMODEL',ipmod1,0,iret1)
  99.        if (iret1.gt.0) then
  100.         CALL QUESUP(ipmod1,IPCHEL,6,0,ISUPCH,IRET)
  101.         IF (ISUPCH.GT.1) RETURN
  102.         goto 7
  103.        endif
  104.       endif
  105.  
  106.       CALL ecrcha('MECANIQUE')
  107.       CALL ecrcha('FORM')
  108.       CALL ecrobj('MMODEL',MMODEL)
  109.       CALL EXIS
  110.       call lirlog(log1,0,iret1)
  111.       if (LOG1) then
  112.       CALL ecrcha('MECANIQUE')
  113.       CALL ecrcha('FORM')
  114.       CALL ecrobj('MMODEL',MMODEL)
  115.       CALL EXTRAI
  116.       CALL lirobj('MMODEL',ipmod2,0,iret2)
  117.        if (iret2.gt.0) then
  118.         CALL QUESUP(ipmod2,IPCHEL,INTTYP,0,ISUP,IRET)
  119.         IF (ISUP.GT.1) RETURN
  120.        endif
  121.       endif
  122.  
  123.  
  124.  
  125.  7    CONTINUE
  126. C---
  127. C Initialisation de la matrice d'ADVECTION (chapeau de l'objet RIGIDITE)
  128. C---
  129.       NRIGEL = 0
  130.       SEGINI,MRIGID
  131.       MTYMAT = 'RIGIDITE'
  132.       ICHOLE = 0
  133.       IMGEO1 = 0
  134.       IMGEO2 = 0
  135.       IFORIG = IFOUR
  136.       ISUPEQ = 0
  137.  
  138. c en cas de besoin
  139.            L1 = 8
  140.            n1 = 1
  141.            segini mmode1
  142.            mchelm = ipchel
  143.            n3 = infche(/2)
  144.            segini mchel1
  145.            mchel1.ifoche = ifoche
  146.            n2 = 1
  147.            segini mcham1
  148.            mchel1.ichaml(1) = mcham1
  149. C---
  150. C BOUCLE SUR LES MODELES ELEMENTAIRES
  151. C---
  152.       NB_OK = 0
  153.       DO 10 III = 1, MMODEL.KMODEL(/1)
  154.         IPINTE = 0
  155.         IPINT1 = 0
  156.         MOMATE = 0
  157.         MOTYPE = 0
  158.  
  159. C- Recuperation du sous-modele et de la zone elementaire associee
  160.         IMODEL = MMODEL.KMODEL(III)
  161.         MOTFOR = IMODEL.FORMOD(1)
  162.         NMAT   = IMODEL.MATMOD(/2)
  163.  
  164. C- Selection uniquement des SOUS-MODELES qui nous interessent
  165.         CALL PLACE(MOTFO1,NFO1,ityp1,MOTFOR)
  166.         IF (ityp1 .EQ. 0) GOTO 10
  167.  
  168.         CASE, ityp1
  169.         WHEN,THERMIQUE,DIFFUSION,MECANIQUE
  170.           CALL PLACE(IMODEL.MATMOD,NMAT,iok3,'ADVECTION       ')
  171.  
  172.         WHENOTHERS
  173.           CALL PLACE(IMODEL.MATMOD,NMAT,iok3,'NLIN            ')
  174.         ENDCASE
  175.         IF (iok3 .EQ. 0) GOTO 10
  176.  
  177.         NB_OK = NB_OK + 1
  178.  
  179. C       Recherche dans MATMOD du mot 'CONSERVATIVE'
  180.         CALL PLACE(IMODEL.MATMOD, NMAT, icons, 'CONSERVATIVE    ')
  181.  
  182. C- Recuperation d'informations sur le maillage elementaire
  183.         IPT1   = IMAMOD
  184.         NBNOE1 = IPT1.NUM(/1)
  185.         NBELE1 = IPT1.NUM(/2)
  186.         IF(NEFMOD.EQ.269 .OR. NEFMOD.EQ.270) THEN
  187.           ITUY = 1
  188.         ELSE
  189.           ITUY = 0
  190.         ENDIF
  191.  
  192. C- Quelques informations et verifications sur le modele elementaire
  193.         CONM  = CONMOD
  194.         CMATE = CMATEE
  195.         MATE  = IMATEE
  196.  
  197. C       Seule l'isotropie est disponible en 2D PLAN et AXISYMETRIQUE
  198.         if(ituy.ne.1 .and.
  199.      &     ityp1.eq.THERMIQUE .or. ityp1.eq.DIFFUSION) then
  200.           IF (MATE.EQ.1) THEN
  201.             IF (IFOMOD.EQ.1) THEN
  202.               WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  203.      &                      'LE CAS FOURIER N''EST PAS PRIS EN COMPTE'
  204.               CALL ERREUR(19)
  205.               GOTO 1999
  206.             ENDIF
  207. C         ELSE
  208. C            WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  209. C     &                   'SEULEMENT LE CAS ISOTROPE EST ENVISAGE'
  210. C            CALL ERREUR(19)
  211. C            GOTO 1999
  212.           ENDIF
  213.         endif
  214. *
  215.         IRET = 1
  216.         CALL IDENT(IPT1,CONM,IPCHEL,0, INFOS,IRET)
  217.         IF (IRET.EQ.0) GOTO 1999
  218. *
  219.         NEF  = NEFMOD
  220.         ICOQ = 0
  221.         CALL PLACE2(INUCOQ,LNUCOQ,ICOQ,NEF)
  222.         IMAS = 0
  223.         CALL PLACE2(INUMA,LINUM,IMAS,NEF)
  224.         IF (IMAS.EQ.0.and . ituy.eq.0) THEN
  225.           WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  226.      &          'SEULS LES ELEMENTS FINIS MASSIFS SONT ENVISAGES'
  227.           CALL ERREUR(19)
  228.           GOTO 1999
  229.         ENDIF
  230.  
  231.  
  232. C- Recuperation des noms des composantes du champ vectoriel (obligatoires)
  233. C  Curieux : On ne tient pas compte de l'AXISYMETRIE et autres cas possibles
  234.         nbrfac = 0
  235.         if(ituy .eq. 0) then
  236. C         Cas des elements MASSIFS
  237.           CASE, ityp1
  238.           WHEN,THERMIQUE
  239. C           Vecteur vitesse + 'RHO' + 'C'
  240.             nbrobl = IDIM + 2
  241.             segini,nomid
  242.             IF     (IDIM.EQ.1) THEN
  243.               lesobl(1) ='VITX'
  244.             ELSEIF (IDIM.EQ.2) THEN
  245.               lesobl(1) ='VITX'
  246.               lesobl(2) ='VITY'
  247.             ELSE
  248.               lesobl(1) ='VITX'
  249.               lesobl(2) ='VITY'
  250.               lesobl(3) ='VITZ'
  251.             ENDIF
  252.             lesobl(IDIM + 1) ='RHO'
  253.             lesobl(IDIM + 2) ='C'
  254.  
  255.           WHEN,DIFFUSION
  256. C           Vecteur vitesse + 'CDIF'
  257.             nbrobl = IDIM + 1
  258.             segini,nomid
  259.             IF     (IDIM.EQ.1) THEN
  260.               lesobl(1) ='VITX'
  261.             ELSEIF (IDIM.EQ.2) THEN
  262.               lesobl(1) ='VITX'
  263.               lesobl(2) ='VITY'
  264.             ELSE
  265.               lesobl(1) ='VITX'
  266.               lesobl(2) ='VITY'
  267.               lesobl(3) ='VITZ'
  268.             ENDIF
  269.             lesobl(IDIM + 1) ='CDIF'
  270.  
  271.           WHEN,NAVIER_STOKES
  272.             nbrobl = 1
  273.             segini,nomid
  274.             lesobl(1) = motadv
  275.  
  276.           WHEN,MECANIQUE
  277. C           Vecteur vitesse + 'RHO'
  278.             nbrobl = IDIM + 1
  279.             segini,nomid
  280.             IF     (IDIM.EQ.1) THEN
  281.               lesobl(1) ='VITX'
  282.             ELSEIF (IDIM.EQ.2) THEN
  283.               lesobl(1) ='VITX'
  284.               lesobl(2) ='VITY'
  285.             ELSE
  286.               lesobl(1) ='VITX'
  287.               lesobl(2) ='VITY'
  288.               lesobl(3) ='VITZ'
  289.             ENDIF
  290.             lesobl(IDIM + 1)='RHO'
  291.           ENDCASE
  292.  
  293.         else
  294. C         Cas des elements TUYAUX : TUY2, TUY3
  295.           CASE, ityp1
  296.           WHEN,THERMIQUE
  297.             nbrobl = 4
  298.             SEGINI,nomid
  299.             lesobl(1)='VITE'
  300.             lesobl(2)='RHO'
  301.             lesobl(3)='C'
  302.             lesobl(4)='SECT'
  303.  
  304.           WHEN,DIFFUSION
  305.             nbrobl = 3
  306.             SEGINI,nomid
  307.             lesobl(1)='VITE'
  308.             lesobl(2)='CDIF'
  309.             lesobl(3)='SECT'
  310.           ENDCASE
  311.        endif
  312.  
  313.        NMATO  = lesobl(/2)
  314.        NMATF  = lesfac(/2)
  315.        NMATT  = NMATO + NMATF
  316.        MOMATE = nomid
  317.  
  318. C      Type des composantes attendues
  319.        nbtype = 1
  320.        SEGINI,notype
  321.        if (ityp1.eq. NAVIER_STOKES) then
  322.          notype.type(1)='POINTEURCHPOINT'
  323.        else
  324.          notype.type(1)='REAL*8'
  325.        endif
  326.        MOTYPE = notype
  327.  
  328.  
  329.        CASE, ityp1
  330.        WHEN,THERMIQUE,DIFFUSION,MECANIQUE
  331. C-       Recuperation d'informations sur l'element fini
  332.          CALL TSHAPE(NEF,'GAUSS',IPINTE)
  333.          IF (IERR.NE.0) GOTO 1999
  334.          MINTE = IPINTE
  335.          SEGACT,MINTE
  336.  
  337. C-       Definition du descripteur IDESCR
  338.          IDESCR = 0
  339.          IF (ICOQ.NE.0 .AND. IF1.NE.0) THEN
  340.            PEAU = '          '
  341.            IF (MATMOD(/1) .NE. 0) PEAU = MATMOD(1)
  342.            CALL TCONV2(ICOQ,PEAU,NBNOE1,IDESCR)
  343.  
  344.          ELSE
  345.            NOMPRI = LNOMID(1)
  346.            NOMDUA = LNOMID(2)
  347.  
  348.            if (ityp1.eq.MECANIQUE) then
  349.              call TARIG1(IMODEL,IDESCR,LRE)
  350.  
  351.            else
  352.              CALL TCOND2(ICOQ,NBNOE1,IDESCR,NOMPRI,NOMDUA)
  353.            endif
  354.  
  355.            descr = IDESCR
  356.            SEGACT,descr
  357.            NLIGRD = LISDUA(/2)
  358.            NLIGRP = LISINC(/2)
  359.            if (ityp1.ne.MECANIQUE) LRE = NLIGRP
  360.          ENDIF
  361.  
  362.  
  363.        WHEN,NAVIER_STOKES
  364.          LRE = 3*NBNOE1
  365.        ENDCASE
  366.  
  367.  
  368. C-     Partionnement si necessaire de la matrice d'ADVECTION
  369. C-     determinant ainsi le nombre d'objets elementaires de MRIGID
  370.        LTRK = oooval(1,4)
  371.        IF (LTRK.EQ.0) LTRK = oooval(1,1)
  372.        LTRK=MAX(LTRK,2**24)
  373.  
  374. *      Ajout a la taille en mots de la matrice des infos du segment
  375.        LSEG   = LRE*LRE*NBELE1 + 16
  376.        NBLPRT = (LSEG-1)/LTRK + 1
  377.        NBLMAX = (NBELE1-1)/NBLPRT + 1
  378.        NBLPRT = (NBELE1-1)/NBLMAX + 1
  379. *      write(ioimp,*) ' tadve1 : nblprt nblmax = ',nblprt,nblmax,nbele1
  380.  
  381. C-     Ajout de la matrice d'ADVECTION a la matrice globale
  382.        NRIGE0 = IRIGEL(/2)
  383.        NRIGEL = NRIGE0 + NBLPRT
  384.        if (ityp1.eq.NAVIER_STOKES) nrigel = nrigel + (idim - 1)*nblprt
  385.        SEGADJ,MRIGID
  386.  
  387.        descr  = IDESCR
  388.        meleme = IPT1
  389.        nbnn   = NBNOE1
  390.        nbelem = NBELE1
  391.        nbsous = 0
  392.        nbref  = 0
  393.  
  394. C====
  395. C Boucle sur les PARTITIONS elementaires de la matrice
  396. C====
  397.         DO 200 irige = 1, NBLPRT
  398.           IF (NBLPRT.GT.1) THEN
  399. C-- Partitionnement du maillage support de la matrice elementaire
  400.             ielem = (irige-1)*NBLMAX
  401.             nbelem = MIN(NBLMAX,NBELE1-ielem)
  402. *            write(ioimp,*) 'tadve1 : creation segment ',nbnn,nbelem
  403.             SEGINI,meleme
  404.             itypel = IPT1.itypel
  405.             DO ielt = 1, nbelem
  406.               jelt = ielt + ielem
  407.               DO inoe = 1, nbnn
  408.                 num(inoe,ielt) = IPT1.NUM(inoe,jelt)
  409.               ENDDO
  410.               icolor(ielt) = IPT1.ICOLOR(jelt)
  411.             ENDDO
  412. C-- Recopie du descripteur
  413.             des1 = IDESCR
  414.             SEGINI,descr=des1
  415.             SEGDES,descr
  416.           ENDIF
  417.           ipmail = meleme
  418.           ipdesc = descr
  419.  
  420. C-- Recuperation des valeurs des caracteristiques necessaires
  421.           IVAMAT = 0
  422.           CALL KOMCHA(IPCHEL,ipmail,CONM,MOMATE,MOTYPE,1,INFOS,3,IVAMAT)
  423.           IF (IERR.NE.0) GOTO 2999
  424.           IF (ISUPCH.EQ.1) THEN
  425.             CALL VALCHE(IVAMAT,NMATT,IPINTE,0,MOMATE,NEF)
  426.             IF (IERR.NE.0) THEN
  427.               ISUPCH = 0
  428.               GOTO 2999
  429.             ENDIF
  430.           ENDIF
  431.  
  432.           if (ityp1.eq.NAVIER_STOKES) then
  433.            segact mmode1*mod
  434.            mmode1.kmodel(1) = imodel
  435.            mchel1.conche(1) = conm
  436.            mchel1.imache(1) = ipmail
  437.            mptval = ivamat
  438.            do jj = 1,n2
  439.              mcham1.nomche(1) = motadv
  440.              mcham1.typche(1) = tyval(1)
  441.              mcham1.ielval(1) = ival(1)
  442.            enddo
  443.  
  444.            ipmons = mmode1
  445.            ipchns = mchel1
  446.            call go2nli(ipmons,ipchns,iprins,4)
  447.            if (ierr.ne.0) return
  448.  
  449.            goto 2999
  450.           endif
  451.  
  452.  
  453. C-- Initialisation de la matrice de rigidite elementaire (xmatri)
  454.           NELRIG = nbelem
  455.           rigrel=0
  456.           SEGINI,xmatri
  457.           ipmatr = xmatri
  458.  
  459. C-- Calcul de la matrice elementaire pour la zone irige (ipmail) et
  460. C-- Remplissage de la matrice globale (ipmatr)
  461.           if(imas.ne.0) then
  462.              CALL TADVE8(NEF,ipmail,IPINTE,MATE,IVAMAT,NMATT,ISYMM,
  463.      &                     ipmatr,LRE,ityp1, icons)
  464.  
  465.           elseif(ituy.ne.0) then
  466.             call adtuy(nef,ipmail,ipinte,mate,ivamat,nmatt, ipmatr,
  467.      &           lre, icons)
  468.  
  469.           else
  470.              call erreur(19)
  471.              return
  472.           endif
  473.  
  474. C--       Un peu de menage dans les segments
  475.  2999     CONTINUE
  476.           IF (ISUPCH.EQ.1 .OR. NBLPRT.NE.1) THEN
  477.             CALL DTMVAL(IVAMAT,3)
  478.           ELSE
  479.             CALL DTMVAL(IVAMAT,1)
  480.           ENDIF
  481.  
  482. C--       Sortie prematuree en cas d'erreur
  483.           IF (IERR.NE.0) GOTO 1999
  484.  
  485.           xmatri = ipmatr
  486.           IF (NBLPRT.GT.1) THEN
  487.             meleme = ipmail
  488.           ENDIF
  489.  
  490.           if (ityp1.eq.NAVIER_STOKES) then
  491.             RI3 = iprins
  492.             segact ri3
  493.             if ( ri3.coerig(/1).ne.idim ) then
  494.               call erreur(5)
  495.               return
  496.             endif
  497.  
  498.             do kige = 1,IDIM
  499.               ipdesc = ri3.IRIGEL(3,kige)
  500.               ipmatr = ri3.IRIGEL(4,kige)
  501.               isymm  = ri3.irigel(7,kige)
  502.  
  503.               jrige = NRIGE0 + kige
  504.               COERIG(jrige)   = ri3.coerig(kige)
  505.               IRIGEL(1,jrige) = ipmail
  506.               IRIGEL(2,jrige) = 0
  507.               IRIGEL(3,jrige) = ipdesc
  508.               IRIGEL(4,jrige) = ipmatr
  509.               IRIGEL(5,jrige) = NIFOUR
  510.               IRIGEL(6,jrige) = 0
  511.               IRIGEL(7,jrige) = 0
  512.               IRIGEL(7,jrige) = ri3.irigel(7,kige)
  513.               IRIGEL(8,jrige) = 0
  514.             enddo
  515.  
  516.           else
  517. C--         Stockage de la matrice
  518.             jrige = NRIGE0 + irige
  519.             COERIG(jrige)   = 1.
  520.             IRIGEL(1,jrige) = ipmail
  521.             IRIGEL(2,jrige) = 0
  522.             IRIGEL(3,jrige) = ipdesc
  523.             IRIGEL(4,jrige) = ipmatr
  524.             IRIGEL(5,jrige) = NIFOUR
  525.             IRIGEL(6,jrige) = 0
  526.             IRIGEL(7,jrige) = 0
  527.             IF (ISYMM.NE.1) IRIGEL(7,jrige) = 2
  528.             xmatri.symre=irigel(7,jrige)
  529.             SEGDES,xmatri
  530.             IRIGEL(8,jrige) = 0
  531.           endif
  532.  200    CONTINUE
  533. C====
  534. C FIN DE LA BOUCLE SUR LES PARTITIONS
  535. C====
  536.  
  537. C- Un peu de menage dans les segments
  538.  1999   CONTINUE
  539.         IF (MOMATE.NE.0) THEN
  540.           nomid = MOMATE
  541.           SEGSUP,nomid
  542.         ENDIF
  543.         IF (MOTYPE.NE.0) THEN
  544.           notype = MOTYPE
  545.           SEGSUP,notype
  546.         ENDIF
  547.         IF (IERR.NE.0) GOTO 999
  548.  10   CONTINUE
  549. C---
  550. C FIN DE LA BOUCLE SUR LES MODELES ELEMENTAIRES
  551. C---
  552.       IF(NB_OK .EQ. 0)THEN
  553.         MOTERR='ADVECTION'
  554.         CALL ERREUR(719)
  555.         RETURN
  556.       ENDIF
  557.  
  558.       IPRIGI = MRIGID
  559.       SEGDES,MRIGID
  560.  
  561.       segsup mmode1,mchel1,mcham1
  562.  
  563.  999  CONTINUE
  564.       RETURN
  565.       END
  566.  
  567.  
  568.  
  569.  
  570.  

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