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

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