C SMTP1 SOURCE CHAT 09/10/09 21:24:23 6519 SUBROUTINE SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK, S ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1) C----------------------------------------------------------------------- C 1) Calcul du second membre du systeme en trace de charge dans le cas C de la résolution des équations de Darcy par EFMH avec le modèle DARCY. C 2) Calcul du second membre du systeme en trace de charge et/ou d'une C contribution matricielle provenant de la convection dans le cas de la C résolution d'une équation de diffusion-convection par EFMH avec le C modèle DARCY. 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/ IPFACE : MELEME de type POI1 des faces C E/ IPRIGI : RIGIDITE de sous type 'DARCY' C E/ THETA : Parametre de discretisation temporelle (Theta-methode) C pour le terme de diffusion C E/ THETAC : Parametre de discretisation temporelle (Theta-methode) C pour le terme de convection C E/ DELTAT : Pas de temps C E/ IPCK : MCHAML donnant pour chaque element Ck|K| C E/ ITP : CHPO faces des traces de charge TH au temps precedent C E/ IH : CHPO centre des charges H au temps precedent C E/ ISOUR : CHPO centre des sources de composante SOUR C E/ IFLUX : CHPO face des flux de composante FLUX C E/ IORIE : MCHAML orientation des normales C /S MCHPOI : CHPOINT face de composante FLUX C /S IRI1 : RIGIDITE de sous type CONVEFMH C C---------------------- C Variables en COMMON : C---------------------- C C E/ IFOMOD : cf CCOPTIO C E/ NIFOUR : cf CCOPTIO C C--------------------- C Variables internes : C--------------------- C C LHS1 = 1 si on doit créer une matrice (Left Hand Side); 0 sinon. C RHS1 = 1 si on doit créer un second membre (Right Hand Side); 0 sinon. 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 12/95 F.DABBENE - Extension à la convection C C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC CCOPTIO -INC CCHAMP -INC SMCOORD -INC SMCHPOI -INC SMCHAML -INC SMELEME -INC SMRIGID * SEGMENT IPMAHY INTEGER MAHYBR(NSOUS) ENDSEGMENT SEGMENT ICCPR INTEGER ICPR(NNGOT) ENDSEGMENT SEGMENT MTRAV REAL*8 RLIGNE(NBDDL),RCOLON(NBDDL),FLUORI(NBDDL) ENDSEGMENT * *- Initialisations * MCHPOI = 0 IRI1 = 0 IF (IPCK.EQ.0) THEN IF (ISOUR.EQ.0) THEN RHS1 = 0 ELSE RHS1 = 1 ENDIF IF (IFLUX.EQ.0) THEN LHS1 = 0 ELSE LHS1 = 1 ENDIF ELSE RHS1 = 1 LHS1 = 0 IF (THETAC.GT.1.D-4 .AND. IFLUX.NE.0) LHS1=1 ENDIF * *- Creation du tableau ICPR pour le MELEME POI1 IPFACE * IK = 0 NNGOT = XCOOR(/1)/(IDIM+1) SEGINI ICCPR MELEME = IPFACE SEGACT MELEME N2 = NUM(/2) DO 10 I2=1,N2 K = NUM(1,I2) IF (ICPR(K).EQ.0) THEN IK = IK + 1 ICPR(K) = IK ENDIF 10 CONTINUE SEGDES MELEME * *- Creation du CHAMPOINT face - On laisse actif le MPOVAL du CHPOINT * IF (RHS1.NE.0) THEN NSOUPO = 1 NAT = 1 SEGINI MCHPOI IRET = MCHPOI MTYPOI = 'FACE ' IFOPOI = IFOMOD JATTRI(1) = 2 NC = 1 SEGINI MSOUPO IPCHP(1) = MSOUPO SEGDES MCHPOI NOHARM(1) = NIFOUR IGEOC = IPFACE NOCOMP(1) = 'FLUX' N = N2 SEGINI MPOVAL IPOVAL = MPOVAL SEGDES MSOUPO ENDIF * *- Recuperation du nombre de zone NBMAIL du MMODEL a partir IPMAHY * SEGACT IPMAHY NBMAIL = MAHYBR(/1) * *- Activation du segment RIGIDITE * MRIGID = IPRIGI SEGACT MRIGID * *- Activation du segment MPOVAL du CHAMPOINT centre de composante SOUR * IF (ISOUR.NE.0) THEN MCHPO1 = ISOUR SEGACT MCHPO1 MSOUP1 = MCHPO1.IPCHP(1) SEGDES MCHPO1 SEGACT MSOUP1 MPOVA1 = MSOUP1.IPOVAL SEGDES MSOUP1 SEGACT MPOVA1 ENDIF * *- Données transitoire * IF (IPCK.NE.0) THEN * * Activation du segment MPOVAL du CHAMPOINT face de traces de charge * MCHPO2 = ITP SEGACT MCHPO2 MSOUP2 = MCHPO2.IPCHP(1) SEGDES MCHPO2 SEGACT MSOUP2 MPOVA2 = MSOUP2.IPOVAL SEGDES MSOUP2 SEGACT MPOVA2 * * Activation du segment MPOVAL du CHAMPOINT centre de charges * MCHPO3 = IH SEGACT MCHPO3 MSOUP3 = MCHPO3.IPCHP(1) SEGDES MCHPO3 SEGACT MSOUP3 MPOVA3 = MSOUP3.IPOVAL SEGDES MSOUP3 SEGACT MPOVA3 * * Activation du MCHAML contenant Ck|K| * MCHELM = IPCK SEGACT MCHELM ENDIF * *- Activation des données pour la convection * IF (IFLUX.NE.0) THEN * * Activation du MCHAML d'orientation des normales * MCHEL1 = IORIE SEGACT MCHEL1 * * Activation du segment MPOVAL du CHAMPOINT flux aux faces * MCHPO4 = IFLUX SEGACT MCHPO4 MSOUP4 = MCHPO4.IPCHP(1) SEGDES MCHPO4 SEGACT MSOUP4 MPOVA4 = MSOUP4.IPOVAL SEGDES MSOUP4 SEGACT MPOVA4 * * Initialisation du chapeau de rigidité * IF (LHS1.NE.0) THEN NRIGE = 7 NRIGEL = NBMAIL SEGINI RI1 IRI1 = RI1 RI1.MTYMAT = 'CONVEFMH' RI1.IFORIG = IFORIG ENDIF ENDIF * *-------------------------------------------------- *= BOUCLE SUR LES MAILLAGES ELEMENTAIRES,ZONE IMAIL *-------------------------------------------------- * ITELEM = 0 DO 110 IMAIL=1,NBMAIL C C- Activation de l'objet maillage ELTFA pour la zone IMAIL C MELEME = MAHYBR(IMAIL) IF (MELEME.EQ.0) GOTO 110 SEGACT MELEME * *- Recuperation des caracteristiques de RIGIDITE de sous type DARCY *- pour la zone IMAIL en cours de traitement * xMATRI = IRIGEL(4,IMAIL) SEGACT xMATRI NBELEM = re(/3) NELRIG = NBELEM * XMATRI = IMATTT(1) * SEGACT XMATRI NBDDL = RE(/1) NLIGRP = NBDDL NLIGRD = NBDDL * SEGDES XMATRI SEGINI MTRAV * *- Activation du segment contenant les valeurs *- du MCHAML Ck|K| pour la zone IMAIL * IF (IPCK.NE.0) THEN MCHAML = ICHAML(IMAIL) SEGACT MCHAML MELVAL = IELVAL(1) SEGACT MELVAL ENDIF * *- Convection * IF (IFLUX.NE.0) THEN * * Activation du segment contenant les valeurs du MCHAML d'orientation * des normales par face pour la zone IMAIL * MCHAM1 = MCHEL1.ICHAML(IMAIL) SEGACT MCHAM1 MELVA1 = MCHAM1.IELVAL(1) SEGACT MELVA1 * * Initialisation du chapeau de rigidité pour la zone IMAIL * IF (LHS1.NE.0) THEN RI1.COERIG(IMAIL) = 1.D0 RI1.IRIGEL(1,IMAIL) = IRIGEL(1,IMAIL) RI1.IRIGEL(3,IMAIL) = IRIGEL(3,IMAIL) SEGINI xMATR1 RI1.IRIGEL(4,IMAIL) = xMATR1 RI1.IRIGEL(7,IMAIL) = 2 ENDIF ENDIF * *------------------------------------------------------- * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL *------------------------------------------------------- * DO 100 IEL=1,NBELEM ITELEM = ITELEM + 1 * *- Orientation du flux et calcul de la source = -div(u*Th) * DIVUTH = 0.D0 IF (IFLUX.NE.0) THEN DO 20 I=1,NBDDL IPOPTS = ICPR(NUM(I,IEL)) FLUORI(I)= MPOVA4.VPOCHA(IPOPTS,1)*MELVA1.VELCHE(I,IEL) DIVUTH = DIVUTH - (FLUORI(I)*MPOVA2.VPOCHA(IPOPTS,1)) 20 CONTINUE ENDIF * *- Calcul de la somme des coefs pour une ligne ; une colonne *- -1 t -1 *- LIGNE = RE * DIV ; COLON = DIV * RE *- -1 t *- Calcul de CONSD = DIV * RE * DIV * CONSD = 0.D0 * XMATRI = IMATTT(IEL) * SEGACT XMATRI DO 40 I=1,NBDDL RCOLON(I) = 0.D0 RLIGNE(I) = 0.D0 DO 30 J=1,NBDDL RCOLON(I) = RCOLON(I) + RE(J,I,iel) RLIGNE(I) = RLIGNE(I) + RE(I,J,iel) 30 CONTINUE CONSD = CONSD + RLIGNE(I) 40 CONTINUE * SEGDES XMATRI * *- Recuperation de la valeur SOUR au centre de l'element *- ET ajout du terme de convection éventuelle (1-thetac)*(-div(u*Th)) * IF (IPCK.EQ.0) THEN VSOU = MPOVA1.VPOCHA(ITELEM,1) ELSE IF (ISOUR.NE.0) THEN VSOU = MPOVA1.VPOCHA(ITELEM,1)+((1.D0-THETAC)*DIVUTH) ELSE VSOU = (1.D0-THETAC) * DIVUTH ENDIF ENDIF * * IF (IPCK.EQ.0) THEN IF (RHS1.NE.0) THEN DO 45 IDDL=1,NBDDL VALFLU = RLIGNE(IDDL) * VSOU / CONSD IPOPTS = ICPR(NUM(IDDL,IEL)) VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - VALFLU 45 CONTINUE ENDIF IF (LHS1.NE.0) THEN * SEGINI XMATR1 * IMATR1.IMATTT(IEL) = XMATR1 DO 55 J=1,NBDDL REMU1 = -1.D0 / CONSD DO 50 I=1,NBDDL XMATR1.RE(I,J,iel) = REMU1*RLIGNE(I)*FLUORI(J) 50 CONTINUE 55 CONTINUE * SEGDES XMATR1 ENDIF ELSE BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT) VH = MPOVA3.VPOCHA(ITELEM,1) DO 70 IDDL=1,NBDDL SMFLU = RLIGNE(IDDL) * VSOU * BETA/CONSD SMP = RLIGNE(IDDL) * VH * (1.D0-BETA) STP = 0.D0 DO 60 J=1,NBDDL VTH = MPOVA2.VPOCHA(ICPR(NUM(J,IEL)),1) STP = STP + RCOLON(J) * VTH 60 CONTINUE SMTP = RLIGNE(IDDL) * STP * (1.D0-THETA)*BETA/CONSD SMDDL = SMFLU + SMP + SMTP IPOPTS = ICPR(NUM(IDDL,IEL)) VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - SMDDL 70 CONTINUE IF (LHS1.NE.0) THEN * SEGINI XMATR1 * IMATR1.IMATTT(IEL) = XMATR1 REMU1 = (-1.D0) * THETAC * BETA / CONSD DO 90 J=1,NBDDL DO 80 I=1,NBDDL XMATR1.RE(I,J,iel)=REMU1*RLIGNE(I) * FLUORI(J) 80 CONTINUE 90 CONTINUE * SEGDES XMATR1 ENDIF ENDIF 100 CONTINUE SEGDES xMATRI SEGDES MELEME IF (IPCK.NE.0) THEN SEGDES MCHAML SEGDES MELVAL ENDIF IF (IFLUX.NE.0) THEN SEGDES MCHAM1 SEGDES MELVA1 IF (LHS1.NE.0) SEGDES xMATR1 ENDIF SEGSUP MTRAV 110 CONTINUE * *- Desactivation des segments * IF (RHS1.NE.0) THEN SEGDES MPOVAL ENDIF IF (IPCK.NE.0) THEN SEGDES MCHELM SEGDES MPOVA2 SEGDES MPOVA3 ENDIF IF (IFLUX.NE.0) THEN SEGDES MCHEL1 SEGDES MPOVA4 IF (LHS1.NE.0) SEGDES RI1 ENDIF IF (ISOUR.NE.0) THEN SEGDES MPOVA1 ENDIF SEGDES IPMAHY SEGDES MRIGID SEGSUP ICCPR * RETURN END