Sources Esope | Cast3M

pre112
C PRE112    SOURCE    OF166741  24/12/13    21:16:57     12097          
C PRE111    SOURCE    BECC      98/08/24    21:15:02     3286
      SUBROUTINE PRE112(ICEN,IFACE,IFACEL,INORM,
     &                  IROC,IVITC,IPC,IGAMC,
     &                  IROF,IVITF,IPF,IGAMF,
     &                  LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PRE111
C
C DESCRIPTION       :  Voir PRE11
C
C                      Cas Deux Dimensions
C
C                          Une Espece
C
C                          1er ordre en espace, 1re ordre en temps
C
C                      Creations des objets MCHAML IROF, IVITF, IPF,
C                      IGAMF
C
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
C
C************************************************************************
C
C
C APPELES (Outils)  :  KRIPAD, LICHT
C
C APPELES (Calcul)  :  AUCUN
C
C
C************************************************************************
C
C ENTREES
C
C     1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
C
C     ICEN    : MELEME de 'POI1' SPG des CENTRES
C
C     IFACE   : MELEME de 'POI1' SPG des FACES
C
C     IFACEL  : MELEME de 'SEG3' avec
C               CENTRE d'Elt "gauche"
C               CENTRE de Face
C               CENTRE d'Elt "droite"
C
C               N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
C
C     INORM   : CHPOINT des cosinus directeurs de normales aux faces
C
C     2) Pointeurs des  CHPOINTs
C
C     IROC    : CHPOINT "CENTRE" contenant la masse volumique RHO
C
C     IVITC   : CHPOINT "CENTRE" contenant la vitesse UX, UY ;
C
C     IPC     : CHPOINT "CENTRE" contenat la pression P;
C
C     IGAMC   : CHPOINT "CENTRE" contenat le "Gamma" du gaz
C
C
C SORTIES
C
C
C     IROF    : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la masse volumique RHO
C
C     IVITF   : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la vitesse UN, UT dans le repaire local
C               (n,t) et defini sur le MELEME de pointeur IFACE,
C               contenant les cosinus directeurs du repere local
C
C     IPF     : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la pression P
C
C     IGAMF   : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant le "gamma" du gaz
C
C     LOGAN   : anomalie detectee (changement de la convention dans
C               la table domaine)
C
C     LOGNEG  : (LOGICAL): si .TRUE. une pression ou une densité
C               negative a été detectée -> en interactif le
C               programme s'arrete  en GIBIANE
C               (erreur stocké en MESERR et VALER)
C
C     LOGBOR  : (LOGICAL): si .TRUE. un gamma a ete detecte
C               dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
C               en VAL1 et en VAL2 on stocke 1.0 et 3.0)
C
C     MESERR
C     VALER
C     VAL1,
C     VAL2    : pour les messages d'erreur
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créée le 11.6.98.
C
C************************************************************************
C
C
C ATTENTION: Cet programme marche que si le MAILLAGE est convex;
C            si non il faut changer l'argoritme de calcul de
C            l'orientation des normales aux faces.
C
C
C************************************************************************
C
C**** Les variables
C
      IMPLICIT INTEGER(I-N)
      INTEGER ICEN, IFACE, IFACEL, IROC, IVITC, IPC ,IGAMC, INORM
     &       , IROF, IVITF, IPF, IGAMF
     &       , IGEOM, NFAC
     &       , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1
     &       , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1
      REAL*8 VALER, VAL1, VAL2,XG,YG,ZG,XC,YC,ZC
     &       ,DXG, DYG, DZG,ORIENT
     &       , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
     &       , ROG, PG, GAMG, UXG, UYG, UZG, UNG, UTG, UVG
     &       , ROD, PD, GAMD, UXD, UYD, UZD, UND, UTD, UVD
      CHARACTER*(40) MESERR
      CHARACTER*(8) TYPE
      LOGICAL LOGAN,LOGNEG, LOGBOR
C
C**** Les Includes
C
-INC SMCOORD
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
      POINTEUR MPROC.MPOVAL, MPVITC.MPOVAL, MPPC.MPOVAL,
     &         MPGAMC.MPOVAL, MPNORM.MPOVAL
-INC SMCHAML
      POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
     &         MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
     &         MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
      POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL
      POINTEUR MELRO.MELVAL, MELP.MELVAL,
     &         MELGAM.MELVAL
-INC SMLENTI
-INC SMELEME
C
C**** Initialisation des parametres d'erreur déjà faite, i.e.
C
C     LOGNEG  = .FALSE.
C     LOGBOR  = .FALSE.
C     MESERR = '                                        '
C     MOTERR(1:40) = MESERR(1:40)
C     VALER = 0.0D0
C     VAL1 = 0.0D0
C     VAL2 = 0.0D0
C
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(ICEN,MLENT1)
C
C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
C
C     Si i est le numero global d'un noeud de ICEN,
C     MLENT1.LECT(i) contient sa position, i.e.
C
C     I              = numero global du noeud centre
C     MLENT1.LECT(i) = numero local du noeud centre
C
C     MLENT1 déjà activé, i.e.
C
C     SEGACT MLENT1
C
C**** Activation de CHPOINTs
C
C     densité
C     vitesse
C     pression
C     gamma
C     cosinus directeurs des normales aux surface
C
      CALL LICHT(IROC ,MPROC ,TYPE,IGEOM)
      CALL LICHT(IVITC,MPVITC,TYPE,IGEOM)
      CALL LICHT(IPC  ,MPPC  ,TYPE,IGEOM)
      CALL LICHT(IGAMC,MPGAMC,TYPE,IGEOM)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
C
C**** MPOVA1 - MPOVA5 sont déjà activés i.e.:
C
C     SEGACT MPROC
C     SEGACT MPVITC
C     SEGACT MPPC
C     SEGACT MPGAMC
C     SEGACT MPNORM
C
C**** Le MELEME FACEL
C
      IPT1 = IFACEL
      IPT2 = IFACE
      SEGACT IPT1
      SEGACT IPT2
      NFAC = IPT1.NUM(/2)
C
C**** Creation de MCHAMLs contenant les etats gauche et droite,
C
C     i.e.:
C
C     vitesse + cosinus directors du repere local
C     densité
C     pression
C     gamma
C
C**** Cosinus directors du repere local et vitesse
C
C     Les cosinus directeurs
C
      N1 = 2
      N3 = 6
      L1 = 28
      SEGINI MCHEL1
      IVITF = MCHEL1
      MCHEL1.TITCHE = 'U  '
      MCHEL1.IMACHE(1) = IFACE
      MCHEL1.IMACHE(2) = IFACEL
      MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)'
      MCHEL1.CONCHE(2) = ' U in (n,t,v)   '
C
C**** Valeurs des cosinus definies par respect au repair global, i.e.
C
      MCHEL1.INFCHE(1,1) = 2
      MCHEL1.INFCHE(1,3) = NIFOUR
      MCHEL1.INFCHE(1,4) = 0
      MCHEL1.INFCHE(1,5) = 0
      MCHEL1.INFCHE(1,6) = 0
      MCHEL1.IFOCHE = IFOUR
C
C**** Valeurs de vitesse definies par respect au repair local, i.e.
C
      MCHEL1.INFCHE(2,1) = 1
      MCHEL1.INFCHE(2,3) = NIFOUR
      MCHEL1.INFCHE(2,4) = 0
      MCHEL1.INFCHE(2,5) = 0
      MCHEL1.INFCHE(2,6) = 0
C
C**** Le cosinus directeurs
C
      N1PTEL = 1
      N1EL = NFAC
      N2PTEL = 0
      N2EL =   0
C
C**** MCHAML a N2 composantes:
C
C     cosinus directeurs du repere local (n,t1,t2)=(n,t,v)
C
C     IDIM = 3 -> 9 composantes
C
      N2 = 9
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NX      '
      MCHAM1.NOMCHE(2) = 'NY      '
      MCHAM1.NOMCHE(3) = 'NZ      '
      MCHAM1.NOMCHE(4) = 'TX      '
      MCHAM1.NOMCHE(5) = 'TY      '
      MCHAM1.NOMCHE(6) = 'TZ      '
      MCHAM1.NOMCHE(7) = 'VX      '
      MCHAM1.NOMCHE(8) = 'VY      '
      MCHAM1.NOMCHE(9) = 'VZ      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      MCHAM1.TYPCHE(5) = 'REAL*8          '
      MCHAM1.TYPCHE(6) = 'REAL*8          '
      MCHAM1.TYPCHE(7) = 'REAL*8          '
      MCHAM1.TYPCHE(8) = 'REAL*8          '
      MCHAM1.TYPCHE(9) = 'REAL*8          '
      SEGINI MELVNX
      SEGINI MELVNY
      SEGINI MELVNZ
      SEGINI MELT1X
      SEGINI MELT1Y
      SEGINI MELT1Z
      SEGINI MELT2X
      SEGINI MELT2Y
      SEGINI MELT2Z
      MCHAM1.IELVAL(1) = MELVNX
      MCHAM1.IELVAL(2) = MELVNY
      MCHAM1.IELVAL(3) = MELVNZ
      MCHAM1.IELVAL(4) = MELT1X
      MCHAM1.IELVAL(5) = MELT1Y
      MCHAM1.IELVAL(6) = MELT1Z
      MCHAM1.IELVAL(7) = MELT2X
      MCHAM1.IELVAL(8) = MELT2Y
      MCHAM1.IELVAL(9) = MELT2Z
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 3 -> 3 composantes
C
      N2 = 3
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'UN      '
      MCHAM1.NOMCHE(2) = 'UT      '
      MCHAM1.NOMCHE(3) = 'UV      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      SEGINI MELVUN
      SEGINI MELVUT
      SEGINI MELVUV
      MCHAM1.IELVAL(1) = MELVUN
      MCHAM1.IELVAL(2) = MELVUT
      MCHAM1.IELVAL(3) = MELVUV
      SEGDES MCHAM1
C
C**** Densite
C
      N1 = 1
      N3 = 6
      L1 = 15
      SEGINI MCHEL2
      IROF = MCHEL2
      MCHEL2.IMACHE(1) = IFACEL
      MCHEL2.TITCHE = 'RO             '
      MCHEL2.CONCHE(1) = '                '
C
C**** Valeurs independente du repére, i.e.
C
      MCHEL2.INFCHE(1,1) = 0
      MCHEL2.INFCHE(1,3) = NIFOUR
      MCHEL2.INFCHE(1,4) = 0
      MCHEL2.INFCHE(1,5) = 0
      MCHEL2.INFCHE(1,6) = 0
      MCHEL2.IFOCHE = IFOUR
      N2 = 1
      SEGINI MCHAM1
      MCHEL2.ICHAML(1) = MCHAM1
      SEGDES MCHEL2
      MCHAM1.NOMCHE(1) = 'SCAL    '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      SEGINI MELRO
      MCHAM1.IELVAL(1) = MELRO
      SEGDES MCHAM1
C
C**** Pression
C
      MCHEL1 = IROF
      SEGINI, MCHEL2 = MCHEL1
      IPF = MCHEL2
      MCHEL2.TITCHE = 'P            '
C
C**** MCHAM1 = MCHAML de la densite
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MELP
      MCHAM2.IELVAL(1) = MELP
      SEGDES MCHAM2
C
C**** Gamma
C
      MCHEL1 = IROF
      SEGINI, MCHEL2 = MCHEL1
      IGAMF = MCHEL2
      MCHEL2.TITCHE = 'GAMMA        '
C
C**** MCHAM1 = MCHAML de la densite
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MELGAM
      MCHAM2.IELVAL(1) = MELGAM
      SEGDES MCHAM2
C
C**** Recapitulatif: les MELVALs et les MPOVALs actives
C
C     MELVNX, MELVNY, MELVNZ
C     MELT1X, MELT1Y, MELT1Z
C     MELT2X, MELT2Y, MELT2Z
C
C     MELVUN, MELVUT , MELVUN-> vitesse
C
C     MELRO  -> densite
C
C     MELP   -> pression
C
C     MELGAM -> gamma
C
C     MPROC  -> densite
C
C     MPVITC -> vitesse
C
C     MPPC   -> pression
C
C     MPGAMC -> gamma
C
C     MPNORM -> normales aux faces
C
C**** Boucle sur le faces
C
      IDIMP1=IDIM+1
      DO NLCF = 1, NFAC
C
C******* NLCF  = numero local du centre de face
C        NGCF  = numero global du centre de face
C        NGCEG = numero global du centre ELT "gauche"
C        NLCEG = numero local du centre ELT "gauche"
C        NGCED = numero global du centre ELT "droite"
C        NLCED = numero local du centre ELT "droite"
C
         NGCEG = IPT1.NUM(1,NLCF)
         NGCF  = IPT1.NUM(2,NLCF)
         NGCED = IPT1.NUM(3,NLCF)
         NLCEG = MLENT1.LECT(NGCEG)
         NLCED = MLENT1.LECT(NGCED)
C
C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
C
         NGCF1 = IPT2.NUM(1,NLCF)
         IF( NGCF1 .NE. NGCF) THEN
            LOGAN = .TRUE.
            MESERR(1:40) = 'PRET, subroutine pre112.eso             '
            GOTO 9999
         ENDIF
C
C******* Cosinus directeurs des NORMALES aux faces
C
C        On impose que les normales sont direct "Gauche" -> "Centre"
C
         INDCEL = (NGCEG-1)*IDIMP1
         XG = XCOOR(INDCEL+1)
         YG = XCOOR(INDCEL+2)
         ZG = XCOOR(INDCEL+3)
         INDCEL = (NGCF-1)*IDIMP1
         XC = XCOOR(INDCEL + 1)
         YC = XCOOR(INDCEL + 2)
         ZC = XCOOR(INDCEL+3)
         DXG = XC - XG
         DYG = YC - YG
         DZG = ZC - ZG
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
         CNX = MPNORM.VPOCHA(NLCF,7)
         CNY = MPNORM.VPOCHA(NLCF,8)
         CNZ = MPNORM.VPOCHA(NLCF,9)
         ORIENT = CNX * DXG + CNY * DYG + CNZ * DZG
         ORIENT = SIGN(1.0D0,ORIENT)
         IF(ORIENT .NE. 1.0D0)THEN
            LOGAN = .TRUE.
            MESERR(1:30)=
     &           'PRET , subroutine pre121.eso. '
            GOTO 9999
         ENDIF
         CNX = CNX * ORIENT
         CNY = CNY * ORIENT
         CNZ = CNZ * ORIENT
C
C********** Cosinus directeurs de tangente 1
C
         CTX = MPNORM.VPOCHA(NLCF,1) * ORIENT
         CTY = MPNORM.VPOCHA(NLCF,2) * ORIENT
         CTZ = MPNORM.VPOCHA(NLCF,3) * ORIENT
C
C********** Cosinus directeurs de tangente 2
C
         CVX = MPNORM.VPOCHA(NLCF,4) * ORIENT
         CVY = MPNORM.VPOCHA(NLCF,5) * ORIENT
         CVZ = MPNORM.VPOCHA(NLCF,6) * ORIENT
C
C******* Les autres MELVALs
C
C
C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
C              GAMMA \in (1,3)
C              Si non il faut le faire, en utlisant LOGBOR,
C              LOGNEG, VALER, VAL1, VAL2
C
C
C
C******* NGCEG = NGCED -> Mur
C
         IF( NGCEG .EQ. NGCED)THEN
            ROG = MPROC.VPOCHA(NLCEG , 1)
            PG = MPPC.VPOCHA(NLCEG, 1)
            GAMG = MPGAMC.VPOCHA(NLCEG, 1)
            UXG = MPVITC.VPOCHA(NLCEG , 1)
            UYG = MPVITC.VPOCHA(NLCEG , 2)
            UZG = MPVITC.VPOCHA(NLCEG , 3)
            UNG = UXG * CNX + UYG * CNY + UZG * CNZ
            UTG = UXG * CTX + UYG * CTY + UZG * CTZ
            UVG = UXG * CVX + UYG * CVY + UZG * CVZ
C
C********** Son etat droite
C
            ROD = ROG
            PD = PG
            GAMD = GAMG
            UND = -1.0D0 * UNG
            UTD = UTG
            UVD = UVG
C
C************* Fin cas mur
C
         ELSE
C
C************* Etat gauche
C
            ROG = MPROC.VPOCHA(NLCEG, 1)
            PG = MPPC.VPOCHA(NLCEG, 1)
            GAMG = MPGAMC.VPOCHA(NLCEG, 1)
            UXG = MPVITC.VPOCHA(NLCEG , 1)
            UYG = MPVITC.VPOCHA(NLCEG , 2)
            UZG = MPVITC.VPOCHA(NLCEG , 3)
            UNG = UXG * CNX + UYG * CNY + UZG * CNZ
            UTG = UXG * CTX + UYG * CTY + UZG * CTZ
            UVG = UXG * CVX + UYG * CVY + UZG * CVZ
C
C********** Etat gauche
C
            ROD = MPROC.VPOCHA(NLCED,1)
            PD = MPPC.VPOCHA(NLCED,1)
            GAMD = MPGAMC.VPOCHA(NLCED,1)
C
C************* On suppose qu'on a déjà controlle ROG, PG > 0
C              Si non il faut le faire!!!
C
            UXD = MPVITC.VPOCHA(NLCED,1)
            UYD = MPVITC.VPOCHA(NLCED,2)
            UZD = MPVITC.VPOCHA(NLCED,3)
            UND = UXD * CNX + UYD * CNY + UZD * CNZ
            UTD = UXD * CTX + UYD * CTY + UZD * CTZ
            UVD = UXD * CVX + UYD * CVY + UZD * CVZ
         ENDIF
C
C************* Les MELVALs
C
         MELRO.VELCHE(1,NLCF) = ROG
         MELRO.VELCHE(3,NLCF) = ROD
         MELP.VELCHE(1,NLCF) = PG
         MELP.VELCHE(3,NLCF) = PD
         MELGAM.VELCHE(1,NLCF) = GAMG
         MELGAM.VELCHE(3,NLCF) = GAMD
         MELVUN.VELCHE(1,NLCF) = UNG
         MELVUN.VELCHE(3,NLCF) = UND
         MELVUT.VELCHE(1,NLCF) = UTG
         MELVUT.VELCHE(3,NLCF) = UTD
         MELVUV.VELCHE(1,NLCF) = UVG
         MELVUV.VELCHE(3,NLCF) = UVD
         MELVNX.VELCHE(1,NLCF) = CNX
         MELVNY.VELCHE(1,NLCF) = CNY
         MELVNZ.VELCHE(1,NLCF) = CNZ
         MELT1X.VELCHE(1,NLCF) = CTX
         MELT1Y.VELCHE(1,NLCF) = CTY
         MELT1Z.VELCHE(1,NLCF) = CTZ
         MELT2X.VELCHE(1,NLCF) = CVX
         MELT2Y.VELCHE(1,NLCF) = CVY
         MELT2Z.VELCHE(1,NLCF) = CVZ
      ENDDO
C
C**** Desactivation des SEGMENTs
C
      SEGDES IPT1
      SEGDES IPT2
C
      SEGDES MPROC
      SEGDES MPVITC
      SEGDES MPPC
      SEGDES MPGAMC
      SEGDES MPNORM
C
      SEGDES MELRO
      SEGDES MELP
      SEGDES MELGAM
      SEGDES MELVUN
      SEGDES MELVUT
      SEGDES MELVUV
      SEGDES MELVNX
      SEGDES MELVNY
      SEGDES MELVNZ
      SEGDES MELT1X
      SEGDES MELT1Y
      SEGDES MELT1Z
      SEGDES MELT2X
      SEGDES MELT2Y
      SEGDES MELT2Z
C
C**** Destruction du MELNTI correspondance local/global
C
      SEGSUP MLENT1
C
 9999 CONTINUE
C
      RETURN
      END