1 : C MATP SOURCE CHAT 11/03/16 21:27:48 6902 2 : SUBROUTINE MATP 3 : C----------------------------------------------------------------------- 4 : C Cette subroutine permet de creer les matrices elementaires du système 5 : C matriciel en trace de charge dans le cadre de la résolution des 6 : C équations de DARCY par une méthode d'éléments finis mixtes hybride. 7 : C On obtient ce systeme en exprimant les inconnues de charge et de 8 : C flux en fonction des traces de charge. 9 : C -1 t -1 10 : C M1 = M2 * ( COEF * DIV * DIV * M2 - Id ) 11 : C 12 : C avec M2 : Matrice massse hybride 13 : C DIV : Matrice uniligne representant la divergence 14 : C COEF : Scalaire valant 1/d dans le cas permanant, THETA*BETA/d 15 : C dans le cas transitoire. 16 : C -1 t 17 : C d : Scalaire valant ( DIV * M2 * DIV ) 18 : C 19 : C THETA : Scalaire de la THETA-méthode | Données 20 : C BETA : Scalaire d*dt/(Ck|K| + THETA*d*dt) | utiles 21 : C dt : Scalaire pas de temps | uniquement 22 : C Ck : Scalaire coefficient d'emmagasinement | en 23 : C |K| : Scalaire aire de l'élément K | transitoire 24 : C 25 : C--------------------------- 26 : C Phrase d'appel (GIBIANE) : 27 : C--------------------------- 28 : C 29 : C MRIGI1 = MATP MMODEL MRIGI2 (TABLE2) ; 30 : C 31 : C------------------------ 32 : C Operandes et resultat : 33 : C------------------------ 34 : C 35 : C MRIGI1 : Matrices elementaires en trace de charge. 36 : C MMODEL : Objet modele specifiant la formulation. 37 : C MRIGI2 : Matrices masses hybrides elementaires creees par MHYB. 38 : C TABLE2 : TABLE DARCY_TRANSITOIRE contenant les infos transitoires. 39 : C 40 : C----------------------------------------------------------------------- 41 : C 42 : C Langage : ESOPE + FORTRAN77 43 : C 44 : C Auteurs : 09/93 F.DABBENE - Cas permanent 45 : C 09/94 X.NOUVELLON - Extension au cas transitoire 46 : C 47 : C----------------------------------------------------------------------- 48 : IMPLICIT INTEGER(I-N) 49 : IMPLICIT REAL*8 (A-H,O-Z) 50 : * 51 : -INC CCOPTIO 52 : -INC SMELEME 53 : -INC SMRIGID 54 : -INC SMMODEL 55 : -INC SMTABLE 56 : -INC SMCHAML 57 : * 58 : SEGMENT IPMAHY 59 : INTEGER MAHYBR(NSOUS) 60 : ENDSEGMENT 61 : * 62 : LOGICAL LOGRE,LOGIN 63 : CHARACTER*6 CHAR6 64 : CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE 65 : * 66 : * Initialisations 67 : * 68 : IVALIN = 0 69 : XVALIN = 0.D0 70 : LOGIN = .TRUE. 71 : IOBIN = 0 72 : TAPIND = 'MOT ' 73 : * 74 : * Lecture du MMODEL 75 : * 75 : CALL LIROBJ('MMODEL',IPMODE,1,IRET1) 76 : IF (IERR.NE.0) RETURN 77 : MMODEL = IPMODE 78 : * 79 : * Lecture de la TABLE domaine et récupération des infos 80 : * 81 : IPTABL = 0 81 : CALL LEKMOD(IPMODE,IPTABL,IRET) 82 : IF (IERR.NE.0) RETURN 83 : CHARIN = 'MAILLAGE' 84 : TYPOBJ = 'MAILLAGE' 84 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 85 : IF (IERR.NE.0) RETURN 86 : IPGEOM = IOBRE 86 : CALL LEKTAB(IPTABL,'ELTFA',IOBRE) 87 : IF (IERR.NE.0) RETURN 88 : IELTFA = IOBRE 89 : * 90 : * Lecture de RIGIDITE 91 : * 91 : CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET1) 92 : IF (IERR.NE.0) RETURN 93 : MRIGID = IPRIGI 94 : * 95 : * Lecture éventuelle de la TABLE DARCY_TRANSITOIRE 96 : * 97 : IPTABL = 0 97 : CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET) 98 : IF (IERR.NE.0) RETURN 99 : IF (IRET.EQ.1) THEN 100 : TYPOBJ = 'FLOTTANT' 100 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN, 101 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 102 : IF (IERR.NE.0) RETURN 103 : THETA = XVALRE 104 : THEMIN = -1.D-12 105 : THEMAX = 1.D0 - THEMIN 106 : IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN 107 : REAERR(1) = REAL(THETA) 108 : REAERR(2) = REAL(0.D0) 109 : REAERR(3) = REAL(1.D0) 110 : MOTERR(1:8) = ' THETA ' 110 : CALL ERREUR(42) 111 : RETURN 112 : ENDIF 113 : IF (THETA.LT.0.D0) THETA=0.D0 114 : IF (THETA.GT.1.D0) THETA=1.D0 114 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN, 115 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 116 : IF (IERR.NE.0) RETURN 117 : DELTAT = XVALRE 118 : IF (DELTAT.LT.0.D0) THEN 119 : REAERR(1) = REAL(DELTAT) 120 : REAERR(2) = REAL(0.D0) 121 : MOTERR(1:8) = ' DELTAT ' 121 : CALL ERREUR(41) 122 : RETURN 123 : ENDIF 124 : TYPOBJ = 'MCHAML ' 124 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN, 125 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 126 : IF (IERR.NE.0) RETURN 127 : IPCK = IOBRE 128 : ELSE 129 : THETA = 0.D0 130 : DELTAT = 0.D0 131 : IPCK = 0 132 : ENDIF 133 : * 134 : * Test de l'ojet MODELE 135 : * Récupération des zones élémentaires ou DARCY est défini 136 : * 137 : SEGACT MMODEL 138 : NSOUS = KMODEL(/1) 139 : SEGINI IPMAHY 140 : IDARCY = 0 141 : IPT2=IPGEOM 142 : IPT1=IPGEOM 143 : DO 10 ISOUS=1,NSOUS 144 : IF(NSOUS.GT.1)THEN 145 : SEGACT IPT2 146 : IPT1=IPT2.LISOUS(ISOUS) 147 : SEGDES IPT2 148 : ENDIF 149 : IMODEL = KMODEL(ISOUS) 150 : SEGACT IMODEL 151 : LETYPE = FORMOD(1) 152 : IF (LETYPE.EQ.'DARCY') THEN 153 : IDARCY = IDARCY + 1 154 : MAHYBR(ISOUS) = IPT1 155 : ENDIF 156 : SEGDES IMODEL 157 : 10 CONTINUE 158 : SEGDES MMODEL 159 : IF (IDARCY.EQ.0) THEN 160 : MOTERR = LETYPE 160 : CALL ERREUR(193) 161 : GOTO 100 162 : ENDIF 163 : * 164 : * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini 165 : * 166 : MELEME = IELTFA 167 : SEGACT MELEME 168 : LZONES = LISOUS(/1) 169 : IF (LZONES.EQ.0) LZONES = 1 170 : IPT1 = IPGEOM 171 : SEGACT IPT1 172 : DO 30 ISOUS=1,NSOUS 173 : IMAMEL = MAHYBR(ISOUS) 174 : IF (IMAMEL.NE.0) THEN 175 : DO 20 ISZ=1,LZONES 176 : IF (LZONES.EQ.1) THEN 177 : IPT2 = IPT1 178 : IPT3 = MELEME 179 : ELSE 180 : IPT2 = IPT1.LISOUS(ISZ) 181 : IPT3 = LISOUS(ISZ) 182 : ENDIF 183 : IF (IPT2.EQ.IMAMEL) THEN 184 : MAHYBR(ISOUS) = IPT3 185 : GOTO 30 186 : ENDIF 187 : 20 CONTINUE 188 : IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN 189 : MOTERR(1:8) = ' MODELE ' 190 : MOTERR(9:16) = ' ELTFA ' 191 : INTERR(1) = ISOUS 191 : CALL ERREUR(698) 192 : SEGDES IPT1 193 : SEGDES MELEME 194 : GOTO 100 195 : ENDIF 196 : ENDIF 197 : 30 CONTINUE 198 : SEGDES IPT1 199 : SEGDES MELEME 200 : * 201 : * Test du sous-type de la matrice de rigidité récupérée 202 : * Controle des MELEME de l'objet RIGIDITE % au ELTFA liée au MMODEL 203 : * 204 : SEGACT MRIGID 205 : LETYPE = MTYMAT 206 : IF (LETYPE.NE.'DARCY') THEN 207 : MOTERR(1:8) = 'RIGIDITE' 208 : MOTERR(9:16) = 'DARCY' 208 : CALL ERREUR(79) 209 : SEGDES MRIGID 210 : GOTO 100 211 : ENDIF 212 : DO 40 ISOUS=1,NSOUS 213 : IF (MAHYBR(ISOUS).NE.0) THEN 214 : IPTEST = IRIGEL(1,ISOUS) 215 : IF (IPTEST.NE.MAHYBR(ISOUS)) THEN 216 : MOTERR(1:8) = 'DARCY' 217 : MOTERR(9:16) = ' ELTFA ' 218 : INTERR(1) = ISOUS 218 : CALL ERREUR(698) 219 : SEGDES MRIGID 220 : GOTO 100 221 : ENDIF 222 : ENDIF 223 : 40 CONTINUE 224 : SEGDES MRIGID 225 : * 226 : * Vérification du support du MCHAML % au MMODEL 227 : * Controle du nombre de composantes reelles 228 : * 229 : IF (IPCK.NE.0) THEN 230 : ISUP = 2 231 : ICOND = 1 231 : CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2) 232 : IF (IRET.GT.1) GOTO 100 233 : MCHELM = IPCK 234 : SEGACT MCHELM 235 : DO 50 ISOUS=1,NSOUS 236 : IF (MAHYBR(ISOUS).NE.0) THEN 237 : MCHAML = ICHAML(ISOUS) 238 : SEGACT MCHAML 239 : N2 = TYPCHE(/2) 240 : IF (N2.GT.1) THEN 240 : CALL ERREUR(320) 241 : SEGDES MCHELM 242 : SEGDES MCHAML 243 : GOTO 100 244 : ENDIF 245 : CHAR6 = TYPCHE(1)(1:6) 246 : IF (CHAR6.NE.'REAL*8') THEN 247 : MOTERR(1:8) = NOMCHE(1) 247 : CALL ERREUR(679) 248 : SEGDES MCHELM 249 : SEGDES MCHAML 250 : GOTO 100 251 : ENDIF 252 : SEGDES MCHAML 253 : ENDIF 254 : 50 CONTINUE 255 : SEGDES MCHELM 256 : ENDIF 257 : SEGDES IPMAHY 258 : * 259 : * Construction de IPRIG1 260 : * 260 : CALL MATP1(NSOUS,IPGEOM,IPMAHY,IPRIGI,THETA,DELTAT,IPCK,IPRIG1) 260 : CALL ECROBJ('RIGIDITE',IPRIG1) 261 : * 262 : * Ménage 263 : * 264 : 100 CONTINUE 265 : SEGSUP IPMAHY 266 : RETURN 267 : END 268 : 269 : 270 : 271 : 272 : 273 : 274 : 275 : 276 : 277 : 278 : 279 : 280 : 281 :
© Cast3M 2003 - All rights reserved.
Disclaimer