C HYBP1 SOURCE CHAT 09/10/09 21:18:48 6519 SUBROUTINE HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE, S ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,IRET) C----------------------------------------------------------------------- C Calcul de la charge au centre de chaque maille C----------------------------------------------------------------------- C C--------------------------- C Parametres Entree/Sortie : C--------------------------- C C E/ IPMAHY : Segment contenant le pointeur vers le meleme des C connectivites elements/faces pour les zones du MMODEL C ou on a defini DARCY. C E/ ICENTR : Pointeur vers l'objet maillage CENTRE C E/ IPFACE : Pointeur vers l'objet maillage FACE -1 C E/ IPDARC : Matrice de sous type DARCY (contient RE ). C E/ ITPN : CHPO face des Traces de charge TH au temps n C E/ THETA : Parametre de discretisation temporelle (Theta-methode) C E/ DELTAT : Pas de temps C E/ IPCK : MCHAML donnant pour chaque element Ck|K| C E/ ITP : CHPO face des traces de charge TH au temps n-1 C E/ IH : CHPO centre des charges H au temps n-1 C E/ ISOUR : CHPO centre SOUR contenant un eventuel terme source C E/ IPFORC : Matrice de sous type MASSE (contient RE ). C E/ IFORC : CHPO centre force contenant la force volumique C /S IRET : CHPO centre H resultat contenant la charge au temps n C C---------------------- C Variables en COMMON : C---------------------- C C E/ IDIM : cf CCOPTIO C E/ IFOMOD : cf CCOPTIO C E/ NIFOUR : cf CCOPTIO C C---------------------- C Tableaux de travail : C---------------------- C C ICPR(I)=J : Le noeud I a le numero local J C Correspondance numerotation globale/locale C NNGOT : Nombre de noeuds total du domaine 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 C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC CCOPTIO -INC SMCHAML -INC SMCHPOI POINTEUR MCHPO5.MCHPOI,MCHPO6.MCHPOI,MSOUP6.MSOUPO -INC SMCOORD -INC SMELEME -INC SMRIGID * CHARACTER*4 MOREFD, MOREFP CHARACTER*1 BLAN1 SEGMENT ICCPR INTEGER ICPR(NNGOT) ENDSEGMENT SEGMENT IORGA REAL*8 VAA(ITES) ENDSEGMENT SEGMENT IPMAHY INTEGER MAHYBR(NSOUS) ENDSEGMENT SEGMENT MTRAV REAL*8 VECT(NBDDL),FLFOR(NBDDL),RFOR(NBDDL) ENDSEGMENT * BLAN1=' ' * * Initialisations * IRET = 0 IPT1 = IPFACE MRIGID = IPDARC RI1 = IPFORC MCHPOI = ITPN MCHELM = IPCK MCHEL1 = IORIE MCHPO1 = ISOUR MCHPO2 = ITP MCHPO3 = IH MCHPO4 = IFORC MCHPO5 = INORM MCHPO6 = ISURF VSOUR = 0.D0 * * Initialisations se servant de IPMAHY et MRIGID * SEGACT IPMAHY NBMAIL = MAHYBR(/1) SEGACT MRIGID DO 5 IMAIL=1,NBMAIL IF (MAHYBR(IMAIL).NE.0) THEN DESCR = IRIGEL(3,IMAIL) SEGACT DESCR MOREFP = LISINC(1) SEGDES DESCR GOTO 10 ENDIF 5 CONTINUE 10 CONTINUE NNGOT = XCOOR(/1)/(IDIM+1) SEGINI ICCPR * *= Creation du tableau ICPR pour le maillage IPT1 * 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 * *= Extraction des TH et réduction sur le support IPT1 du CHPO ITPN *= Les valeurs du CHPOINT en TH sont dans VAA avec l'ordre IPT1 * SEGINI IORGA SEGACT MCHPOI NSOUPO = IPCHP(/1) DO 30 ISOU=1,NSOUPO MSOUPO = IPCHP(ISOU) SEGACT MSOUPO NC = NOCOMP(/2) NUMHAR = NOHARM(1) ITTP = 0 DO 20 I=1,NC IF (NOCOMP(I).EQ.MOREFP) THEN ITTP = I ENDIF 20 CONTINUE IF (ITTP.NE.0) THEN MELEME = IGEOC SEGACT MELEME N2 = NUM(/2) MPOVAL = IPOVAL SEGACT MPOVAL 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 MPOVAL, MELEME, MSOUPO SEGDES MCHPOI, MRIGID, IPMAHY SEGSUP ICCPR, IORGA RETURN ELSE VAA(ICPR(K)) = VPOCHA(I2,ITTP) ENDIF 25 CONTINUE SEGDES MPOVAL, MELEME ENDIF SEGDES MSOUPO 30 CONTINUE SEGDES MCHPOI * *= Création des chapeaux du CHPO centre P résultat * NSOUPO = 1 NAT = 1 SEGINI MCHPOI IRET = MCHPOI DO 31 ITIT=1,72 MOCHDE(ITIT:ITIT)=BLAN1 31 CONTINUE MTYPOI = 'CENTRE ' IFOPOI = IFOMOD JATTRI(1) = 1 NC = 1 SEGINI MSOUPO IPCHP(1) = MSOUPO SEGDES MCHPOI NOHARM(1) = NUMHAR IGEOC = ICENTR NOCOMP(1) = 'H ' IPT2 = ICENTR SEGACT IPT2 N = IPT2.NUM(/2) SEGDES IPT2 SEGINI MPOVAL IPOVAL = MPOVAL SEGDES MSOUPO * * Activation du MPOVAL du CHPO centre SOUR * IF (ISOUR.NE.0) THEN SEGACT MCHPO1 MSOUP1 = MCHPO1.IPCHP(1) SEGDES MCHPO1 SEGACT MSOUP1 MPOVA1 = MSOUP1.IPOVAL SEGDES MSOUP1 SEGACT MPOVA1 ENDIF * * Activation du MPOVAL du CHPO face trace de charge au temps (n-1) * IF (ITP.NE.0) THEN SEGACT MCHPO2 MSOUP2 = MCHPO2.IPCHP(1) SEGDES MCHPO2 SEGACT MSOUP2 MPOVA2 = MSOUP2.IPOVAL SEGDES MSOUP2 SEGACT MPOVA2 ENDIF * * Activation du MPOVAL du CHPO centre charge au temps (n-1) * IF (IH.NE.0) THEN SEGACT MCHPO3 MSOUP3 = MCHPO3.IPCHP(1) SEGDES MCHPO3 SEGACT MSOUP3 MPOVA3 = MSOUP3.IPOVAL SEGDES MSOUP3 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 MCHEL1 * * Activation du MRIGI de la matrice masse hybride * SEGACT RI1 ENDIF * * Activation de Ck|K| * IF (IPCK.NE.0) SEGACT MCHELM * *--------------------------------------- *= Boucle sur les maillages elementaires *--------------------------------------- * ITELEM = 0 DO 110 IMAIL=1,NBMAIL IF (MAHYBR(IMAIL).EQ.0) GOTO 110 * * Récupération des infos pour la zone IMAIL ou Darcy est définie * MELEME = IRIGEL(1,IMAIL) xMATRI = IRIGEL(4,IMAIL) SEGACT MELEME NBDDL = NUM(/1) NBELEM = NUM(/2) SEGACT xMATRI SEGINI MTRAV * * Activation des objets necessaires à la prise en compte des forces de volumes * IF (IFORC.NE.0) THEN MCHAM1 = MCHEL1.ICHAML(IMAIL) SEGACT MCHAM1 MELVA1 = MCHAM1.IELVAL(1) SEGDES MCHAM1 SEGACT MELVA1 xMATR1 = RI1.IRIGEL(4,IMAIL) SEGACT xMATR1 ELSE DO 35 I=1,NBDDL RFOR(I)=0.D0 35 CONTINUE ENDIF * C C Activation du MELVAL du MCHAML Ck|K| C IF (IPCK.NE.0) THEN MCHAML = ICHAML(IMAIL) SEGACT MCHAML MELVAL = IELVAL(1) SEGDES MCHAML SEGACT MELVAL ENDIF * *------------------------------------------------------- * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL *------------------------------------------------------- * DO 100 IEL=1,NBELEM ITELEM = ITELEM + 1 * *- Calcul de VECT = DIV * RE , CONSD = VECT * DIV * * XMATRI = IMATTT(IEL) * SEGACT XMATRI CONSD = 0.D0 DO 45 J=1,NBDDL VECT(J) = 0.D0 DO 40 I=1,NBDDL VECT(J) = VECT(J) + RE(I,J,iel) 40 CONTINUE CONSD = CONSD + VECT(J) 45 CONTINUE * SEGDES XMATRI * *- Recuperation de SOUR et de H(n-1) * IF (ISOUR.NE.0) VSOUR=MPOVA1.VPOCHA(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,IEL)) DO 50 I=1,IDIM FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I) $ *MELVA1.VELCHE(IDDL,IEL) * MPOVA4 $ .VPOCHA(ITELEM,I) *MPOVA6.VPOCHA(IPOPTS,1) 50 CONTINUE 55 CONTINUE * *- Construction du tableau aux faces M.FORCE * * XMATR1 = IMATR1.IMATTT(IEL) * SEGACT XMATR1 DO 65 I=1,NBDDL RFOR(I)=0.D0 DO 60 J=1,NBDDL RFOR(I) = RFOR(I) + XMATR1.RE(I,J,iel)*FLFOR(J) 60 CONTINUE 65 CONTINUE * SEGDES XMATR1 ENDIF * *- Calcul de p pour l'element IEL * VALP = 0.D0 IF (IPCK.EQ.0) THEN DO 80 I=1,NBDDL NUMFA = ICPR(NUM(I,IEL)) VALP = VALP + VECT(I) * (VAA(NUMFA) - RFOR(I)) 80 CONTINUE VPOCHA(ITELEM,1) = (VALP + VSOUR) / CONSD ELSE BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT) DO 90 I=1,NBDDL NUMFA = ICPR(NUM(I,IEL)) VALP = VALP + VECT(I) * ( THETA*VAA(NUMFA) S + (1.D0-THETA)*MPOVA2.VPOCHA(NUMFA,1)) 90 CONTINUE VALP = (VALP + VSOUR) * BETA / CONSD VH = MPOVA3.VPOCHA(ITELEM,1) VPOCHA(ITELEM,1) = VALP + (1.D0-BETA)*VH ENDIF 100 CONTINUE IF (IFORC.NE.0) SEGDES xMATR1, MELVA1 SEGDES xMATRI, MELEME IF (IPCK.NE.0) SEGDES MELVAL SEGSUP MTRAV 110 CONTINUE SEGDES MPOVAL SEGDES MRIGID, IPMAHY IF (IPCK.NE.0) SEGDES MCHELM, MPOVA2, MPOVA3 IF (ISOUR.NE.0) SEGDES MPOVA1 IF (IFORC.NE.0) S SEGDES MPOVA4, MPOVA5, MPOVA6, S MCHEL1, RI1 SEGSUP ICCPR , IORGA * RETURN END