C HDEBI2 SOURCE PV090527 26/04/30 21:15:38 12529 C HDEBI1 SOURCE CHAT 97/12/23 22:20:50 3021 SUBROUTINE HDEBI2(IPMAHY,IPFACE,IPDARC,ICHP0,IPCHEL, S ICHP1,ICHP2,ITTH,INORM,IORIE,ISURF,IPFORC,IFORC,IRET) C----------------------------------------------------------------------- C Calcul le débit aux faces et orientation suivant le sens de la normale C lorsqu'on connait la concentrations au centre C----------------------------------------------------------------------- C C--------------------------- C Parametres Entree/Sortie : C--------------------------- C C E/ IPMAHY : MELEME des connectivités éléments/faces pour Darcy C E/ IPFACE : MELEME des points FACEs -1 C E/ IPDARC : RIGIDITE de sous type DARCY (contient RE ). C E/ ICHP0 : CHPO face des traces de concentration au temps n C E/ IPCHEL : MCHAML des orientations de normale (1=out,-1=in) C E/ ICHP1 : CHPO concentration au centre au temps n C C Parametre optionnel C E/ ICHP2 : CHPO face flux de vitesse C C /S IRET : CHPO face des flux les noms des composantes sont C ceux des composantes de ITPN et IPCH1. C Si ICHP2 est donné on ajoute le flux convectif C C---------------------- C Variables en COMMON : C---------------------- C C E/ IFOUR : cf CCOPTIO C E/ IDIM : cf CCOPTIO C C---------------------- C Tableaux de travail : C---------------------- C C ICPR(I)=J : Le noeud I a le numero J dans le MELEME des faces C Correspondance numérotation globale/locale C ITES : Nombre de noeuds FACE C NNGOT : Nombre de noeuds total du domaine C IVAA(I) : indice du CHAMPOINT au Ieme noeud global C C----------------------------------------------------------------------- C C Langage : ESOPE + FORTRAN77 C C Auteurs : 09/93 F.DABBENE - Cas permanent C 09/94 X.NOUVELLON - Extension au cas transitoire C 02/96 L.V.BENET - Prise en compte de forces de volume C 05/98 F.AURIOL expression en fonction des concentrations C (charges) et traces de concentrations (traces C de concentrations) possibilités de champs C à plusieurs composantes, en transitoire. C C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMCHPOI POINTEUR MCHPO5.MCHPOI,MCHPO6.MCHPOI,MSOUP6.MSOUPO -INC SMCHAML -INC SMRIGID -INC SMCOORD * CHARACTER*(LOCOMP) MOREFD,MOREFP SEGMENT IPMAHY INTEGER MAHYBR(NSOUS) ENDSEGMENT SEGMENT ICCPR INTEGER ICPR(NNGOT) ENDSEGMENT SEGMENT IORGA INTEGER IVAA(ITES), IVBB(ITES) ENDSEGMENT SEGMENT ITRAV REAL*8 RLIGNE(NBDDL) REAL*8 FLFOR(NBDDL),RFOR(NBDDL) ENDSEGMENT C C *' Initialisations' C MRIGID = IPDARC RI1 = IPFORC MCHELM = IPCHEL IPT1 = IPFACE MCHPO1 = ICHP0 MCHPO2 = ICHP1 MCHPO3 = ICHP2 IRET = 0 MCHEL2 = IORIE MCHPO4 = IFORC MCHPO5 = INORM MCHPO6 = ISURF C SEGACT IPMAHY NBMAIL = MAHYBR(/1) SEGACT MRIGID * *' Creation du tableau ICPR pour le maillage IPT1' * NNGOT = nbpts SEGINI ICCPR SEGACT IPT1 N2 = IPT1.NUM(/2) IK = 0 DO 15 I2=1,N2 K = IPT1.NUM(1,I2) IF (ICPR(K).EQ.0) THEN IK = IK + 1 ICPR(K) = IK ENDIF 15 CONTINUE SEGDES IPT1 ITES = IK * *' Activation du MPOVAL du CHPO de traces de concentrations * SEGACT MCHPO1 MSOUP1 = MCHPO1.IPCHP(1) SEGDES MCHPO1 SEGACT MSOUP1 MPOVA1 = MSOUP1.IPOVAL SEGACT MPOVA1 NBCOMP= MPOVA1.VPOCHA(/2) * * Activation du MPOVAL du CHPO des concentrations au centre * SEGACT MCHPO2 MSOUP2 = MCHPO2.IPCHP(1) SEGDES MCHPO2 SEGACT MSOUP2 MPOVA2 = MSOUP2.IPOVAL SEGACT MPOVA2 * * Activation du MPOVAL du CHPO flux de vitesses aux faces * IF (ICHP2.NE.0) THEN SEGACT MCHPO3 MSOUP3 = MCHPO3.IPCHP(1) SEGDES MCHPO3 SEGACT MSOUP3 MPOVA3 = MSOUP3.IPOVAL SEGACT MPOVA3 ENDIF * * activation des objets liés à la présence d'une force volumique * IF (IFORC.NE.0) THEN * * Activation du MPOVAL du CHPO force appuyé au centre des éléments volumiques * SEGACT MCHPO4 MSOUP4 = MCHPO4.IPCHP(1) SEGDES MCHPO4 SEGACT MSOUP4 MPOVA4 = MSOUP4.IPOVAL SEGDES MSOUP4 SEGACT MPOVA4 * * Activation du MPOVAL du CHPO des vecteurs normales appuyé au centre des faces * SEGACT MCHPO5 MSOUP5 = MCHPO5.IPCHP(1) SEGDES MCHPO5 SEGACT MSOUP5 MPOVA5 = MSOUP5.IPOVAL SEGDES MSOUP5 SEGACT MPOVA5 * * Activation du MPOVAL du CHPO des surfaces appuyé au centre des faces * SEGACT MCHPO6 MSOUP6 = MCHPO6.IPCHP(1) SEGDES MCHPO6 SEGACT MSOUP6 MPOVA6 = MSOUP6.IPOVAL SEGDES MSOUP6 SEGACT MPOVA6 * * Activation du MCHEL des orientations des faces * SEGACT MCHEL2 * * Activation du MRIGI de la matrice masse hybride * SEGACT RI1 ENDIF * * On recherche l ordre des traces de concentrations par rapport à IPT1 * SEGINI IORGA MELEME = MSOUP1.IGEOC SEGACT MELEME N2 = NUM(/2) DO 25 I2=1,N2 K = NUM(1,I2) IF (ICPR(K).EQ.0) THEN INTERR(1) = K MOTERR(1:8) = 'FACE ' CALL ERREUR(64) SEGDES MELEME, MSOUP1 SEGDES MCHPO2, MRIGID, IPMAHY SEGSUP ICCPR, IORGA RETURN ELSE IVAA(ICPR(K)) = I2 ENDIF 25 CONTINUE SEGDES MELEME CALL INITI( IVBB,ITES,0) * * Construction de CHPOIN resultat les composantes ont les noms * de celles des concentrations au centre ( ou aux faces) * SEGACT IPT1 NPN=IPT1.NUM(/2) SEGDES IPT1 NSOUPO=1 NAT=1 SEGINI MCHPOI MTYPOI=' ' MOCHDE=' CHPOIN CREE PAR HDEBI1 ' IFOPOI=IFOUR JATTRI(1)=2 NC=NBCOMP SEGINI MSOUPO IPCHP(1)=MSOUPO DO 5 L=1,NBCOMP NOCOMP(L)=MSOUP1.NOCOMP(L) NOHARM(L)=MSOUP1.NOHARM(L) 5 CONTINUE IGEOC=IPFACE N=NPN SEGINI MPOVAL IPOVAL=MPOVAL NB=N*NC CALL INITD(VPOCHA,NB,0.D0) IF(ITTH.EQ.1) THEN C C cas des traces de charges récupération du nom des composantes C NBMAIL = MAHYBR(/1) DO 27 IMAIL=1,NBMAIL IF (MAHYBR(IMAIL).NE.0) THEN DESCR = IRIGEL(3,IMAIL) SEGACT DESCR MOREFD = LISDUA(1) MOREFP = LISINC(1) SEGDES DESCR GOTO 30 ENDIF 27 CONTINUE 30 CONTINUE NOCOMP(1)=MOREFD ENDIF * * * C C-------------------------------------------------- *' Boucle 310 sur les OBJETS RIGIDITES ELEMENTAIRES' C-------------------------------------------------- C ITELEM = 0 SEGACT MCHELM DO 310 IRI=1,NBMAIL C C Récupération MELEME ou Darcy est défini C MELEME = MAHYBR(IRI) IF (MELEME.EQ.0) GOTO 310 SEGACT MELEME N1 = NUM(/1) N2 = NUM(/2) C C Récupération des infos pour la zone IRI dans le chapeau MRIGID C DESCR = IRIGEL(3,IRI) SEGACT DESCR NBDDL = LISDUA(/2) SEGDES DESCR SEGINI ITRAV xMATRI = IRIGEL(4,IRI) SEGACT xMATRI C C Activation du MELVAL du MCHAML d'orientation C MCHAML = ICHAML(IRI) SEGACT MCHAML MELVAL = IELVAL(1) SEGDES MCHAML SEGACT MELVAL * * Activation des objets necessaires à la prise en compte des forces de volumes * IF (IFORC.NE.0) THEN MCHAM2 = MCHEL2.ICHAML(IRI) SEGACT MCHAM2 MELVA2 = MCHAM2.IELVAL(1) SEGDES MCHAM2 SEGACT MELVA2 xMATR1 = RI1.IRIGEL(4,IRI) SEGACT xMATR1 ELSE DO 35 I=1,NBDDL RFOR(I)=0.D0 35 CONTINUE ENDIF C C------------------------------------------ *' Boucle 300 sur les MATRICES ELEMENTAIRES.' C------------------------------------------ C DO 300 I2=1,N2 ITELEM = ITELEM + 1 IF (IFORC.NE.0) THEN * *- calcul des flux de forces aux faces de l'element * DO 55 IDDL=1,NBDDL FLFOR(IDDL)= 0.D0 IPOPTS = ICPR(NUM(IDDL,I2)) DO 50 I=1,IDIM FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I) * S MELVA2.VELCHE(IDDL,I2) * MPOVA4.VPOCHA(ITELEM,I) * S MPOVA6.VPOCHA(IPOPTS,1) 50 CONTINUE 55 CONTINUE * *- Construction du tableau aux faces M.FORCE * * XMATR1 = IMATR1.IMATTT(I2) * SEGACT XMATR1 DO 65 I=1,NBDDL RFOR(I)=0.D0 DO 60 J=1,NBDDL RFOR(I) = RFOR(I) + XMATR1.RE(I,J,i2)*FLFOR(J) 60 CONTINUE 65 CONTINUE * SEGDES XMATR1 ENDIF * * Recuperation de la matrice elementaire * * XMATRI = IMATTT(I2) * SEGACT XMATRI * *- De la somme des coefs pour une ligne *- -1 t *- LIGNE = RE * DIV *- -1 t DO 140 I=1,NBDDL RLIGNE(I) = 0.D0 DO 130 J=1,NBDDL RLIGNE(I) = RLIGNE(I) + RE(I,J,i2) 130 CONTINUE 140 CONTINUE C C Calcul du flux aux faces C DO 200 IN=1,NBDDL NUMFA = ICPR(NUM(IN,I2)) IF (IVBB(NUMFA).EQ.0) THEN VVV= 0.D0 IF(ICHP2.NE.0)THEN VVV=MPOVA3.VPOCHA(NUMFA,1) ENDIF DO 180 K=1,NBCOMP VA1 = 0.D0 VA2 = 0.D0 DO 190 JN=1,NBDDL VA1 = VA1+RE(IN,JN,i2)*(MPOVA1.VPOCHA( S IVAA(ICPR(NUM(JN,I2))),K)-RFOR(JN)) 190 CONTINUE VA2=RLIGNE(IN)*MPOVA2.VPOCHA(ITELEM,K) VA3= VVV*MPOVA1.VPOCHA(NUMFA,K) VPOCHA(NUMFA,K) = (VA2 -VA1+ VA3 ) * VELCHE(IN,I2) 180 CONTINUE IVBB(NUMFA)=1 ENDIF 200 CONTINUE * SEGDES XMATRI 300 CONTINUE SEGDES MELVAL, xMATRI, MELEME SEGSUP ITRAV IF (IFORC.NE.0) THEN SEGDES MELVA2, xMATR1 ENDIF 310 CONTINUE C C Nettoyage final C SEGDES MCHELM, MRIGID, IPMAHY, MSOUPO, MPOVAL, MCHEL2 SEGDES MSOUP1, MPOVA1 SEGDES MSOUP2, MPOVA2 IF (MCHPO3.NE.0) SEGDES MPOVA3,MSOUP3 IF (IFORC.NE.0) THEN SEGDES RI1 SEGDES MSOUP4, MPOVA4 SEGDES MSOUP5, MPOVA5 SEGDES MSOUP6, MPOVA6 ENDIF C SEGDES MCHPOI IRET = MCHPOI C SEGSUP IORGA, ICCPR C RETURN END