1 : C CCGMUS SOURCE GOUNAND 06/12/08 21:15:18 5597 2 : SUBROUTINE CCGMUS(LCOF,NOMLOI, 3 : $ FC, 4 : $ IMPR,IRET) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : IMPLICIT INTEGER (I-N) 7 : C*********************************************************************** 8 : C NOM : CCGMUS 9 : C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : 10 : C Taille d'un élément suivant une direction donnée 11 : C 12 : C LANGAGE : ESOPE 13 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 14 : C mél : gounand@semt2.smts.cea.fr 15 : C*********************************************************************** 16 : C APPELES : 17 : C APPELE PAR : 18 : C*********************************************************************** 19 : C ENTREES : 20 : C ENTREES/SORTIES : 21 : C SORTIES : - 22 : C TRAVAIL : 23 : C*********************************************************************** 24 : C VERSION : v1, 30/08/06, version initiale 25 : C HISTORIQUE : v1, 30/08/06, création 26 : C HISTORIQUE : 27 : C HISTORIQUE : 28 : C*********************************************************************** 29 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 30 : C en cas de modification de ce sous-programme afin de faciliter 31 : C la maintenance ! 32 : C*********************************************************************** 33 : -INC CCOPTIO 34 : -INC CCREEL 35 : CBEGININCLUDE SMCHAEL 36 : SEGMENT MCHAEL 37 : POINTEUR IMACHE(N1).MELEME 38 : POINTEUR ICHEVA(N1).MCHEVA 39 : ENDSEGMENT 40 : SEGMENT MCHEVA 41 : REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM) 42 : ENDSEGMENT 43 : SEGMENT LCHEVA 44 : POINTEUR LISCHE(NBCHE).MCHEVA 45 : ENDSEGMENT 46 : CENDINCLUDE SMCHAEL 47 : INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1 48 : POINTEUR FC.MCHEVA 49 : POINTEUR LCOF.LCHEVA 50 : POINTEUR JMAJAC.MCHEVA 51 : POINTEUR JMIJAC.MCHEVA 52 : POINTEUR JDTJAC.MCHEVA 53 : POINTEUR JMAREG.MCHEVA 54 : POINTEUR JDIAMA.MCHEVA 55 : POINTEUR JPC.MCHEVA 56 : POINTEUR JVIT.MCHEVA 57 : POINTEUR JRHO.MCHEVA 58 : POINTEUR JMU.MCHEVA 59 : POINTEUR JPEC.MCHEVA 60 : CHARACTER*8 NOMLOI 61 : INTEGER ICOF 62 : * 63 : -INC TMXMAT 64 : POINTEUR A.MXMAT 65 : POINTEUR JMA.MXMAT 66 : POINTEUR KJMA.MXMAT 67 : POINTEUR JM1.MXMAT 68 : POINTEUR K.MXMAT 69 : * 70 : SEGMENT MVIT 71 : POINTEUR MVCOMP(IDIM).MCHEVA 72 : ENDSEGMENT 73 : LOGICAL LCRIT 74 : * 75 : INTEGER IMPR,IRET 76 : * 77 : * Executable statements 78 : * 79 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgmus' 80 : XPETI=SQRT(XPETIT) 80 : CALL CH2INT(NOMLOI(7:7),IMETH,IMPR,IRET) 81 : IF (IRET.NE.0) GOTO 9999 81 : CALL CH2INT(NOMLOI(8:8),IDIMI,IMPR,IRET) 82 : IF (IRET.NE.0) GOTO 9999 83 : NLFC=FC.VELCHE(/6) 84 : NPFC=FC.VELCHE(/5) 85 : ICOF=0 86 : * 87 : * Récupération de RHO et MU 88 : * 89 : ICOF=ICOF+1 90 : JRHO=LCOF.LISCHE(ICOF) 91 : NLRO=JRHO.VELCHE(/6) 92 : NPRO=JRHO.VELCHE(/5) 93 : ICOF=ICOF+1 94 : JMU=LCOF.LISCHE(ICOF) 95 : NLMU=JMU.VELCHE(/6) 96 : NPMU=JMU.VELCHE(/5) 97 : * 98 : * Récupération des coefficients de la metrique 99 : * 100 : ICOF1=ICOF+1 101 : SEGINI MVIT 102 : DO IIDIM=1,IDIM 103 : ICOF=ICOF+1 104 : JVIT=LCOF.LISCHE(ICOF) 105 : IF (ICOF.EQ.ICOF1) THEN 106 : NLJV=JVIT.VELCHE(/6) 107 : NPJV=JVIT.VELCHE(/5) 108 : ELSE 109 : NLJV2=JVIT.VELCHE(/6) 110 : NPJV2=JVIT.VELCHE(/5) 111 : IF (NLJV2.NE.NLJV.OR.NPJV2.NE.NPJV) THEN 112 : WRITE(IOIMP,*) 'Erreur grave dims JVIT' 113 : GOTO 9999 114 : ENDIF 115 : ENDIF 116 : MVIT.MVCOMP(IIDIM)=JVIT 117 : ENDDO 118 : * 119 : * Récupération du Péclet critique 120 : * 121 : ICOF=ICOF+1 122 : JPEC=LCOF.LISCHE(ICOF) 123 : NLPEC=JPEC.VELCHE(/6) 124 : NPPEC=JPEC.VELCHE(/5) 125 : * 126 : ICOF=ICOF+1 127 : JMAJAC=LCOF.LISCHE(ICOF) 128 : C NLJA=JMAJAC.VELCHE(/6) 129 : C NPJA=JMAJAC.VELCHE(/5) 130 : IREF=JMAJAC.VELCHE(/4) 131 : IREL=JMAJAC.VELCHE(/3) 132 : * 133 : IF (IREL.NE.IDIM.OR.IREF.NE.IDIM) THEN 134 : WRITE(IOIMP,*) 'Erreur dims JMAJAC' 135 : GOTO 9999 136 : ENDIF 137 : * 138 : ICOF=ICOF+1 139 : JMIJAC=LCOF.LISCHE(ICOF) 140 : IF (JMIJAC.EQ.0) THEN 141 : WRITE(IOIMP,*) 'Erreur JMIJAC=0' 142 : GOTO 9999 143 : ENDIF 144 : NLJI=JMIJAC.VELCHE(/6) 145 : NPJI=JMIJAC.VELCHE(/5) 146 : IREL2=JMIJAC.VELCHE(/4) 147 : IREF2=JMIJAC.VELCHE(/3) 148 : * 149 : IF (IREL2.NE.IREL.OR.IREF2.NE.IREF) THEN 150 : WRITE(IOIMP,*) 'Erreur dims JMIJAC' 151 : GOTO 9999 152 : ENDIF 153 : * 154 : ICOF=ICOF+1 155 : ICOF=ICOF+1 156 : JMAREG=LCOF.LISCHE(ICOF) 157 : NLJR=JMAREG.VELCHE(/6) 158 : NPJR=JMAREG.VELCHE(/5) 159 : I1 =JMAREG.VELCHE(/4) 160 : I2 =JMAREG.VELCHE(/3) 161 : IF ((NLJR.NE.1).OR.(NPJR.NE.1).OR.(I1.NE.IREF).OR.(I2.NE.IREF)) 162 : $ THEN 163 : WRITE(IOIMP,*) 'Erreur dims JMAREG' 164 : GOTO 9999 165 : ENDIF 166 : ICOF=ICOF+1 167 : JDIAMA=LCOF.LISCHE(ICOF) 168 : NLJD=JDIAMA.VELCHE(/6) 169 : NPJD=JDIAMA.VELCHE(/5) 170 : I1 =JDIAMA.VELCHE(/4) 171 : I2 =JDIAMA.VELCHE(/3) 172 : IF ((NLJD.NE.1).OR.(NPJD.NE.1).OR.(I1.NE.1).OR.(I2.NE.1)) 173 : $ THEN 174 : WRITE(IOIMP,*) 'Erreur dims JDIAMA' 175 : GOTO 9999 176 : ENDIF 177 : XDIAMA=JDIAMA.VELCHE(1,1,1,1,1,1) 178 : * 179 : * Objets temporaires 180 : * 181 : LDIM1=IREL 182 : LDIM2=1 183 : SEGINI,A 184 : LDIM1=IREF 185 : LDIM2=1 186 : SEGINI,JMA 187 : LDIM1=IREF 188 : LDIM2=1 189 : SEGINI,KJMA 190 : LDIM1=IREF 191 : LDIM2=IREL 192 : SEGINI,JM1 193 : LDIM1=IREF 194 : LDIM2=IREF 195 : SEGINI,K 196 : * 197 : * Calcul de la métrique des éléments réguliers 198 : * 198 : CALL MAMA(JMAREG.VELCHE,IREF,IREF, 199 : $ 'COPIE ',K.XMAT,IREF,IREF, 200 : $ IMPR,IRET) 201 : IF (IRET.NE.0) GOTO 9999 202 : * SEGPRT,H 203 : * 204 : DO ILFC=1,NLFC 205 : IF (NLRO.EQ.1) THEN 206 : ILRO=1 207 : ELSE 208 : ILRO=ILFC 209 : ENDIF 210 : IF (NLMU.EQ.1) THEN 211 : ILMU=1 212 : ELSE 213 : ILMU=ILFC 214 : ENDIF 215 : IF (NLJV.EQ.1) THEN 216 : ILJV=1 217 : ELSE 218 : ILJV=ILFC 219 : ENDIF 220 : IF (NLPEC.EQ.1) THEN 221 : ILPEC=1 222 : ELSE 223 : ILPEC=ILFC 224 : ENDIF 225 : C IF (NLJA.EQ.1) THEN 226 : C ILJA=1 227 : C ELSE 228 : C ILJA=ILFC 229 : C ENDIF 230 : IF (NLJI.EQ.1) THEN 231 : ILJI=1 232 : ELSE 233 : ILJI=ILFC 234 : ENDIF 235 : * 236 : DO IPFC=1,NPFC 237 : IF (NPRO.EQ.1) THEN 238 : IPRO=1 239 : ELSE 240 : IPRO=IPFC 241 : ENDIF 242 : IF (NPMU.EQ.1) THEN 243 : IPMU=1 244 : ELSE 245 : IPMU=IPFC 246 : ENDIF 247 : IF (NPJV.EQ.1) THEN 248 : IPJV=1 249 : ELSE 250 : IPJV=IPFC 251 : ENDIF 252 : IF (NPPEC.EQ.1) THEN 253 : IPPEC=1 254 : ELSE 255 : IPPEC=IPFC 256 : ENDIF 257 : C IF (NPJA.EQ.1) THEN 258 : C IPJA=1 259 : C ELSE 260 : C IPJA=IPFC 261 : C ENDIF 262 : IF (NPJI.EQ.1) THEN 263 : IPJI=1 264 : ELSE 265 : IPJI=IPFC 266 : ENDIF 267 : XRHO=JRHO.VELCHE(1,1,1,1,IPRO,ILRO) 268 : XMU=JMU.VELCHE(1,1,1,1,IPMU,ILMU) 269 : XPEC=JPEC.VELCHE(1,1,1,1,IPPEC,ILPEC) 270 : * 271 : * Copie de la vitesse 272 : * 273 : DO IIDIM=1,IDIM 274 : JVIT=MVIT.MVCOMP(IIDIM) 275 : A.XMAT(IIDIM,1)=JVIT.VELCHE(1,1,1,1,IPJV,ILJV) 276 : ENDDO 277 : * SEGPRT,A 278 : * 279 : * Copie de l'inverse du jacobien ref->reel 280 : * 280 : CALL MAMA(JMIJAC.VELCHE(1,1,1,1,IPJI,ILJI),IREF,IREL, 281 : $ 'COPIE ', 282 : $ JM1.XMAT,IREF,IREL, 283 : $ IMPR,IRET) 284 : IF (IRET.NE.0) GOTO 9999 285 : C SEGPRT,JM1 286 : * 287 : * Vecteur A dans le repère de référence 288 : * 288 : CALL MAMAMA(JM1.XMAT,IREF,IREL,A.XMAT,IREL,1, 289 : $ 'FOIS ',JMA.XMAT,IREF,1, 290 : $ IMPR,IRET) 291 : IF (IRET.NE.0) GOTO 9999 292 : * 293 : * Vecteur A dans le repère de l'élément régulier 294 : * 294 : CALL MAMAMA(K.XMAT,IREF,IREF,JMA.XMAT,IREF,1, 295 : $ 'FOIS ',KJMA.XMAT,IREF,1, 296 : $ IMPR,IRET) 297 : IF (IRET.NE.0) GOTO 9999 298 : * 299 : * Normes de A 300 : * 301 : XAREG2=XZERO 302 : DO IIREF=1,IREF 303 : XAREG2=XAREG2+(KJMA.XMAT(IIREF,1)**2) 304 : ENDDO 305 : XA2=XZERO 306 : DO IIDIM=1,IDIM 307 : XA2=XA2+(A.XMAT(IIDIM,1)**2) 308 : ENDDO 309 : * 310 : IF (XA2.LT.XPETI) THEN 311 : CONTRI=XZERO 312 : ELSE 313 : XV=SQRT(XA2) 314 : * Pondération par un diamètre externe ! 315 : * XH=SQRT(XA2/XAREG2) 316 : XH=SQRT(XA2/XAREG2)*XDIAMA 317 : * 318 : IF (IMETH.EQ.1) THEN 319 : * Diffusion artificielle UPWIND 320 : CONTRI=((XRHO*XV*XH)/XPEC) 321 : CONTRI=SQRT(CONTRI) 322 : IF (IDIMI.GT.0) THEN 323 : CONTRI=CONTRI* 324 : $ (A.XMAT(IDIMI,1)/XV) 325 : ENDIF 326 : ELSEIF (IMETH.EQ.2) THEN 327 : * Diffusion artificielle comme Alberto (equiv critical approx) 328 : * Doit-on décentrer ? 329 : LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU)) 330 : IF (LCRIT) THEN 331 : CONTRI=((XRHO*XV*XH)/XPEC)-XMU 332 : CONTRI=SQRT(CONTRI) 333 : IF (IDIMI.GT.0) THEN 334 : CONTRI=CONTRI* 335 : $ (A.XMAT(IDIMI,1)/XV) 336 : ENDIF 337 : ELSE 338 : CONTRI=XZERO 339 : ENDIF 340 : ELSEIF (IMETH.EQ.3) THEN 341 : * Diffusion artificielle (equiv doubly asymptotic) 342 : LCRIT=((XRHO*XV*XH).GT.(XPEC*XMU*3.D0)) 343 : IF (LCRIT) THEN 344 : CONTRI=((XRHO*XV*XH)/XPEC) 345 : ELSE 346 : CONTRI=(((XRHO*XV*XH)/XPEC)**2)/(3.D0*XMU) 347 : ENDIF 348 : CONTRI=SQRT(CONTRI) 349 : IF (IDIMI.GT.0) THEN 350 : CONTRI=CONTRI* 351 : $ (A.XMAT(IDIMI,1)/XV) 352 : ENDIF 353 : ELSE 354 : WRITE(IOIMP,*) 'IMETH=',IMETH,' unknown' 355 : GOTO 9999 356 : ENDIF 357 : ENDIF 358 : FC.VELCHE(1,1,1,1,IPFC,ILFC)= 359 : $ FC.VELCHE(1,1,1,1,IPFC,ILFC)+ 360 : $ CONTRI 361 : ENDDO 362 : ENDDO 363 : SEGSUP,K 364 : SEGSUP,JM1 365 : SEGSUP,KJMA 366 : SEGSUP,JMA 367 : SEGSUP,A 368 : SEGSUP,MVIT 369 : * 370 : * Normal termination 371 : * 372 : IRET=0 373 : RETURN 374 : * 375 : * Format handling 376 : * 377 : * 378 : * Error handling 379 : * 380 : 9999 CONTINUE 381 : IRET=1 382 : WRITE(IOIMP,*) 'An error was detected in subroutine ccgmus' 383 : RETURN 384 : * 385 : * End of subroutine CCGMUS 386 : * 387 : END 388 : 389 : 390 : 391 :
© Cast3M 2003 - All rights reserved.
Disclaimer