1 : C DMMU SOURCE CHAT 11/03/16 21:19:16 6902 2 : SUBROUTINE DMMU 3 : C----------------------------------------------------------------------- 4 : C -1 t 5 : C Calcul du produit D (M -U) B CHPface 6 : C -1 t 7 : C ou du produit D M B CHPface 8 : C le résultat est un CHPcentre 9 : C----------------------------------------------------------------------- 10 : C 11 : C--------------------------- 12 : C Phrase d'appel (GIBIANE) : 13 : C--------------------------- 14 : C 15 : C CHP3 = 'DMMU' MMODEL RIG1 (CHP1) CHP2 ; 16 : C 17 : C------------------------ 18 : C Operandes et resultat : 19 : C------------------------ 20 : C 21 : C MMODEL : modèle formulation DARCY. 22 : C RIG1 : Matrices hybrides elementaires de DARCY crees par MHYB. 23 : C CHP1 : flux aux faces 24 : C CHP2 : CHPO face à plusieur composantes. 25 : C CHP3 : CHPO centre à plusieur composantes dont le support géometrique 26 : C est le maillage centre du domaine les noms des composantes 27 : C sont ceux de CHP2 28 : C 29 : C----------------------------------------------------------------------- 30 : C 31 : C Langage : ESOPE + FORTRAN77 32 : C 33 : C 34 : C----------------------------------------------------------------------- 35 : IMPLICIT INTEGER(I-N) 36 : IMPLICIT REAL*8 (A-H,O-Z) 37 : * 38 : -INC CCOPTIO 39 : -INC SMCOORD 40 : -INC SMCHPOI 41 : -INC SMCHAML 42 : -INC SMELEME 43 : -INC SMRIGID 44 : -INC SMTABLE 45 : -INC SMMODEL 46 : * 47 : SEGMENT ICCPR 48 : INTEGER ICPR(NNGOT) 49 : ENDSEGMENT 50 : SEGMENT IORGA 51 : INTEGER IAA(ITES) 52 : ENDSEGMENT 53 : * 54 : LOGICAL LOGRE,LOGIN 55 : CHARACTER*8 TAPIND,TYPOBJ,CHARRE,LETYPE 56 : CHARACTER*4 NOMTOT(1) 57 : * 58 : * Initialisations 59 : * 60 : IVALIN = 0 61 : XVALIN = 0.D0 62 : LOGIN = .TRUE. 63 : IOBIN = 0 64 : TAPIND = 'MOT ' 65 : * 66 : * 67 : * Lecture du MMODEL 68 : * 68 : CALL LIROBJ('MMODEL',IPMODE,1,IRET) 69 : IF(IERR.NE.0)RETURN 70 : MMODEL=IPMODE 71 : SEGACT MMODEL 72 : N1=KMODEL(/1) 73 : DO 7 I=1,N1 74 : IMODEL=KMODEL(I) 75 : SEGACT IMODEL 76 : IF(FORMOD(1).NE.'DARCY')THEN 77 : MOTERR(1:16) = 'DARCY ' 77 : CALL ERREUR(719) 78 : RETURN 79 : ENDIF 80 : 7 CONTINUE 81 : C on récupère la table DOMAINE 82 : IPTABL = 0 82 : CALL LEKMOD(MMODEL,IPTABL,IRET) 83 : IF (IERR.NE.0) RETURN 84 : TYPOBJ = 'MAILLAGE' 84 : CALL LEKTAB(IPTABL,'ELTFA',IOBRE) 85 : IF (IERR.NE.0) RETURN 86 : IELTFA = IOBRE 86 : CALL LEKTAB(IPTABL,'CENTRE',IOBRE) 87 : IF (IERR.NE.0) RETURN 88 : ICENTR = IOBRE 88 : CALL LEKTAB(IPTABL,'FACE',IOBRE) 89 : IF (IERR.NE.0) RETURN 90 : IFACE = IOBRE 91 : * 92 : * Lecture de RIGIDITE 93 : * 93 : CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) 94 : IF (IERR.NE.0) RETURN 95 : MRIGID = IPRIGI 96 : * 97 : * 98 : * 99 : * Test du sous-type de la matrice de rigiditée récupérée 100 : * 101 : SEGACT MRIGID 102 : LETYPE = MTYMAT 103 : IF (LETYPE.NE.'DARCY') THEN 104 : MOTERR(1:8) = 'RIGIDITE' 105 : MOTERR(9:16) = 'DARCY ' 105 : CALL ERREUR(79) 106 : SEGDES MRIGID 107 : GOTO 100 108 : ENDIF 109 : * 110 : * Controle des pointeurs de MELEME de la rigidité 111 : * 112 : NRIGEL=IRIGEL(/2) 113 : C write(6,*)' NRIGEL----- ' ,nrigel 114 : MELEME=IELTFA 115 : SEGACT MELEME 116 : NBSOUS=LISOUS(/1) 117 : IF(NBSOUS.EQ.0)THEN 118 : IF((NRIGEL.NE.1).OR.(IRIGEL(1,1).NE.MELEME))THEN 119 : MOTERR(1:8) = 'DARCY ' 120 : MOTERR(9:16) = 'ELTFA ' 121 : INTERR(1) = 1 121 : CALL ERREUR(698) 122 : SEGDES MRIGID 123 : GOTO 100 124 : ENDIF 125 : ELSE 126 : IF(NRIGEL.NE.NBSOUS)THEN 127 : MOTERR(1:8) = 'DARCY ' 128 : MOTERR(9:16) = 'ELTFA ' 129 : INTERR(1) = 1 129 : CALL ERREUR(698) 130 : SEGDES MRIGID 131 : GOTO 100 132 : ENDIF 133 : DO 10 ISOUS=1,NBSOUS 134 : IF (LISOUS(ISOUS).NE.IRIGEL(1,ISOUS)) THEN 135 : MOTERR(1:8) = 'DARCY ' 136 : MOTERR(9:16) = 'ELTFA ' 137 : INTERR(1) = ISOUS 137 : CALL ERREUR(698) 138 : SEGDES MRIGID 139 : GOTO 100 140 : ENDIF 141 : 10 CONTINUE 142 : ENDIF 143 : * 144 : * LECTURE DES FLUX OU 145 : * DU CHPOIN FACE A MULTIPLIER 146 : * 146 : CALL LIROBJ('CHPOINT',IPCHF,1,IRET) 147 : IF(IRET.NE.1) RETURN 148 : MCHPO1=IPCHF 149 : NOMTOT(1)=' ' 150 : INDIC=0 151 : NBCOMP=-1 151 : CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT) 152 : IF(IERR.NE.0)RETURN 153 : SEGACT MCHPO1 154 : MSOUP1=MCHPO1.IPCHP(1) 155 : SEGACT MSOUP1 156 : IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN 157 : MPOVA1=MSOUP1.IPOVAL 158 : SEGACT MPOVA1 159 : * 160 : * 160 : CALL LIROBJ('CHPOINT',IPCHF,1,IRET) 161 : IF(IRET.NE.1) RETURN 162 : MCHPO2=IPCHF 163 : NOMTOT(1)=' ' 164 : INDIC=0 165 : NBCOMP=-1 165 : CALL QUEPOI(MCHPO2,IFACE,INDIC,NBCOMP,NOMTOT) 166 : IF(IERR.NE.0)RETURN 167 : SEGACT MCHPO2 168 : MSOUP2=MCHPO2.IPCHP(1) 169 : SEGACT MSOUP2 170 : MPOVA2=MSOUP2.IPOVAL 171 : SEGACT MPOVA2 172 : ELSE 173 : MCHPO2=MCHPO1 174 : MCHPO1=0 175 : MSOUP2=MSOUP1 176 : MPOVA2=MSOUP2.IPOVAL 177 : SEGACT MPOVA2 177 : CALL LIROBJ('CHPOINT',IPCHF,0,IRET) 178 : IF(IRET.EQ.1) THEN 179 : MCHPO1=IPCHF 180 : NOMTOT(1)=' ' 181 : INDIC=0 182 : NBCOMP=-1 182 : CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT) 183 : IF(IERR.NE.0)RETURN 184 : SEGACT MCHPO1 185 : MSOUP1=MCHPO1.IPCHP(1) 186 : SEGACT MSOUP1 187 : IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN 188 : MPOVA1=MSOUP1.IPOVAL 189 : SEGACT MPOVA1 190 : ELSE 190 : CALL ERREUR(21) 191 : RETURN 192 : ENDIF 193 : ENDIF 194 : ENDIF 195 : * 196 : * Construction de MCHPOI 197 : * 198 : * 199 : IPT1=ICENTR 200 : SEGACT IPT1 201 : NPN=IPT1.NUM(/2) 202 : NSOUPO=1 203 : NAT=1 204 : SEGINI MCHPOI 205 : MTYPOI=' ' 206 : MOCHDE=' CHPOIN CREE PAR DMMU ' 207 : IFOPOI=IFOUR 208 : JATTRI(1)=2 209 : NC=NBCOMP 210 : SEGINI MSOUPO 211 : IPCHP(1)=MSOUPO 212 : DO 5 L=1,NBCOMP 213 : NOCOMP(L)=MSOUP2.NOCOMP(L) 214 : NOHARM(L)=MSOUP2.NOHARM(L) 215 : 5 CONTINUE 216 : IGEOC=ICENTR 217 : N=NPN 218 : SEGINI MPOVAL 219 : IPOVAL=MPOVAL 220 : NB=N*NC 220 : CALL INITD(VPOCHA,NB,0.D0) 221 : * 222 : * Creation du tableau ICPR 223 : * 224 : IK = 0 225 : NNGOT = XCOOR(/1)/(IDIM+1) 226 : SEGINI ICCPR 227 : MELEME = IFACE 228 : SEGACT MELEME 229 : N2 = NUM(/2) 230 : DO 15 I2=1,N2 231 : K = NUM(1,I2) 232 : IF (ICPR(K).EQ.0) THEN 233 : IK = IK + 1 234 : ICPR(K) = IK 235 : ENDIF 236 : 15 CONTINUE 237 : ITES=IK 238 : SEGDES MELEME 239 : SEGINI IORGA 240 : MELEME = MSOUP2.IGEOC 241 : SEGACT MELEME 242 : N2 = NUM(/2) 243 : DO 25 I2=1,N2 244 : K = NUM(1,I2) 245 : IF (ICPR(K).EQ.0) THEN 246 : INTERR(1) = K 247 : MOTERR(1:8) = 'FACE ' 247 : CALL ERREUR(64) 248 : SEGSUP ICCPR, IORGA 249 : RETURN 250 : ELSE 251 : IAA(ICPR(K)) = I2 252 : ENDIF 253 : 25 CONTINUE 254 : SEGDES MELEME 254 : CALL LEKTAB(IPTABL,'XXNORMAE',IORIE) 255 : MCHELM=IORIE 256 : SEGACT MCHELM 257 : C 258 : C Calcul du produit 259 : C 260 : MELEME=IELTFA 261 : SEGACT MELEME 262 : IPT1=MELEME 263 : ITELEM=0 264 : IF(MCHPO1.NE.0)THEN 265 : C 266 : C Cas avec vitesses 267 : C 268 : DO 50 ISOUS=1,NRIGEL 269 : MCHAML= ICHAML(ISOUS) 270 : SEGACT MCHAML 271 : MELVAL=IELVAL(1) 272 : SEGACT MELVAL 273 : IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS) 274 : SEGACT IPT1 275 : xMATRI=IRIGEL(4,ISOUS) 276 : SEGACT xMATRI 277 : NELRIG=re(/3) 278 : C write(6,*)' ISOUS NELRIG ',ISOUS,NELRIG 279 : DO 60 IEL=1,NELRIG 280 : ITELEM=ITELEM+1 281 : * XMATRI=IMATTT(IEL) 282 : * SEGACT XMATRI 283 : NLIGRD=RE(/1) 284 : NLIGRP=RE(/2) 285 : C write(6,*) ' NLIGRD NLIGRP ',NLIGRD,NLIGRP 286 : DO 40 I=1,NLIGRD 287 : RLIGN=0.D0 288 : DO 30 J=1,NLIGRP 289 : RLIGN=RLIGN+RE(I,J,iel) 290 : 30 CONTINUE 291 : IPOPTS=ICPR(IPT1.NUM(I,IEL)) 292 : III=IAA(IPOPTS) 293 : C write(6,*)' dmmu ',itelem,i,rlign,MPOVA1.VPOCHA(IPOPTS,1) 294 : C * , VELCHE(I,IEL) 295 : RLIGN=RLIGN-MPOVA1.VPOCHA(IPOPTS,1)*VELCHE(I,IEL) 296 : DO 20 K=1,NBCOMP 297 : VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN* 298 : * MPOVA2.VPOCHA(III,K) 299 : 20 CONTINUE 300 : 301 : 40 CONTINUE 302 : * SEGDES XMATRI 303 : 60 CONTINUE 304 : SEGDES xMATRI,MCHAML,MELVAL 305 : 50 CONTINUE 306 : ELSE 307 : C 308 : C Cas ou il n'y a pas de vitesse 309 : C On a dupliqué toute la boucle 50 pour ne pas faire 310 : C NRIGEL*NELRIG*NLIGRD fois le meme test 311 : C 312 : DO 51 ISOUS=1,NRIGEL 313 : MCHAML= ICHAML(ISOUS) 314 : SEGACT MCHAML 315 : MELVAL=IELVAL(1) 316 : SEGACT MELVAL 317 : IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS) 318 : SEGACT IPT1 319 : xMATRI=IRIGEL(4,ISOUS) 320 : SEGACT xMATRI 321 : NELRIG=re(/3) 322 : DO 61 IEL=1,NELRIG 323 : ITELEM=ITELEM+1 324 : * XMATRI=IMATTT(IEL) 325 : * SEGACT XMATRI 326 : NLIGRD=RE(/1) 327 : NLIGRP=RE(/2) 328 : DO 41 I=1,NLIGRD 329 : RLIGN=0.D0 330 : DO 31 J=1,NLIGRP 331 : RLIGN=RLIGN+RE(I,J,iel) 332 : 31 CONTINUE 333 : IPOPTS=ICPR(IPT1.NUM(I,IEL)) 334 : III=IAA(IPOPTS) 335 : DO 21 K=1,NBCOMP 336 : VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN* 337 : * MPOVA2.VPOCHA(III,K) 338 : 21 CONTINUE 339 : 340 : 41 CONTINUE 341 : * SEGDES XMATRI 342 : 61 CONTINUE 343 : SEGDES xMATRI,MCHAML,MELVAL 344 : 51 CONTINUE 345 : ENDIF 346 : SEGDES MRIGID 346 : CALL ECROBJ('CHPOINT',MCHPOI) 347 : * 348 : * Ménage 349 : * 350 : SEGDES MELEME 351 : SEGSUP ICCPR ,IORGA 352 : SEGDES MPOVAL,MSOUPO,MCHPOI 353 : SEGDES MPOVA2,MSOUP2,MCHPO2 354 : IF(MCHPO1.NE.0) THEN 355 : SEGDES MPOVA1,MSOUP1,MCHPO1 356 : ENDIF 357 : 100 CONTINUE 358 : RETURN 359 : END 360 : 361 : 362 : 363 : 364 : 365 : 366 : 367 : 368 : 369 : 370 : 371 : 372 : 373 : 374 : 375 : 376 : 377 :
© Cast3M 2003 - All rights reserved.
Disclaimer