Sources Esope | Cast3M

pre311
C PRE311    SOURCE    OF166741  24/10/03    21:15:33     12022          
      SUBROUTINE PRE311(ICEN,IFACE,IFACEL,INORM,IROC,IVITC,IPC,
     &                  IYC,ISCAC,
     &                  IROF,IVITF,IPF,IYF,ISCAF,
     &                  LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PRE311
C
C DESCRIPTION       :  Voir PRE31
C
C                      Cas Deux Dimensions
C
C                          Mono/MultiEspeces
C
C                          1er ordre en espace, 1re ordre en temps
C
C                      Creations des objets MCHAML IROF, IVITF, IPF,IYF
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     IYC     : CHPOINT "CENTRE" contenat les fractions massiques
C               (ou 0 dans le cas monoespece);
C
C     ISCAC   : CHPOINT "CENTRE" contenat les scalaires passifs
C               (ou 0);
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     IYF     : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant les fractions massiques (ou 0 dans le cas
C               monoespece);
C
C     ISCAF   : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant les scalaire passifs (ou 0)
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 Y 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 18.12.98.
C
C
C               17.02.2000: transport des scalaires passifs
C
C************************************************************************
C
C
C ATTENTION: Cet programme marche 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 ,IYC, ISCAC, INORM
     &       , IROF, IVITF, IPF, IYF, NESP, ISCAF, NSCA
     &       , IGEOM, NFAC
     &       , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1
     &       , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1, I1, NMA
      REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, DXG, DYG
     &       , CNX, CNY, CTX, CTY, ORIENT
     &       , ROG, PG, UXG, UYG, UNG, UTG
     &       , ROD, PD, UXD, UYD, UND, UTD
      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,
     &         MPNORM.MPOVAL, MPYC.MPOVAL, MPSC.MPOVAL
-INC SMCHAML
      POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
     &         MELT1X.MELVAL, MELT1Y.MELVAL
      POINTEUR MELVUN.MELVAL, MELVUT.MELVAL
      POINTEUR MELRO.MELVAL, MELP.MELVAL
      POINTEUR MCHAMY.MCHAML, MCHAMS.MCHAML
-INC SMLENTI
-INC SMELEME
C
C**** Segments des fractions massiques gauche et droit
C
      SEGMENT FRAMAS
         REAL*8 FRAMG(NMA), FRAMD(NMA)
      ENDSEGMENT
      POINTEUR SCALPA.FRAMAS
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     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(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 MPNORM
C
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     fractions massiques
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) in (x,y) '
      MCHEL1.CONCHE(2) = ' U in (n,t)     '
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) = 1
      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) = 1
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)
C
C     IDIM = 2 -> 4 composantes
C
      N2 = 4
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NX      '
      MCHAM1.NOMCHE(2) = 'NY      '
      MCHAM1.NOMCHE(3) = 'TX      '
      MCHAM1.NOMCHE(4) = 'TY      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      SEGINI MELVNX
      SEGINI MELVNY
      SEGINI MELT1X
      SEGINI MELT1Y
      MCHAM1.IELVAL(1) = MELVNX
      MCHAM1.IELVAL(2) = MELVNY
      MCHAM1.IELVAL(3) = MELT1X
      MCHAM1.IELVAL(4) = MELT1Y
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 2 -> 2 composantes
C
      N2 = 2
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'UN      '
      MCHAM1.NOMCHE(2) = 'UT      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      SEGINI MELVUN
      SEGINI MELVUT
      MCHAM1.IELVAL(1) = MELVUN
      MCHAM1.IELVAL(2) = MELVUT
      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) = 1
      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**** Les fractions massiques: le CHPOINT et le relative CHAMELEM
C
      IF(IYC .NE. 0)THEN
         MCHPO1 = IYC
         SEGACT MCHPO1
         MSOUP1 = MCHPO1.IPCHP(1)
         SEGDES MCHPO1
         SEGACT MSOUP1
         NESP = MSOUP1.NOCOMP(/2)
         MPYC = MSOUP1.IPOVAL
         SEGACT MPYC
C
         MCHEL1 = IROF
         SEGINI, MCHEL2 = MCHEL1
         IYF = MCHEL2
         MCHEL2.TITCHE = 'Y            '
         N2 = NESP
         SEGINI MCHAMY
         MCHEL2.ICHAML(1) = MCHAMY
         SEGDES MCHEL2
         N1EL = NFAC
         N1PTEL = 3
         N2EL = 0
         N2PTEL = 0
         DO I1 = 1, NESP
            SEGINI MELVA1
            MCHAMY.IELVAL(I1) = MELVA1
            MCHAMY.NOMCHE(I1) = MSOUP1.NOCOMP(I1)
            MCHAMY.TYPCHE(I1) = 'REAL*8          '
         ENDDO
C
         SEGDES MSOUP1
         NMA = NESP
         SEGINI FRAMAS
      ELSE
         IYF = 0
         NESP = 0
      ENDIF
C
C**** Les scalaires passifs: le CHPOINT et le relative CHAMELEM
C
      IF(ISCAC .NE. 0)THEN
         MCHPO1 = ISCAC
         SEGACT MCHPO1
         MSOUP1 = MCHPO1.IPCHP(1)
         SEGDES MCHPO1
         SEGACT MSOUP1
         NSCA = MSOUP1.NOCOMP(/2)
         MPSC = MSOUP1.IPOVAL
         SEGACT MPSC
C
         MCHEL1 = IROF
         SEGINI, MCHEL2 = MCHEL1
         ISCAF = MCHEL2
         MCHEL2.TITCHE = 'SCALPASS     '
         N2 = NSCA
         SEGINI MCHAMS
         MCHEL2.ICHAML(1) = MCHAMS
         SEGDES MCHEL2
         N1EL = NFAC
         N1PTEL = 3
         N2EL = 0
         N2PTEL = 0
         DO I1 = 1, NSCA
            SEGINI MELVA1
            MCHAMS.IELVAL(I1) = MELVA1
            MCHAMS.NOMCHE(I1) = MSOUP1.NOCOMP(I1)
            MCHAMS.TYPCHE(I1) = 'REAL*8          '
         ENDDO
C
         SEGDES MSOUP1
         NMA = NSCA
         SEGINI SCALPA
C
C**** On laisse actives les segments pointes par
C     MPYC, MCHAMY,FRAMAS, et le MELVALs relatifs aux
C     fractions massiques
C
C
      ELSE
         ISCAF = 0
         NSCA = 0
      ENDIF
C
C**** Recapitulatif
C
C     MELVNX, MELVNY
C     MELT1X, MELT1Y -> normales et tangentes aux faces
C
C     MELVUN, MELVUT -> vitesse
C
C     MELRO  -> densite
C
C     MELP   -> pression
C
C     MPROC  -> densite
C
C     MPVITC -> vitesse
C
C     MPPC   -> pression
C
C     MPNORM -> normales aux faces
C
C**** Boucle sur le faces
C
      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 pre311.eso             '
            GOTO 9999
         ENDIF
C
C******* Cosinus directeurs des NORMALES aux faces
C
C        On impose que les normales sont direct "Gauche" -> "Centre"
C
         XG = XCOOR((NGCEG-1)*(IDIM+1)+1)
         YG = XCOOR((NGCEG-1)*(IDIM+1)+2)
         XC = XCOOR((NGCF-1)*(IDIM+1)+1)
         YC = XCOOR((NGCF-1)*(IDIM+1)+2)
         DXG = XC - XG
         DYG = YC - YG
 
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
         CNX = MPNORM.VPOCHA(NLCF,1)
         CNY = MPNORM.VPOCHA(NLCF,2)
         ORIENT = CNX * DXG + CNY * DYG
         ORIENT = SIGN(1.0D0,ORIENT)
         IF(ORIENT .NE. 1.0D0)THEN
            LOGAN = .TRUE.
            MESERR(1:30)=
     &           'PRET , subroutine pre311.eso. '
            GOTO 9999
         ENDIF
         CNX = CNX * ORIENT
         CNY = CNY * ORIENT
C
C********** Cosinus directeurs de tangent 2D
C
         CTX = -1.0D0 * CNY
         CTY = CNX
C
C
C******* Les autres MELVALs
C
C
C******* N.B.: On suppose qu'on a déjà controlle RO, P, > 0
C              Y \in (1,3)
C
C******* NGCEG = NGCED -> Mur
C
         IF( NGCEG .EQ. NGCED)THEN
            ROG = MPROC.VPOCHA(NLCEG , 1)
            PG = MPPC.VPOCHA(NLCEG, 1)
            UXG = MPVITC.VPOCHA(NLCEG , 1)
            UYG = MPVITC.VPOCHA(NLCEG , 2)
            UNG = UXG * CNX + UYG * CNY
            UTG = UXG * CTX + UYG * CTY
C
C********** Son etat droite
C
            ROD = ROG
            PD = PG
            UND = -1.0D0 * UNG
            UTD = UTG
C
C********** Les fractiones massiques
C
            DO I1 = 1, NESP, 1
               FRAMAS.FRAMG(I1) =  MPYC.VPOCHA(NLCEG,I1)
               FRAMAS.FRAMD(I1) =  FRAMAS.FRAMG(I1)
            ENDDO
C
C********** Les scalaires passifs
C
            DO I1 = 1, NSCA, 1
               SCALPA.FRAMG(I1) =  MPSC.VPOCHA(NLCEG,I1)
               SCALPA.FRAMD(I1) =  SCALPA.FRAMG(I1)
            ENDDO
C
C************* Fin cas mur
C
         ELSE
C
C************* Etat gauche
C
            ROG = MPROC.VPOCHA(NLCEG, 1)
            PG = MPPC.VPOCHA(NLCEG, 1)
            UXG = MPVITC.VPOCHA(NLCEG , 1)
            UYG = MPVITC.VPOCHA(NLCEG , 2)
            UNG = UXG * CNX + UYG * CNY
            UTG = UXG * CTX + UYG * CTY
C
C********** Etat droit
C
            ROD = MPROC.VPOCHA(NLCED,1)
            PD = MPPC.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)
            UND = UXD * CNX + UYD * CNY
            UTD = UXD * CTX + UYD * CTY
C
C********** Les fractions massiques
C
            DO I1 = 1, NESP, 1
               FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1)
               FRAMAS.FRAMD(I1) = MPYC.VPOCHA(NLCED,I1)
            ENDDO
C
C********** Les scalaires passifs
C
            DO I1 = 1, NSCA, 1
               SCALPA.FRAMG(I1) = MPSC.VPOCHA(NLCEG,I1)
               SCALPA.FRAMD(I1) = MPSC.VPOCHA(NLCED,I1)
            ENDDO
         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
         MELVUN.VELCHE(1,NLCF) = UNG
         MELVUN.VELCHE(3,NLCF) = UND
         MELVUT.VELCHE(1,NLCF) = UTG
         MELVUT.VELCHE(3,NLCF) = UTD
         MELVNX.VELCHE(1,NLCF) = CNX
         MELVNY.VELCHE(1,NLCF) = CNY
         MELT1X.VELCHE(1,NLCF) = CTX
         MELT1Y.VELCHE(1,NLCF) = CTY
         DO I1 = 1, NESP, 1
            MELVA1 = MCHAMY.IELVAL(I1)
            MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1)
            MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1)
         ENDDO
         DO I1 = 1, NSCA, 1
            MELVA1 = MCHAMS.IELVAL(I1)
            MELVA1.VELCHE(1,NLCF) = SCALPA.FRAMG(I1)
            MELVA1.VELCHE(3,NLCF) = SCALPA.FRAMD(I1)
         ENDDO
      ENDDO
C
C**** Desactivation des SEGMENTs
C
      SEGDES IPT1
      SEGDES IPT2
C
      SEGDES MPROC
      SEGDES MPVITC
      SEGDES MPPC
      SEGDES MPNORM
C
      SEGDES MELRO
      SEGDES MELP
      SEGDES MELVUN
      SEGDES MELVUT
      SEGDES MELVNX
      SEGDES MELVNY
      SEGDES MELT1X
      SEGDES MELT1Y
C
      IF(NESP .GT. 0)THEN
         SEGDES MPYC
         DO I1 = 1, NESP
            MELVA1 = MCHAMY.IELVAL(I1)
            SEGDES MELVA1
         ENDDO
         SEGDES MCHAMY
         SEGSUP FRAMAS
      ENDIF
      IF(NSCA .GT. 0)THEN
         SEGDES MPSC
         DO I1 = 1, NSCA
            MELVA1 = MCHAMS.IELVAL(I1)
            SEGDES MELVA1
         ENDDO
         SEGDES MCHAMS
         SEGSUP SCALPA
      ENDIF
C
C**** Destruction du MELNTI correspondance local/global
C
      SEGSUP MLENT1
C
 9999 CONTINUE
C
      RETURN
      END