Sources Esope | Cast3M

divufn
C DIVUFN    SOURCE    FANDEUR   22/01/03    21:15:12     11136          
      SUBROUTINE DIVUFN(ICHP2,ICLIM,IPFACE,IFACEL,IRE1,IRE2,IPFONC)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
 
C-----------------------------------------------------------------------
C Calcul d'un flux decentre.
C Le CHAMPOINT résultat est de support FACE.
C-----------------------------------------------------------------------
C
C---------------------------
C Parametres Entree/Sortie :
C---------------------------
 
C E/  ICHP2    : CHPOIN DES VALEURS F(THETA)
C E/  ICLIM    : CHPOIN DES CONDITIONS AUX LIMITES IMPOSEES
C E/  IPFACE   : MELEME DES POINTS FACE
C E/  IFACEL   : MELEME DES POINTS FACE POUR LES C.L.
C E/  IRE1     : Champoint de type FLUX
C E/  IRE2     : Mchaml des orientation de normale (1=out,-1=in)
C S/  IPFONC   : CHAMPOIN RESULTAT DES F(\THETA) DECENTRE
C
C----------------------
C Tableaux de travail :
C----------------------
C
C
C----------------------
C Variables en COMMON :
C----------------------
C
C  IFOUR : cf CCOPTIO.INC
C
C-----------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN77
C
C Auteurs : C. LE POTIER ET F. AURIOL 20/00
C
C-----------------------------------------------------------------------
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCHPOI
-INC SMCHAML
-INC SMCOORD
*
      SEGMENT ICCPR
         INTEGER ICPR(NNGOT)
      ENDSEGMENT
 
      SEGMENT ICCPR1
         INTEGER ICPR1(NNGOT)
      ENDSEGMENT
 
 
C
C= INITIALISATIONS
C
      MCHPO1 = IRE1
      MCHELM = IRE2
      IPT1   = IPFACE
 
      NNGOT  = nbpts
      SEGINI ICCPR1
 
*
*= Creation des tableaux ICPR et INUI pour le maillage IPT1 des FACES
*
*      WRITE(6,*) 'AVANT SEGACT'
      SEGACT IPT1
*      WRITE(6,*) 'ON A PASSE LE PREMER SEGACT'
      N2 = IPT1.NUM(/2)
      IK = 0
      DO 109 I2=1,N2
         K = IPT1.NUM(1,I2)
         IF (ICPR1(K).EQ.0) THEN
            IK      = IK + 1
            ICPR1(K) = IK
         ENDIF
 109   CONTINUE
      SEGDES IPT1
 
*      WRITE(6,*) 'APRES SEDDES IPT1'
C
C- Récupération du pointeur MPOVAL des flux
C
      SEGACT MCHPO1
      MSOUP1 = MCHPO1.IPCHP(1)
      SEGACT MSOUP1
      MPOVA1 = MSOUP1.IPOVAL
      SEGDES MSOUP1
      SEGDES MCHPO1
 
*      WRITE(6,*) 'APRES RECUPERATION DU FLUX'
C
C     DEFINITION DU CHPOIN RESULTAT
C
      NAT=1
      NSOUPO=1
      SEGINI MCHPOI
*      WRITE(6,*) 'MCHPOI'
      IPFONC=MCHPOI
      IFOPOI    = IFOUR
      JATTRI(1)=2
      NC=1
      SEGINI MSOUPO
*      WRITE(6,*) 'MSOUPO'
      NOCOMP(1)='SCAL'
      IPCHP(1)=MSOUPO
      IGEOC=IPFACE
      IPT1=IPFACE
      SEGACT IPT1
      N=IPT1.NUM(/2)
      SEGINI MPOVAL
      IPOVAL=MPOVAL
      NOHARM(1)=NIFOUR
*      SEGDES MCHPOI
*      SEGDES MSOUPO
*      SEGDES IPT1
*      WRITE(6,*) 'DEFINITION DU CHAMPOIN'
C
      IPT3=IFACEL
      SEGACT IPT3
      NBFACE=IPT3.NUM(/2)
      MCHPO2=ICHP2
      SEGACT MCHPO2
      MSOUP2=MCHPO2.IPCHP(1)
      SEGACT MSOUP2
      MPOVA2=MSOUP2.IPOVAL
      SEGACT MPOVA2
      NPCENT=MPOVA2.VPOCHA(/1)
      IPT2=MSOUP2.IGEOC
      SEGACT IPT2
C  On sait que le support de MCHPO2 est le maillage IPCENT (déja vérifié)
      NNGOT=nbpts
      SEGINI ICCPR
      DO 10 I=1,NPCENT
         K=IPT2.NUM(1,I)
         ICPR(K)=I
   10 CONTINUE
 
 
*      WRITE(6,*) 'BOUCLE SUR LES ELEMENTS'
C
C------------------------------------------------
C= Boucle  sur les ZONES ELEMENTAIRES du MCHAML
C------------------------------------------------
C
      ITELEM = 0
      SEGACT MCHELM
      SEGACT MPOVA1
      NRIGEL = IMACHE(/1)
      DO 409 IRI=1,NRIGEL
C
C Recuperation du MELEME et activation
C
         MELEME = IMACHE(IRI)
         SEGACT MELEME
         N1 = NUM(/1)
         N2 = NUM(/2)
C
C Récupération du pointeur MELVAL du MCHAML d'orientation
C
         MCHAML = ICHAML(IRI)
         SEGACT MCHAML
         MELVAL = IELVAL(1)
         SEGDES MCHAML
         SEGACT MELVAL
C
C------------------------------
C= Boucle 30 sur les ELEMENTs.
C------------------------------
C
C CALCUL DE f(THETA) DECENTRE
         DO 309 I2=1,N2
            ITELEM = ITELEM + 1
            DO 209 IN=1,N1
             VALIN1 = MPOVA1.VPOCHA(ICPR1(NUM(IN,I2)),1)*VELCHE(IN,I2)
             IFACE = ICPR1(NUM(IN,I2))
             IP  = ICPR(IPT3.NUM(1,IFACE))
             ID  = ICPR(IPT3.NUM(3,IFACE))
*             write(6,*) 'I2=',I2,'IP=',IP, 'ID=', ID
             IF (ID.EQ.I2) THEN
             ID = IP
             IP = I2
             ENDIF
             IF (VALIN1.LT.0) THEN
             VPOCHA(IFACE,1) = MPOVA2.VPOCHA(IP,1)
             ELSE
             VPOCHA(IFACE,1) = MPOVA2.VPOCHA(ID,1)
             ENDIF
 209         CONTINUE
 309      CONTINUE
         SEGDES MELVAL, MELEME
 409   CONTINUE
 
 
      IF(ICLIM.NE.0)THEN
      MCHPO4=ICLIM
      SEGACT MCHPO4
      NSOUP4=MCHPO4.IPCHP(/1)
      CALL INITI(ICPR,NNGOT,0)
      IPT4=IPFACE
      DO 30 I=1,NBFACE
          K=IPT4.NUM(1,I)
          ICPR(K)=I
   30 CONTINUE
      DO 40 I=1,NSOUP4
          MSOUP4=MCHPO4.IPCHP(I)
          SEGACT MSOUP4
          IPT5=MSOUP4.IGEOC
          MPOVA5=MSOUP4.IPOVAL
          SEGACT IPT5,MPOVA5
          NBP5=IPT5.NUM(/2)
          DO 50 J=1,NBP5
              NUMP=IPT5.NUM(1,J)
              VPOCHA(ICPR(NUMP),1)=MPOVA5.VPOCHA(J,1)
   50     CONTINUE
      SEGDES IPT5,MPOVA5,MSOUP4
   40 CONTINUE
      SEGDES MCHPO4,IPT2
 
      ENDIF
 
 
      SEGDES MCHELM
      SEGDES MPOVA1
      SEGDES MPOVAL
 
      SEGSUP ICCPR
      SEGSUP ICCPR1
C
      RETURN
      END