Télécharger la plaquette Cast3M


fr en

Télécharger ytoimp.eso

Retour à la liste

Numérotation des lignes : 

ytoimp

  1. C YTOIMP    SOURCE    GOUNAND   22/10/17    21:15:03     11483          
  2.       SUBROUTINE YTOIMP
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. C***********************************************************************
  6. C     -----------------------------------------------------------
  7. C     ---------  TOIMP   ----------------------------------------
  8. C     -----------------------------------------------------------
  9. C     ---------  PARAMETRI DELLO OPERATORE (NELLO ORDINE) : -----
  10. C     -----------------------------------------------------------
  11. C     ---------  TENSION ( tau et pression )            ---------
  12. C     -----------------------------------------------------------
  13. C
  14. C     SYNTAXE :
  15. C
  16. C       TOIMP  (tau pression)
  17. C
  18. C       1------2
  19. C       (R1,AL1)  LEF   FLUIDE  NOEUDS 1 2
  20. C
  21. C
  22. C
  23. C   CAS TRIDIMENSIONNEL
  24. C                       4 ________ 3
  25. C                       / FLUIDE /
  26. C                    1 /________/2
  27. C
  28. C
  29. C***********************************************************************
  30.  
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. -INC CCGEOME
  34. -INC CCREEL
  36. -INC SIZFFB
  37.       PARAMETER (LRV=64)
  38.       SEGMENT PETROV
  39.       REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG)
  40.       REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG)
  41.       REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV)
  42.       REAL*8 AJK(LRV,IDIM,IDIM,NPG)
  43.       ENDSEGMENT
  44.  
  45. -INC SMCOORD
  46. -INC SMLENTI
  47.       POINTEUR IPADU.MLENTI
  48. -INC SMELEME
  49.       POINTEUR MELEM1.MELEME
  50. -INC SMCHPOI
  51.       POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL
  52.       POINTEUR IZTU1.MPOVAL,MZMU.MPOVAL,MZUN.MPOVAL
  53.       POINTEUR MZTO.MPOVAL,MZDT.MPOVAL
  54.  
  55. -INC SMLMOTS
  56.       POINTEUR LINCO.MLMOTS
  57.       CHARACTER*8 TYPE,TYPC
  58.       CHARACTER*(LOCOMP)NOM0,NOMI,NOM4(3)
  59.       PARAMETER (NTB=1)
  60.       CHARACTER*8 LTAB(NTB)
  61.       DIMENSION KTAB(NTB),IXV(3)
  62.       SAVE IPAS
  63.       DATA LTAB/'KIZX    '/,IPAS/0/
  64. C*****************************************************************************
  65. CTOIMP
  66. C     write(6,*)' Debut TOIMP '
  67.       segact mcoord
  68.       CALL LITABS(LTAB,KTAB,NTB,1,IRET)
  69.       IF (IERR.NE.0) RETURN
  70.       MTABX=KTAB(1)
  71. C
  72. C- Récupération de la table EQEX (pointeur MTAB1)
  73. C
  74.       CALL LEKTAB(MTABX,'EQEX',MTAB1)
  75.       IF(MTAB1.EQ.0)THEN
  76. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  77.          MOTERR( 1: 8) = '  EQEX  '
  78.          MOTERR( 9:16) = '  EQEX  '
  79.          MOTERR(17:24) = '  KIZX  '
  80.          CALL ERREUR(786)
  81.       RETURN
  82.       ENDIF
  83.       CALL ACME(MTAB1,'NAVISTOK',NASTOK)
  84.       IF(NASTOK.EQ.0)THEN
  85.       CALL ZZOIMP(MTABX,MTAB1)
  86.       RETURN
  87.       ENDIF
  88. C
  89. C- Récupération de la table INCO (pointeur KINC)
  90. C
  91.       CALL LEKTAB(MTAB1,'INCO',KINC)
  92.       IF(KINC.EQ.0)THEN
  93. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  94.          MOTERR( 1: 8) = '  INCO  '
  95.          MOTERR( 9:16) = '  INCO  '
  96.          MOTERR(17:24) = '  EQEX  '
  97.          CALL ERREUR(786)
  98.       RETURN
  99.       ENDIF
  100.  
  101. C*****************************************************************************
  102. C OPTIONS
  103. C KIMPL = 0 -> EXPL     1 -> IMPL     2 -> SEMI
  104. C KFORM = 0 -> SI       1 -> EF       2 -> VF    3 -> EFMC
  105. C IDCEN = 0-> rien  1-> CENTREE  2-> SUPGDC  3-> SUPG  4-> TVISQUEU  5-> CNG
  106. C Pour IDCEN on autorise 1,2 et 3
  107.  
  108.       IAXI=0
  109.       IF(IFOMOD.EQ.0)IAXI=2
  110. C
  111. C- Récupération de la table des options KOPT (pointeur KOPTI)
  112. C
  113.       CALL LEKTAB(MTABX,'KOPT',KOPTI)
  114.       IF (KOPTI.EQ.0) THEN
  115. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  116.          MOTERR( 1: 8) = '  KOPT  '
  117.          MOTERR( 9:16) = '  KOPT  '
  118.          MOTERR(17:24) = '  KIZX  '
  119.          CALL ERREUR(786)
  120.          RETURN
  121.       ENDIF
  122.  
  123.       IKOMP=0
  124.       CALL ACME(KOPTI,'KIMPL',KIMPL)
  125.       CALL ACME(KOPTI,'KFORM',KFORM)
  126.       CALL ACMF(KOPTI,'CMD',CMD)
  127.       CALL ACME(KOPTI,'IDCEN',IDCEN)
  128.  
  129.       KDIM=1
  130.       IF(IDCEN.EQ.2)KDIM=IDIM
  131.       KDCEN=(IDCEN-1)*(IDCEN-2)*(IDCEN-3)
  132.       IF(KDCEN.NE.0)THEN
  133.             MOTERR( 1: 8) = 'IDCEN'
  134.             CALL ERREUR(803)
  135.             RETURN
  136.         ENDIF
  137.  
  138.       IF(KFORM.NE.0.AND.KFORM.NE.1)THEN
  139. C        Option %m1:8 incompatible avec les données
  140.             MOTERR( 1: 8) = 'EF/EFM1 '
  141.             CALL ERREUR(803)
  142.             RETURN
  143.         ENDIF
  144.       IF (IERR.NE.0) RETURN
  145.  
  146. C     write(6,*)' Apres les options '
  147. C*****************************************************************************
  148. C
  149. C- Récupération de la table DOMAINE associée au domaine local
  150. C
  151.       CALL LEKTAB(MTABX,'DOMZ',MTABZ)
  152.       IF(MTABZ.EQ.0)THEN
  153. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  154.          MOTERR( 1: 8) = '  DOMZ  '
  155.          MOTERR( 9:16) = '  DOMZ  '
  156.          MOTERR(17:24) = '  KIZX  '
  157.          CALL ERREUR(786)
  158.          RETURN
  159.       ENDIF
  160.  
  161.       CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  162.       CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  163.       IF (IERR.NE.0) RETURN
  164.  
  165.       SEGACT MELEME
  166.  
  167. C*************************************************************************
  168. C VERIFICATIONS SUR LES INCONNUES
  169. C
  170. C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI
  171. C
  172.       TYPE='LISTMOTS'
  173.       CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
  174.       IF (IERR.NE.0) RETURN
  175.       SEGACT LINCO
  176.       NBINC=LINCO.MOTS(/2)
  177.       IF(NBINC.NE.1)THEN
  178. C        Indice %m1:8 : contient plus de %i1 %m9:16
  179.          MOTERR( 1:8) = 'LISTINCO'
  180.          INTERR(1) = 1
  181.          MOTERR(9:16) = ' MOTS   '
  182.          CALL ERREUR(799)
  183.       RETURN
  184.       ENDIF
  185.  
  186. C --> 1 ere Inconnue
  187.  
  188.       NOMI=LINCO.MOTS(1)
  189.       DO 15 I=1,IDIM
  190.       WRITE(NOM4(I),FMT='(I1)')I
  191.       NOM4(I)=NOM4(I)(1:1)//NOMI(1:LOCOMP-1)
  192.  15   CONTINUE
  193.  
  194.       TYPE=' '
  195.       CALL ACMO(KINC,NOMI,TYPE,MCHPOI)
  196.       IF(TYPE.NE.'CHPOINT ')THEN
  197. C        Indice %m1:8 : ne contient pas un objet de type %m9:16
  198.          MOTERR='INC '//NOMI
  199.          MOTERR( 9:16) = 'CHPOINT '
  200.          CALL ERREUR(800)
  201.       RETURN
  202.       ELSE
  203.       CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1)
  204.       ENDIF
  205.       NPTU=IZTU1.VPOCHA(/2)
  206.  
  207. C*************************************************************************
  208. C Le domaine de definition est donne par le SPG  de la premiere inconnue
  209. C Les inconnues suivantes devront posseder ce meme pointeur
  210. C On verifie que les points de la zone sont tous inclus dans ce SPG
  211.  
  212.       CALL KRIPAD(MELEM1,MLENTI)
  213.       IPADU=MLENTI
  214.  
  215.       IF(IPAS.EQ.0)THEN
  216.       CALL VERPAD(MLENTI,MELEME,IRET)
  217.       IF(IRET.NE.0)THEN
  218. C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  219.          MOTERR= 'INC '//NOMI
  220.          MOTERR(9:16) = 'CHPOINT '
  221.          CALL ERREUR(788)
  222.       IPAS=0
  223.       RETURN
  224.       ENDIF
  225.       ENDIF
  226.  
  227. C*************************************************************************
  228. C Lecture du coefficient
  229. C Type du coefficient :
  230. C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur
  231.  
  232.       CALL ACME(MTABX,'IARG',IARG)
  233.       IF(IARG.NE.1.AND.IDCEN.EQ.1.AND.
  234.      &   IARG.NE.3.AND.(IDCEN.EQ.2.OR.IDCEN.EQ.2)
  235.      &)THEN
  236.       WRITE(6,*)'Opérateur TOIMP : nombre d''arguments incorrect'
  237. C           Indice %m1:8 : nombre d'arguments incorrect
  238.             MOTERR(1:8) = 'IARG    '
  239.             CALL ERREUR(804)
  240.       RETURN
  241.       ENDIF
  242.  
  243.       MZMU=IZTU1
  244.       MZUN=IZTU1
  245.       MZDT=IZTU1
  246.  
  247.       IXV(1)=-MELEMC
  248.       IXV(2)=0
  249.       IXV(3)=1
  250.       CALL LEKCOF('Opérateur TOIMP :',
  251.      & MTABX,KINC,1,IXV,MTO,MZTO,NTAU,NC1,IKS,IRET)
  252.       IF(IRET.EQ.0)RETURN
  253.  
  254.       IF(IKS.EQ.2)IKS=1
  255.  
  256.  
  257.       IF(IARG.GE.2)THEN
  258.  
  259. C 2ème coefficient  UN , champ de vitesse transportant
  260.       IXV(1)=-MELEM1
  261.       IXV(2)=0
  262.       IXV(3)=1
  263.       CALL LEKCOF('Opérateur TOIM :',
  264.      & MTABX,KINC,2,IXV,MUN,MZUN,NPTU,NC2,IKU,IRET)
  265.       IF(IRET.EQ.0)GO TO 90
  266.       IF(IKU.EQ.2)IKU=1
  267.  
  268. C 3ème coefficient  viscosité
  269.       IXV(1)=0
  270.       IXV(2)=1
  271.       IXV(3)=0
  272.       CALL LEKCOF('Opérateur TOIM :',
  273.      & MTABX,KINC,3,IXV,MMU,MZMU,NPT3,NC3,IKM,IRET)
  274.       IF(IRET.EQ.0)GO TO 90
  275.       ELSE
  276.       CALL LEKTAB(MTABX,'ARGS3',MMU)
  277.       CALL LICHT(MMU,MZMU,TYPC,IGEOM0)
  278.       MZMU.VPOCHA(1,1)=1.D0
  279.       IKM=1
  280.  
  281.       ENDIF
  282.  
  283.       IF(IARG.EQ.4)THEN
  284. C 4ème coefficient  Dt
  285.       IXV(1)=0
  286.       IXV(2)=1
  287.       IXV(3)=0
  288.       CALL LEKCOF('Opérateur TOIM :',
  289.      & MTABX,KINC,4,IXV,MDT,MZDT,NPT4,NC4,IKT,IRET)
  290.       IF(IRET.EQ.0)RETURN
  291.       DT=MZDT.VPOCHA(1,1)
  292.       ELSE
  293.       DT=0.D0
  294.  
  295.       ENDIF
  296.  
  297. C     write(6,*)' Fin lecture Arguments '
  298. C Fin lecture Arguments ************************************************
  299.  
  300.  
  301.       IF(KIMPL.EQ.0)THEN
  302.       NC=IZTU1.VPOCHA(/2)
  303.       TYPE='SOMMET'
  304.       CALL KRCHPT(TYPE,MELEM1,NC,IZG1,NOM4)
  305.  
  306.       ELSE
  307.  
  308.       NAT=2
  309.       NSOUPO=1
  310.       SEGACT MELEM1
  311.       N=MELEM1.NUM(/2)
  312.       NC=IZTU1.VPOCHA(/2)
  313.       NINKO=NC
  314.       SEGINI MCHPO1,MSOUP1,MPOVA1
  315.       MCHPO1.IFOPOI=IFOUR
  316.       MCHPO1.MOCHDE=TITREE
  317.       MCHPO1.MTYPOI='SMBR'
  318.       MCHPO1.JATTRI(1)=2
  319.       MCHPO1.IPCHP(1)=MSOUP1
  320.       DO 177 N=1,NINKO
  321.       WRITE(NOM0,FMT='(I1)')N
  322.       MSOUP1.NOCOMP(N)=NOM0(1:1)//NOMI(1:LOCOMP-1)
  323.  177  CONTINUE
  324.       MSOUP1.IGEOC=MELEM1
  325.       MSOUP1.IPOVAL=MPOVA1
  326.       IZG1=MCHPO1
  327.  
  328.       ENDIF
  329.  
  330.       CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)
  331.  
  332.       IF(IGEOM.NE.MELEM1)THEN
  333.       WRITE(6,*)' Opérateur TOIM'
  334.       WRITE(6,*)' Incompatibilité de SPG entre 1ères inconnues'
  335.       RETURN
  336.       ENDIF
  337.  
  338.       SEGACT MELEME
  339.       NBSOUS=LISOUS(/1)
  340.       IF(NBSOUS.EQ.0)NBSOUS=1
  341.  
  342.       NPTD=IZGG1.VPOCHA(/1)
  343.       K0=0
  344.       DO 1 LS=1,NBSOUS
  345.       IPT1=MELEME
  346.       IF(NBSOUS.NE.1)IPT1=LISOUS(LS)
  347.       SEGACT IPT1
  348.  
  349.       NOM0 = NOMS(IPT1.ITYPEL)//'    '
  350.       CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  351.       SEGACT IZFFM
  352.       IZHR=KZHR(1)
  353.       SEGACT IZHR*MOD
  354.       NES=GR(/1)
  355.       NPG=GR(/3)
  356.  
  357.       NP  =IPT1.NUM(/1)
  358.       NBEL=IPT1.NUM(/2)
  359.  
  360.       SEGINI PETROV
  361.  
  362.       DEUPI=1.D0
  363.       IF(IAXI.NE.0)DEUPI=2.D0*XPI
  364.  
  365. C Calcul du nombre de paquets de LRV éléments
  366. C
  367.       NNN=MOD(NBEL,LRV)
  368.       IF(NNN.EQ.0) THEN
  369.          NPACK=NBEL/LRV
  370.       ELSE
  371.          NPACK=1+(NBEL-NNN)/LRV
  372.       ENDIF
  373.       KPACKD=1
  374.       KPACKF=NPACK
  375. C
  376. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  377. C
  378.       DO 7001 KPACK=KPACKD,KPACKF
  379. C
  380. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  381. C
  382. C 1. Calcul des limites du paquet courant.
  383.          KDEB=1+(KPACK-1)*LRV
  384.          KFIN=MIN(NBEL,KDEB+LRV-1)
  385. C
  386.          DO 7002 K=KDEB,KFIN
  387.             KP=K-KDEB+1
  388.             NK=K0+K
  389.             NKM=(1-IKM)*(NK-1)+1
  390.             ANUK(KP)=MZMU.VPOCHA(NKM,1)
  391.  7002    CONTINUE
  392.  
  393. C     CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN
  394.          DT=0.D0
  395.          IF(IDCEN.EQ.2)THEN
  396.             DO 108 NC=1,IDIM
  397.                CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  398.      &              NES,NP,NPG,IAXI,XCOOR,
  399.      &              IPT1.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  400.      &              MZUN.VPOCHA(1,NC),IPADU.LECT,MZUN.VPOCHA,IPADU.LECT
  401.      $              ,NPTU,ANUK,WT(1,1,1,NC),WS(1,1,1,NC),HK,PGSK,RPGK
  402.      $              ,AJK,AIRE,UIL,DUIL,DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  403.  108        CONTINUE
  404.          ELSE
  405.             CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  406.      &           NES,NP,NPG,IAXI,XCOOR,
  407.      &           IPT1.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  408.      &           MZUN.VPOCHA,IPADU.LECT,MZUN.VPOCHA,IPADU.LECT,NPTU,ANUK
  409.      $           ,WT,WS,HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL,DTM1,DT,DTT1,DTT2
  410.      $           ,DIAEL,NUEL)
  411.          ENDIF
  412.  
  413.          IF(NES.NE.IDIM)THEN
  414.             DO 7003 K=KDEB,KFIN
  415.                KP=K-KDEB+1
  416.                NK=K0+K
  417.                KA=1+(1-IKS)*(NK-1)
  418.                DO    I=1,NP
  419.                   DO    N=1,IDIM
  420.                      N1=1
  421.                      IF(IDCEN.EQ.2) N1=N
  422.                      FF1=0.D0
  423.                      DO 51 L=1,NPG
  424.                         DO 52 M=1,IDIM
  425.                            FF1=FF1 + MZTO.VPOCHA(KA,M)*AJK(KP,N,M,L)
  426.      $                          *WT(KP,I,L,N1)*PGSK(KP,L)*DEUPI*RPGK(KP
  427.      $                          ,L)
  428.  52                     CONTINUE
  429.  51                  CONTINUE
  430.                      NF=LECT(IPT1.NUM(I,K))
  431.                      IZGG1.VPOCHA(NF,N)=IZGG1.VPOCHA(NF,N)+FF1
  432.                   ENDDO
  433.                ENDDO
  434.  7003       CONTINUE
  435.          ELSEIF(NES.EQ.IDIM)THEN
  436.             DO 7004 K=KDEB,KFIN
  437.                KP=K-KDEB+1
  438.                NK=K0+K
  439.                KA=1+(1-IKS)*(NK-1)
  440.                DO    I=1,NP
  441.                   DO    N=1,IDIM
  442.                      N1=1
  443.                      IF(IDCEN.EQ.2)N1=N
  444.                      FF1=0.D0
  445.                      DO 61 L=1,NPG
  446.                         FF1=FF1 + MZTO.VPOCHA(KA,N)*WT(KP,I,L,N1)
  447.      $                       *PGSK(KP,L)*DEUPI*RPGK(KP,L)
  448.  61                  CONTINUE
  449.                      NF=LECT(IPT1.NUM(I,K))
  450.                      IZGG1.VPOCHA(NF,N)=IZGG1.VPOCHA(NF,N)+FF1
  451.                   ENDDO
  452.                ENDDO
  453.  7004       CONTINUE
  454.          ENDIF
  455.  7001 CONTINUE
  456.       SEGDES IPT1
  457.       K0=K0+NBEL
  458.       SEGSUP PETROV
  459.  1    CONTINUE
  460.       SEGDES MZTO
  461.       SEGDES MELEME
  462.  
  463.       NRIGE=7
  464.       NKID =9
  465.       NKMT =7
  466.       NMATRI=0
  467.       SEGINI MATRIK
  468.       CALL ECROBJ('MATRIK',MATRIK)
  469.       CALL ECROBJ('CHPOINT',IZG1)
  470.       SEGDES MELEME
  471.       SEGDES IZTU1
  472.       SEGDES IZG1,IZGG1
  473.       SEGDES LINCO
  474.       SEGSUP MLENTI
  475.       IPAS=1
  476. C     write(6,*)' FIN TOIMP '
  477.       RETURN
  478.  90   CONTINUE
  479.       WRITE(6,*)' Interuption anormale de TOIMP '
  480. C        Option %m1:8 incompatible avec les données
  481.             MOTERR( 1: 8) = '  EF    '
  482.       CALL ERREUR(803)
  483.       RETURN
  484.  1002 FORMAT(10(1X,1PE11.4))
  485.       END
  486.  
  487.  

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