C MATP SOURCE CHAT 11/03/16 21:27:48 6902 SUBROUTINE MATP C----------------------------------------------------------------------- C Cette subroutine permet de creer les matrices elementaires du système C matriciel en trace de charge dans le cadre de la résolution des C équations de DARCY par une méthode d'éléments finis mixtes hybride. C On obtient ce systeme en exprimant les inconnues de charge et de C flux en fonction des traces de charge. C -1 t -1 C M1 = M2 * ( COEF * DIV * DIV * M2 - Id ) C C avec M2 : Matrice massse hybride C DIV : Matrice uniligne representant la divergence C COEF : Scalaire valant 1/d dans le cas permanant, THETA*BETA/d C dans le cas transitoire. C -1 t C d : Scalaire valant ( DIV * M2 * DIV ) C C THETA : Scalaire de la THETA-méthode | Données C BETA : Scalaire d*dt/(Ck|K| + THETA*d*dt) | utiles C dt : Scalaire pas de temps | uniquement C Ck : Scalaire coefficient d'emmagasinement | en C |K| : Scalaire aire de l'élément K | transitoire C C--------------------------- C Phrase d'appel (GIBIANE) : C--------------------------- C C MRIGI1 = MATP MMODEL MRIGI2 (TABLE2) ; C C------------------------ C Operandes et resultat : C------------------------ C C MRIGI1 : Matrices elementaires en trace de charge. C MMODEL : Objet modele specifiant la formulation. C MRIGI2 : Matrices masses hybrides elementaires creees par MHYB. C TABLE2 : TABLE DARCY_TRANSITOIRE contenant les infos transitoires. 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 C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC CCOPTIO -INC SMELEME -INC SMRIGID -INC SMMODEL -INC SMTABLE -INC SMCHAML * SEGMENT IPMAHY INTEGER MAHYBR(NSOUS) ENDSEGMENT * LOGICAL LOGRE,LOGIN CHARACTER*6 CHAR6 CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE * * Initialisations * IVALIN = 0 XVALIN = 0.D0 LOGIN = .TRUE. IOBIN = 0 TAPIND = 'MOT ' * * Lecture du MMODEL * CALL LIROBJ('MMODEL',IPMODE,1,IRET1) IF (IERR.NE.0) RETURN MMODEL = IPMODE * * Lecture de la TABLE domaine et récupération des infos * IPTABL = 0 CALL LEKMOD(IPMODE,IPTABL,IRET) IF (IERR.NE.0) RETURN CHARIN = 'MAILLAGE' TYPOBJ = 'MAILLAGE' CALL LEKTAB(IPTABL,CHARIN,IOBRE) IF (IERR.NE.0) RETURN IPGEOM = IOBRE CALL LEKTAB(IPTABL,'ELTFA',IOBRE) IF (IERR.NE.0) RETURN IELTFA = IOBRE * * Lecture de RIGIDITE * CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET1) IF (IERR.NE.0) RETURN MRIGID = IPRIGI * * Lecture éventuelle de la TABLE DARCY_TRANSITOIRE * IPTABL = 0 CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET) IF (IERR.NE.0) RETURN IF (IRET.EQ.1) THEN TYPOBJ = 'FLOTTANT' CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN, S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) IF (IERR.NE.0) RETURN THETA = XVALRE THEMIN = -1.D-12 THEMAX = 1.D0 - THEMIN IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN REAERR(1) = REAL(THETA) REAERR(2) = REAL(0.D0) REAERR(3) = REAL(1.D0) MOTERR(1:8) = ' THETA ' CALL ERREUR(42) RETURN ENDIF IF (THETA.LT.0.D0) THETA=0.D0 IF (THETA.GT.1.D0) THETA=1.D0 CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN, S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) IF (IERR.NE.0) RETURN DELTAT = XVALRE IF (DELTAT.LT.0.D0) THEN REAERR(1) = REAL(DELTAT) REAERR(2) = REAL(0.D0) MOTERR(1:8) = ' DELTAT ' CALL ERREUR(41) RETURN ENDIF TYPOBJ = 'MCHAML ' CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN, S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) IF (IERR.NE.0) RETURN IPCK = IOBRE ELSE THETA = 0.D0 DELTAT = 0.D0 IPCK = 0 ENDIF * * Test de l'ojet MODELE * Récupération des zones élémentaires ou DARCY est défini * SEGACT MMODEL NSOUS = KMODEL(/1) SEGINI IPMAHY IDARCY = 0 IPT2=IPGEOM IPT1=IPGEOM DO 10 ISOUS=1,NSOUS IF(NSOUS.GT.1)THEN SEGACT IPT2 IPT1=IPT2.LISOUS(ISOUS) SEGDES IPT2 ENDIF IMODEL = KMODEL(ISOUS) SEGACT IMODEL LETYPE = FORMOD(1) IF (LETYPE.EQ.'DARCY') THEN IDARCY = IDARCY + 1 MAHYBR(ISOUS) = IPT1 ENDIF SEGDES IMODEL 10 CONTINUE SEGDES MMODEL IF (IDARCY.EQ.0) THEN MOTERR = LETYPE CALL ERREUR(193) GOTO 100 ENDIF * * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini * MELEME = IELTFA SEGACT MELEME LZONES = LISOUS(/1) IF (LZONES.EQ.0) LZONES = 1 IPT1 = IPGEOM SEGACT IPT1 DO 30 ISOUS=1,NSOUS IMAMEL = MAHYBR(ISOUS) IF (IMAMEL.NE.0) THEN DO 20 ISZ=1,LZONES IF (LZONES.EQ.1) THEN IPT2 = IPT1 IPT3 = MELEME ELSE IPT2 = IPT1.LISOUS(ISZ) IPT3 = LISOUS(ISZ) ENDIF IF (IPT2.EQ.IMAMEL) THEN MAHYBR(ISOUS) = IPT3 GOTO 30 ENDIF 20 CONTINUE IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN MOTERR(1:8) = ' MODELE ' MOTERR(9:16) = ' ELTFA ' INTERR(1) = ISOUS CALL ERREUR(698) SEGDES IPT1 SEGDES MELEME GOTO 100 ENDIF ENDIF 30 CONTINUE SEGDES IPT1 SEGDES MELEME * * Test du sous-type de la matrice de rigidité récupérée * Controle des MELEME de l'objet RIGIDITE % au ELTFA liée au MMODEL * SEGACT MRIGID LETYPE = MTYMAT IF (LETYPE.NE.'DARCY') THEN MOTERR(1:8) = 'RIGIDITE' MOTERR(9:16) = 'DARCY' CALL ERREUR(79) SEGDES MRIGID GOTO 100 ENDIF DO 40 ISOUS=1,NSOUS IF (MAHYBR(ISOUS).NE.0) THEN IPTEST = IRIGEL(1,ISOUS) IF (IPTEST.NE.MAHYBR(ISOUS)) THEN MOTERR(1:8) = 'DARCY' MOTERR(9:16) = ' ELTFA ' INTERR(1) = ISOUS CALL ERREUR(698) SEGDES MRIGID GOTO 100 ENDIF ENDIF 40 CONTINUE SEGDES MRIGID * * Vérification du support du MCHAML % au MMODEL * Controle du nombre de composantes reelles * IF (IPCK.NE.0) THEN ISUP = 2 ICOND = 1 CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2) IF (IRET.GT.1) GOTO 100 MCHELM = IPCK SEGACT MCHELM DO 50 ISOUS=1,NSOUS IF (MAHYBR(ISOUS).NE.0) THEN MCHAML = ICHAML(ISOUS) SEGACT MCHAML N2 = TYPCHE(/2) IF (N2.GT.1) THEN CALL ERREUR(320) SEGDES MCHELM SEGDES MCHAML GOTO 100 ENDIF CHAR6 = TYPCHE(1)(1:6) IF (CHAR6.NE.'REAL*8') THEN MOTERR(1:8) = NOMCHE(1) CALL ERREUR(679) SEGDES MCHELM SEGDES MCHAML GOTO 100 ENDIF SEGDES MCHAML ENDIF 50 CONTINUE SEGDES MCHELM ENDIF SEGDES IPMAHY * * Construction de IPRIG1 * CALL MATP1(NSOUS,IPGEOM,IPMAHY,IPRIGI,THETA,DELTAT,IPCK,IPRIG1) CALL ECROBJ('RIGIDITE',IPRIG1) * * Ménage * 100 CONTINUE SEGSUP IPMAHY RETURN END