Sources Esope | Cast3M

pente3
C PENTE3    SOURCE    OF166741  24/12/13    21:16:55     12097          
      SUBROUTINE PENTE3(NFAC,MELEFL,MLECEN,NCOMP,MPOCHP,
     &     MPOGRA,MPOMIN,MPOMAX,MPOALP)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PENTE3
C
C DESCRIPTION       :  Cette subroutine est appellée par la subroutine
C                      PENTE1 (calcul du gradient d'un CHPOINT 2D de type
C                      CENTRE)
C                      Elle contient la partie du calcul du limiteur.
C
C LANGUAGE          :  ESOPE 2000 avec extensions CISI
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
C AUTEUR (Modif.)   :  R. MOREL, DRN/DMT/SEMT/TTMF
C
C************************************************************************
C
C
C APPELES (E/S)     : none
C
C APPELES (Calcul)  : none
C
C
C************************************************************************
C
C ENTREES  : NFAC : nombre de faces
C
C            MELEFL : pointeur du  MELEME  'FACEL'
C
C            MLECEN : pointeur de MLENTI qui contient la table
C                     numerotation global/local de CENTREs
C
C            NCOMP  : nombre de composantes de CHPOINT dont on veut
C                     calculer les pentes
C
C            MPOCHP : pointeur de MPOVAL de CHPOINT dont on  veut le
C                     gradient
C
C            MPOGRA : pointeur de MPOVAL  du  gradient
C
C            MPOMIN : pointeur de MPOVAL du  minimum  sur le stencil
C
C            MPOMAX : pointeur de MPOVAL du  maximum  sur le stencil
C
C SORTIES  : MPOALP : pointeur de MELVAL du limiteur
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE : Cree le 6-4-1998
C
C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
-INC PPARAM
-INC CCOPTIO
C
      INTEGER NFAC,NCOMP
     &     ,NLCF,NGCF,NGCEG,NLCEG,NGCED,NLCED
     &     ,INDCEL,I1,IGR,IDIMP1
C
C**** Pour l'instant 2D
C
      REAL*8 XG,YG,XD,YD,XC,YC,DXG,DYG,DXD,DYD,DG,DD,DT
     &     ,PHI,PHIMAX,PHIMIN,DPHIST,GRADX,GRADY,DPHI0,DPHI1,ALPHA
     &     ,VALEUR
C
C**** Extension 3D
C
      REAL*8 ZG,ZD,ZC,DZG,DZD,GRADZ
C
C
-INC SMCOORD
-INC SMCHPOI
-INC SMELEME
-INC SMLENTI
      POINTEUR MPOMIN.MPOVAL, MPOMAX.MPOVAL, MPOALP.MPOVAL,
     &         MPOCHP.MPOVAL, MPOGRA.MPOVAL
      POINTEUR MELEFL.MELEME, MLECEN.MLENTI
C
C**** N.B. Tous les pointeurs ici sont déjà activés!
C
C**** Maillage FACEL
C     Boucle sur les faces
C
C
C**** NGCEG = Numero global centre d'elt "gauche"
C     NLCEG = Numero local centre d'elt "gauche"
C     NGCF  =   "    global centre de face
C     NLCF  = numero local  centre de face
C     NGCEG = Numero global centre d'elt "gauche"
C     NLCEG = Numero local centre d'elt "gauche"
C
      IDIMP1 = IDIM+1
C
C Cas de la dimension 2
C
      IF (IDIM .EQ. 2) THEN
         DO NLCF = 1, NFAC
            NGCEG = MELEFL.NUM(1,NLCF)
            NLCEG = MLECEN.LECT(NGCEG)
            NGCF  = MELEFL.NUM(2,NLCF)
            NGCED = MELEFL.NUM(3,NLCF)
            NLCED = MLECEN.LECT(NGCED)
            IF(NGCEG .EQ. NGCED)THEN
C
C********** Limiteur dans le cas mur
C
C
C********** Coordonees et parametres geometriques
C
C
C                 XCOOR   : table de coordonnes de points;
C                           pour le point de numero global NG
C                            X_NG = XCOOR((NG-1)*(IDIM+1)+1)
C                            Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
C                            Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
C                            (dans notre cas IDIM=2, i.e. 2D)
C
               INDCEL = (NGCEG-1)*IDIMP1
               XG = XCOOR(INDCEL+1)
               YG = XCOOR(INDCEL+2)
               INDCEL = (NGCF-1)*IDIMP1
               XC = XCOOR(INDCEL+1)
               YC = XCOOR(INDCEL+2)
               DXG = XC - XG
               DYG = YC - YG
C
C********** Boucle sur le composantes
C
               DO I1 = 1, NCOMP
                  IGR = 2*I1 - 1
C
C************* Limiteur
C
                  PHI = MPOCHP.VPOCHA(NLCEG,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
                  DPHI0 = GRADX * DXG + GRADY * DYG
                  ALPHA = MPOALP.VPOCHA(NLCEG,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*0.5D0
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*0.5D0
                  ENDIF
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCEG,I1) = ALPHA
               ENDDO
            ELSE
C
C******* NGCEG != NGCED
C
               INDCEL = (NGCEG-1)*IDIMP1
               XG = XCOOR(INDCEL+1)
               YG = XCOOR(INDCEL+2)
               INDCEL = (NGCED-1)*IDIMP1
               XD = XCOOR(INDCEL+1)
               YD = XCOOR(INDCEL+2)
               INDCEL = (NGCF-1)*IDIMP1
               XC = XCOOR(INDCEL+1)
               YC = XCOOR(INDCEL+2)
C
               DXG = XC - XG
               DYG = YC - YG
               DXD = XC - XD
               DYD = YC - YD
               DG = DXG * DXG + DYG * DYG
               DG = SQRT(DG)
               DD = DXD * DXD + DYD * DYD
               DD = SQRT(DD)
               DT = DG + DD
C
C********** Boucle sur le composantes
C
               DO I1 = 1, NCOMP
                  IGR = 2*I1 - 1
C
C************* Limiteur a gauche
C
                  PHI = MPOCHP.VPOCHA(NLCEG,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
                  DPHI0 = GRADX * DXG + GRADY * DYG
                  ALPHA = MPOALP.VPOCHA(NLCEG,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*DG/DT
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*DG/DT
                  ENDIF
 
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCEG,I1) = ALPHA
C
C************* Limiteur a droite
C
                  PHI = MPOCHP.VPOCHA(NLCED,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCED,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
                  DPHI0 = GRADX * DXD + GRADY * DYD
                  ALPHA = MPOALP.VPOCHA(NLCED,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*DD/DT
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*DD/DT
                  ENDIF
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCED,I1) = ALPHA
               ENDDO
            ENDIF
         ENDDO
C
C Cas de la dimension 3
C
      ELSE
         DO NLCF = 1, NFAC
            NGCEG = MELEFL.NUM(1,NLCF)
            NLCEG = MLECEN.LECT(NGCEG)
            NGCF  = MELEFL.NUM(2,NLCF)
            NGCED = MELEFL.NUM(3,NLCF)
            NLCED = MLECEN.LECT(NGCED)
            IF(NGCEG .EQ. NGCED)THEN
C
C********** Limiteur dans le cas mur
C
C
C********** Coordonees et parametres geometriques
C
C
C                 XCOOR   : table de coordonnes de points;
C                           pour le point de numero global NG
C                            X_NG = XCOOR((NG-1)*(IDIM+1)+1)
C                            Y_NG = XCOOR((NG-1)*(IDIM+1)+2)
C                            Z_NG = XCOOR((NG-1)*(IDIM+1)+3)
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********** Boucle sur le composantes
C
               DO I1 = 1, NCOMP
                  IGR = IDIM*(I1-1)+1
C
C************* Limiteur
C
                  PHI = MPOCHP.VPOCHA(NLCEG,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
                  GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
                  DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
                  ALPHA = MPOALP.VPOCHA(NLCEG,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*0.5D0
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*0.5D0
                  ENDIF
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCEG,I1) = ALPHA
               ENDDO
            ELSE
C
C******* NGCEG != NGCED
C
               INDCEL = (NGCEG-1)*IDIMP1
               XG = XCOOR(INDCEL+1)
               YG = XCOOR(INDCEL+2)
               ZG = XCOOR(INDCEL+3)
               INDCEL = (NGCED-1)*IDIMP1
               XD = XCOOR(INDCEL+1)
               YD = XCOOR(INDCEL+2)
               ZD = XCOOR(INDCEL+3)
               INDCEL = (NGCF-1)*IDIMP1
               XC = XCOOR(INDCEL+1)
               YC = XCOOR(INDCEL+2)
               ZC = XCOOR(INDCEL+3)
C
               DXG = XC - XG
               DYG = YC - YG
               DZG = ZC - ZG
               DXD = XC - XD
               DYD = YC - YD
               DZD = ZC - ZD
               DG = DXG * DXG + DYG * DYG + DZG * DZG
               DG = SQRT(DG)
               DD = DXD * DXD + DYD * DYD + DZD * DZD
               DD = SQRT(DD)
               DT = DG + DD
C
C********** Boucle sur le composantes
C
               DO I1 = 1, NCOMP
                  IGR = IDIM*(I1-1)+1
C
C************* Limiteur a gauche
C
                  PHI = MPOCHP.VPOCHA(NLCEG,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCEG,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCEG,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCEG,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCEG,IGR+1)
                  GRADZ = MPOGRA.VPOCHA(NLCEG,IGR+2)
                  DPHI0 = GRADX * DXG + GRADY * DYG + GRADZ * DZG
                  ALPHA = MPOALP.VPOCHA(NLCEG,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*DG/DT
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*DG/DT
                  ENDIF
 
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCEG,I1) = ALPHA
C
C************* Limiteur a droite
C
                  PHI = MPOCHP.VPOCHA(NLCED,I1)
                  PHIMAX = MPOMAX.VPOCHA(NLCED,I1)
                  PHIMIN = MPOMIN.VPOCHA(NLCED,I1)
                  DPHIST = PHIMAX - PHIMIN
                  GRADX = MPOGRA.VPOCHA(NLCED,IGR)
                  GRADY = MPOGRA.VPOCHA(NLCED,IGR+1)
                  GRADZ = MPOGRA.VPOCHA(NLCED,IGR+2)
                  DPHI0 = GRADX * DXD + GRADY * DYD + GRADZ * DZD
                  ALPHA = MPOALP.VPOCHA(NLCED,I1)
                  IF(ABS(DPHI0) .LE. (DPHIST*1.0D-8))THEN
                     VALEUR = 1.0D0
                  ELSEIF(DPHI0 .GT. 0.0D0)THEN
                     DPHI1 = PHIMAX - PHI
                     VALEUR = DPHI1/DPHI0*DD/DT
                  ELSE
                     DPHI1 = PHIMIN - PHI
                     VALEUR = DPHI1/DPHI0*DD/DT
                  ENDIF
                  ALPHA = MIN(ALPHA, VALEUR)
                  MPOALP.VPOCHA(NLCED,I1) = ALPHA
               ENDDO
            ENDIF
         ENDDO
      ENDIF
C
      RETURN
      END