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    CB215821  26/03/06    21:15:08     12485          
  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- Recuperation d'informations sur le maillage elementaire
  180.         IPT1   = IMAMOD
  181.         NBNOE1 = IPT1.NUM(/1)
  182.         NBELE1 = IPT1.NUM(/2)
  183.         IF(NEFMOD.EQ.269 .OR. NEFMOD.EQ.270) THEN
  184.           ITUY = 1
  185.         ELSE
  186.           ITUY = 0
  187.         ENDIF
  188.  
  189. C- Quelques informations et verifications sur le modele elementaire
  190.         CONM  = CONMOD
  191.         CMATE = CMATEE
  192.         MATE  = IMATEE
  193.  
  194. C       Seule l'isotropie est disponible en 2D PLAN et AXISYMETRIQUE
  195.         if(ituy.ne.1 .and.
  196.      &     ityp1.eq.THERMIQUE .or. ityp1.eq.DIFFUSION) then
  197.           IF (MATE.EQ.1) THEN
  198.             IF (IFOMOD.EQ.1) THEN
  199.               WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  200.      &                      'LE CAS FOURIER N''EST PAS PRIS EN COMPTE'
  201.               CALL ERREUR(19)
  202.               GOTO 1999
  203.             ENDIF
  204. C         ELSE
  205. C            WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  206. C     &                   'SEULEMENT LE CAS ISOTROPE EST ENVISAGE'
  207. C            CALL ERREUR(19)
  208. C            GOTO 1999
  209.           ENDIF
  210.         endif
  211. *
  212.         IRET = 1
  213.         CALL IDENT(IPT1,CONM,IPCHEL,0, INFOS,IRET)
  214.         IF (IRET.EQ.0) GOTO 1999
  215. *
  216.         NEF  = NEFMOD
  217.         ICOQ = 0
  218.         CALL PLACE2(INUCOQ,LNUCOQ,ICOQ,NEF)
  219.         IMAS = 0
  220.         CALL PLACE2(INUMA,LINUM,IMAS,NEF)
  221.         IF (IMAS.EQ.0.and . ituy.eq.0) THEN
  222.           WRITE(IOIMP,*) '### ERREUR DANS ADVE : ',
  223.      &          'SEULS LES ELEMENTS FINIS MASSIFS SONT ENVISAGES'
  224.           CALL ERREUR(19)
  225.           GOTO 1999
  226.         ENDIF
  227.  
  228.  
  229. C- Recuperation des noms des composantes du champ vectoriel (obligatoires)
  230. C  Curieux : On ne tient pas compte de l'AXISYMETRIE et autres cas possibles
  231.         nbrfac = 0
  232.         if(ituy .eq. 0) then
  233. C         Cas des elements MASSIFS
  234.           CASE, ityp1
  235.           WHEN,THERMIQUE
  236. C           Vecteur vitesse + 'RHO' + 'C'
  237.             nbrobl = IDIM + 2
  238.             segini,nomid
  239.             IF     (IDIM.EQ.1) THEN
  240.               lesobl(1) ='VITX'
  241.             ELSEIF (IDIM.EQ.2) THEN
  242.               lesobl(1) ='VITX'
  243.               lesobl(2) ='VITY'
  244.             ELSE
  245.               lesobl(1) ='VITX'
  246.               lesobl(2) ='VITY'
  247.               lesobl(3) ='VITZ'
  248.             ENDIF
  249.             lesobl(IDIM + 1) ='RHO'
  250.             lesobl(IDIM + 2) ='C'
  251.  
  252.           WHEN,DIFFUSION
  253. C           Vecteur vitesse + 'CDIF'
  254.             nbrobl = IDIM + 1
  255.             segini,nomid
  256.             IF     (IDIM.EQ.1) THEN
  257.               lesobl(1) ='VITX'
  258.             ELSEIF (IDIM.EQ.2) THEN
  259.               lesobl(1) ='VITX'
  260.               lesobl(2) ='VITY'
  261.             ELSE
  262.               lesobl(1) ='VITX'
  263.               lesobl(2) ='VITY'
  264.               lesobl(3) ='VITZ'
  265.             ENDIF
  266.             lesobl(IDIM + 1) ='CDIF'
  267.  
  268.           WHEN,NAVIER_STOKES
  269.             nbrobl = 1
  270.             segini,nomid
  271.             lesobl(1) = motadv
  272.  
  273.           WHEN,MECANIQUE
  274. C           Vecteur vitesse + 'RHO'
  275.             nbrobl = IDIM + 1
  276.             segini,nomid
  277.             IF     (IDIM.EQ.1) THEN
  278.               lesobl(1) ='VITX'
  279.             ELSEIF (IDIM.EQ.2) THEN
  280.               lesobl(1) ='VITX'
  281.               lesobl(2) ='VITY'
  282.             ELSE
  283.               lesobl(1) ='VITX'
  284.               lesobl(2) ='VITY'
  285.               lesobl(3) ='VITZ'
  286.             ENDIF
  287.             lesobl(IDIM + 1)='RHO'
  288.           ENDCASE
  289.  
  290.         else
  291. C         Cas des elements TUYAUX : TUY2, TUY3
  292.           CASE, ityp1
  293.           WHEN,THERMIQUE
  294.             nbrobl = 4
  295.             SEGINI,nomid
  296.             lesobl(1)='VITE'
  297.             lesobl(2)='RHO'
  298.             lesobl(3)='C'
  299.             lesobl(4)='SECT'
  300.  
  301.           WHEN,DIFFUSION
  302.             nbrobl = 3
  303.             SEGINI,nomid
  304.             lesobl(1)='VITE'
  305.             lesobl(2)='CDIF'
  306.             lesobl(3)='SECT'
  307.           ENDCASE
  308.        endif
  309.  
  310.        NMATO  = lesobl(/2)
  311.        NMATF  = lesfac(/2)
  312.        NMATT  = NMATO + NMATF
  313.        MOMATE = nomid
  314.  
  315. C      Type des composantes attendues
  316.        nbtype = 1
  317.        SEGINI,notype
  318.        if (ityp1.eq. NAVIER_STOKES) then
  319.          notype.type(1)='POINTEURCHPOINT'
  320.        else
  321.          notype.type(1)='REAL*8'
  322.        endif
  323.        MOTYPE = notype
  324.  
  325.  
  326.        CASE, ityp1
  327.        WHEN,THERMIQUE,DIFFUSION,MECANIQUE
  328. C-       Recuperation d'informations sur l'element fini
  329.          CALL TSHAPE(NEF,'GAUSS',IPINTE)
  330.          IF (IERR.NE.0) GOTO 1999
  331.          MINTE = IPINTE
  332.          SEGACT,MINTE
  333.  
  334. C-       Definition du descripteur IDESCR
  335.          IDESCR = 0
  336.          IF (ICOQ.NE.0 .AND. IF1.NE.0) THEN
  337.            PEAU = '          '
  338.            IF (MATMOD(/1) .NE. 0) PEAU = MATMOD(1)
  339.            CALL TCONV2(ICOQ,PEAU,NBNOE1,IDESCR)
  340.  
  341.          ELSE
  342.            NOMPRI = LNOMID(1)
  343.            NOMDUA = LNOMID(2)
  344.  
  345.            if (ityp1.eq.MECANIQUE) then
  346.              call TARIG1(IMODEL,IDESCR,LRE)
  347.  
  348.            else
  349.              CALL TCOND2(ICOQ,NBNOE1,IDESCR,NOMPRI,NOMDUA)
  350.            endif
  351.  
  352.            descr = IDESCR
  353.            SEGACT,descr
  354.            NLIGRD = LISDUA(/2)
  355.            NLIGRP = LISINC(/2)
  356.            if (ityp1.ne.MECANIQUE) LRE = NLIGRP
  357.          ENDIF
  358.  
  359.  
  360.        WHEN,NAVIER_STOKES
  361.          LRE = 3*NBNOE1
  362.        ENDCASE
  363.  
  364.  
  365. C-     Partionnement si necessaire de la matrice d'ADVECTION
  366. C-     determinant ainsi le nombre d'objets elementaires de MRIGID
  367.        LTRK = oooval(1,4)
  368.        IF (LTRK.EQ.0) LTRK = oooval(1,1)
  369.        LTRK=MAX(LTRK,2**24)
  370.  
  371. *      Ajout a la taille en mots de la matrice des infos du segment
  372.        LSEG   = LRE*LRE*NBELE1 + 16
  373.        NBLPRT = (LSEG-1)/LTRK + 1
  374.        NBLMAX = (NBELE1-1)/NBLPRT + 1
  375.        NBLPRT = (NBELE1-1)/NBLMAX + 1
  376. *      write(ioimp,*) ' tadve1 : nblprt nblmax = ',nblprt,nblmax,nbele1
  377.  
  378. C-     Ajout de la matrice d'ADVECTION a la matrice globale
  379.        NRIGE0 = IRIGEL(/2)
  380.        NRIGEL = NRIGE0 + NBLPRT
  381.        if (ityp1.eq.NAVIER_STOKES) nrigel = nrigel + (idim - 1)*nblprt
  382.        SEGADJ,MRIGID
  383.  
  384.        descr  = IDESCR
  385.        meleme = IPT1
  386.        nbnn   = NBNOE1
  387.        nbelem = NBELE1
  388.        nbsous = 0
  389.        nbref  = 0
  390.  
  391. C====
  392. C Boucle sur les PARTITIONS elementaires de la matrice
  393. C====
  394.         DO 200 irige = 1, NBLPRT
  395.           IF (NBLPRT.GT.1) THEN
  396. C-- Partitionnement du maillage support de la matrice elementaire
  397.             ielem = (irige-1)*NBLMAX
  398.             nbelem = MIN(NBLMAX,NBELE1-ielem)
  399. *            write(ioimp,*) 'tadve1 : creation segment ',nbnn,nbelem
  400.             SEGINI,meleme
  401.             itypel = IPT1.itypel
  402.             DO ielt = 1, nbelem
  403.               jelt = ielt + ielem
  404.               DO inoe = 1, nbnn
  405.                 num(inoe,ielt) = IPT1.NUM(inoe,jelt)
  406.               ENDDO
  407.               icolor(ielt) = IPT1.ICOLOR(jelt)
  408.             ENDDO
  409. C-- Recopie du descripteur
  410.             des1 = IDESCR
  411.             SEGINI,descr=des1
  412.             SEGDES,descr
  413.           ENDIF
  414.           ipmail = meleme
  415.           ipdesc = descr
  416.  
  417. C-- Recuperation des valeurs des caracteristiques necessaires
  418.           IVAMAT = 0
  419.           CALL KOMCHA(IPCHEL,ipmail,CONM,MOMATE,MOTYPE,1,INFOS,3,IVAMAT)
  420.           IF (IERR.NE.0) GOTO 2999
  421.           IF (ISUPCH.EQ.1) THEN
  422.             CALL VALCHE(IVAMAT,NMATT,IPINTE,0,MOMATE,NEF)
  423.             IF (IERR.NE.0) THEN
  424.               ISUPCH = 0
  425.               GOTO 2999
  426.             ENDIF
  427.           ENDIF
  428.  
  429.           if (ityp1.eq.NAVIER_STOKES) then
  430.            segact mmode1*mod
  431.            mmode1.kmodel(1) = imodel
  432.            mchel1.conche(1) = conm
  433.            mchel1.imache(1) = ipmail
  434.            mptval = ivamat
  435.            do jj = 1,n2
  436.              mcham1.nomche(1) = motadv
  437.              mcham1.typche(1) = tyval(1)
  438.              mcham1.ielval(1) = ival(1)
  439.            enddo
  440.  
  441.            ipmons = mmode1
  442.            ipchns = mchel1
  443.            call go2nli(ipmons,ipchns,iprins,4)
  444.            if (ierr.ne.0) return
  445.  
  446.            goto 2999
  447.           endif
  448.  
  449.  
  450. C-- Initialisation de la matrice de rigidite elementaire (xmatri)
  451.           NELRIG = nbelem
  452.           SEGINI,xmatri
  453.           ipmatr = xmatri
  454.  
  455. C-- Calcul de la matrice elementaire pour la zone irige (ipmail) et
  456. C-- Remplissage de la matrice globale (ipmatr)
  457.           if(imas.ne.0) then
  458.              CALL TADVE8(NEF,ipmail,IPINTE,MATE,IVAMAT,NMATT,ISYMM,
  459.      &                     ipmatr,LRE,ityp1)
  460.  
  461.           elseif(ituy.ne.0) then
  462.             call adtuy(nef,ipmail,ipinte,mate,ivamat,nmatt, ipmatr,
  463.      &           lre)
  464.  
  465.           else
  466.              call erreur(19)
  467.              return
  468.           endif
  469.  
  470. C--       Un peu de menage dans les segments
  471.  2999     CONTINUE
  472.           IF (ISUPCH.EQ.1 .OR. NBLPRT.NE.1) THEN
  473.             CALL DTMVAL(IVAMAT,3)
  474.           ELSE
  475.             CALL DTMVAL(IVAMAT,1)
  476.           ENDIF
  477.  
  478. C--       Sortie prematuree en cas d'erreur
  479.           IF (IERR.NE.0) GOTO 1999
  480.  
  481.           xmatri = ipmatr
  482.           IF (NBLPRT.GT.1) THEN
  483.             meleme = ipmail
  484.           ENDIF
  485.  
  486.           if (ityp1.eq.NAVIER_STOKES) then
  487.             RI3 = iprins
  488.             segact ri3
  489.             if ( ri3.coerig(/1).ne.idim ) then
  490.               call erreur(5)
  491.               return
  492.             endif
  493.  
  494.             do kige = 1,IDIM
  495.               ipdesc = ri3.IRIGEL(3,kige)
  496.               ipmatr = ri3.IRIGEL(4,kige)
  497.               isymm  = ri3.irigel(7,kige)
  498.  
  499.               jrige = NRIGE0 + kige
  500.               COERIG(jrige)   = ri3.coerig(kige)
  501.               IRIGEL(1,jrige) = ipmail
  502.               IRIGEL(2,jrige) = 0
  503.               IRIGEL(3,jrige) = ipdesc
  504.               IRIGEL(4,jrige) = ipmatr
  505.               IRIGEL(5,jrige) = NIFOUR
  506.               IRIGEL(6,jrige) = 0
  507.               IRIGEL(7,jrige) = 0
  508.               IRIGEL(7,jrige) = ri3.irigel(7,kige)
  509.               IRIGEL(8,jrige) = 0
  510.             enddo
  511.  
  512.           else
  513. C--         Stockage de la matrice
  514.             jrige = NRIGE0 + irige
  515.             COERIG(jrige)   = 1.
  516.             IRIGEL(1,jrige) = ipmail
  517.             IRIGEL(2,jrige) = 0
  518.             IRIGEL(3,jrige) = ipdesc
  519.             IRIGEL(4,jrige) = ipmatr
  520.             IRIGEL(5,jrige) = NIFOUR
  521.             IRIGEL(6,jrige) = 0
  522.             IRIGEL(7,jrige) = 0
  523.             IF (ISYMM.NE.1) IRIGEL(7,jrige) = 2
  524.             xmatri.symre=irigel(7,jrige)
  525.             SEGDES,xmatri
  526.             IRIGEL(8,jrige) = 0
  527.           endif
  528.  200    CONTINUE
  529. C====
  530. C FIN DE LA BOUCLE SUR LES PARTITIONS
  531. C====
  532.  
  533. C- Un peu de menage dans les segments
  534.  1999   CONTINUE
  535.         IF (MOMATE.NE.0) THEN
  536.           nomid = MOMATE
  537.           SEGSUP,nomid
  538.         ENDIF
  539.         IF (MOTYPE.NE.0) THEN
  540.           notype = MOTYPE
  541.           SEGSUP,notype
  542.         ENDIF
  543.         IF (IERR.NE.0) GOTO 999
  544.  10   CONTINUE
  545. C---
  546. C FIN DE LA BOUCLE SUR LES MODELES ELEMENTAIRES
  547. C---
  548.       IF(NB_OK .EQ. 0)THEN
  549.         MOTERR='ADVECTION'
  550.         CALL ERREUR(719)
  551.         RETURN
  552.       ENDIF
  553.  
  554.       IPRIGI = MRIGID
  555.       SEGDES,MRIGID
  556.  
  557.       segsup mmode1,mchel1,mcham1
  558.  
  559.  999  CONTINUE
  560.       RETURN
  561.       END
  562.  
  563.  

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