1 : C CADXFF SOURCE GOUNAND 05/12/21 21:15:03 5281 2 : SUBROUTINE CADXFF(JAC,NPJA,NLJA,IREL,IREF, 3 : $ JSDETG,IG,IGH,IGHIG, 4 : $ NOMLOI, 5 : $ FC,NPFC,NLFC, 6 : $ IMPR,IRET) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : IMPLICIT INTEGER (I-N) 9 : C*********************************************************************** 10 : C NOM : CADXFF 11 : C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : 12 : C Partie en Fortran 77 13 : C 14 : C LANGAGE : ESOPE 15 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 16 : C mél : gounand@semt2.smts.cea.fr 17 : C*********************************************************************** 18 : C APPELES : 19 : C APPELE PAR : 20 : C*********************************************************************** 21 : C ENTREES : 22 : C ENTREES/SORTIES : 23 : C SORTIES : - 24 : C TRAVAIL : 25 : C*********************************************************************** 26 : C VERSION : v1, 04/08/04, version initiale 27 : C HISTORIQUE : v1, 04/08/04, création 28 : C HISTORIQUE : 29 : C HISTORIQUE : 30 : C*********************************************************************** 31 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 32 : C en cas de modification de ce sous-programme afin de faciliter 33 : C la maintenance ! 34 : C*********************************************************************** 35 : -INC CCOPTIO 36 : REAL*8 JAC(IREL,IREF,NPJA,NLJA) 37 : REAL*8 JSDETG(NPJA,NLJA) 38 : REAL*8 IG(IREF,IREF,NPJA,NLJA) 39 : REAL*8 IGH(IREF,IREF,NPJA,NLJA) 40 : REAL*8 IGHIG(IREF,IREF,NPJA,NLJA) 41 : CHARACTER*8 NOMLOI 42 : REAL*8 FC(NPFC,NLFC) 43 : * 44 : INTEGER IMPR,IRET 45 : * 46 : * Executable statements 47 : * 48 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans cadxff' 49 : * 50 : DO ILFC=1,NLFC 51 : IF (NLJA.EQ.1) THEN 52 : ILJA=1 53 : ELSE 54 : ILJA=ILFC 55 : ENDIF 56 : DO IPFC=1,NPFC 57 : IF (NPJA.EQ.1) THEN 58 : IPJA=1 59 : ELSE 60 : IPJA=IPFC 61 : ENDIF 62 : SDETG=JSDETG(IPJA,ILJA) 63 : * 64 : 64 : CALL MARE(IGH(1,1,IPJA,ILJA),IREF,IREF,'TRACE ', 65 : $ TRIGH,IMPR,IRET) 66 : IF (IRET.NE.0) GOTO 9999 67 : * 68 : * Calcul de la fonctionnelle 69 : * 70 : IF (NOMLOI.EQ.'ADCF ') THEN 71 : CONTRI=0.5D0*TRIGH*SDETG 72 : * 73 : * Calcul du résidu 74 : * 75 : ELSEIF (NOMLOI(1:4).EQ.'ADCR') THEN 75 : CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 76 : IF (IRET.NE.0) GOTO 9999 76 : CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 77 : IF (IRET.NE.0) GOTO 9999 78 : TR1=0.D0 79 : DO I=1,IREF 80 : DO J=1,IREF 81 : TR1=TR1+ 82 : $ (IG(I,J,IPJA,ILJA)*JAC(IDIM1,J,IPJA,ILJA) 83 : $ *IGH(IDIM2,I,IPJA,ILJA)) 84 : ENDDO 85 : ENDDO 86 : TR2=0.D0 87 : DO I=1,IREF 88 : DO L=1,IREF 89 : TR2=TR2+ 90 : $ (IG(I,IDIM2,IPJA,ILJA)*JAC(IDIM1,L,IPJA,ILJA) 91 : $ *IGH(L,I,IPJA,ILJA)) 92 : ENDDO 93 : ENDDO 94 : CONTR1=-0.5D0*(TR1+TR2)*SDETG 95 : TR3=0.D0 96 : DO J=1,IREF 97 : TR3=TR3+ 98 : $ (IG(IDIM2,J,IPJA,ILJA)*JAC(IDIM1,J,IPJA,ILJA)) 99 : ENDDO 100 : TR4=0.D0 101 : DO I=1,IREF 102 : TR4=TR4+ 103 : $ (IG(I,IDIM2,IPJA,ILJA)*JAC(IDIM1,I,IPJA,ILJA)) 104 : ENDDO 105 : CONTR2=0.25D0*(TR3+TR4)*TRIGH*SDETG 106 : CONTRI=CONTR1+CONTR2 107 : * CONTRI=CONTR1 108 : C WRITE(IOIMP,*) 'ADCR' 109 : C WRITE(IOIMP,*) 'TR1=',TR1 110 : C WRITE(IOIMP,*) 'TR2=',TR2 111 : C WRITE(IOIMP,*) 'TR3=',TR3 112 : C WRITE(IOIMP,*) 'TR4=',TR4 113 : * 114 : * Calcul du jacobien 115 : * 116 : ELSEIF (NOMLOI(1:4).EQ.'ADCJ') THEN 116 : CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 117 : IF (IRET.NE.0) GOTO 9999 117 : CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 118 : IF (IRET.NE.0) GOTO 9999 118 : CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) 119 : IF (IRET.NE.0) GOTO 9999 119 : CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) 120 : IF (IRET.NE.0) GOTO 9999 121 : TR0=0.D0 122 : DO K=1,IREF 123 : TR0=TR0+ 124 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,IDIM2,IPJA,ILJA)) 125 : ENDDO 126 : * 127 : TR1=0.D0 128 : DO I=1,IREF 129 : DO J=1,IREF 130 : TR1=TR1+ 131 : $ (IG(I,J,IPJA,ILJA)*JAC(IDIM3,J,IPJA,ILJA) 132 : $ *IGH(IDIM4,I,IPJA,ILJA)) 133 : ENDDO 134 : ENDDO 135 : TR2=0.D0 136 : DO I=1,IREF 137 : DO M=1,IREF 138 : TR2=TR2+ 139 : $ (IG(I,IDIM4,IPJA,ILJA)*JAC(IDIM3,M,IPJA,ILJA) 140 : $ *IGH(M,I,IPJA,ILJA)) 141 : ENDDO 142 : ENDDO 143 : * 144 : * 145 : TR3=0.D0 146 : DO K=1,IREF 147 : DO L=1,IREF 148 : TR3=TR3+ 149 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,L,IPJA,ILJA) 150 : $ *JAC(IDIM3,L,IPJA,ILJA)*IGHIG(IDIM4,IDIM2,IPJA 151 : $ ,ILJA)) 152 : ENDDO 153 : ENDDO 154 : TR3=TR3*(-1.D0) 155 : * 156 : TR4=0.D0 157 : DO K=1,IREF 158 : DO IO=1,IREF 159 : TR4=TR4+ 160 : $ (JAC(IDIM1,K,IPJA,ILJA)*IGHIG(K,IDIM4,IPJA 161 : $ ,ILJA)*JAC(IDIM3,IO,IPJA,ILJA)*IG(IO,IDIM2 162 : $ ,IPJA,ILJA)) 163 : ENDDO 164 : ENDDO 165 : TR4=TR4*(-1.D0) 166 : * 167 : TR5=0.D0 168 : IF (IDIM1.EQ.IDIM3) THEN 169 : TR5=TR5+ 170 : $ (IGHIG(IDIM4,IDIM2,IPJA,ILJA)) 171 : ENDIF 172 : * 173 : TR6=0.D0 174 : DO K=1,IREF 175 : DO M=1,IREF 176 : TR6=TR6+ 177 : $ (JAC(IDIM1,K,IPJA,ILJA)*IGHIG(K,M,IPJA,ILJA) 178 : $ *JAC(IDIM3,M,IPJA,ILJA)*IG(IDIM4,IDIM2,IPJA 179 : $ ,ILJA)) 180 : ENDDO 181 : ENDDO 182 : TR6=TR6*(-1.D0) 183 : * 184 : TR7=0.D0 185 : DO K=1,IREF 186 : DO N=1,IREF 187 : TR7=TR7+ 188 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,IDIM4,IPJA,ILJA) 189 : $ *JAC(IDIM3,N,IPJA,ILJA)*IGHIG(N,IDIM2,IPJA 190 : $ ,ILJA)) 191 : ENDDO 192 : ENDDO 193 : TR7=TR7*(-1.D0) 194 : CONTR1=(TR1+TR2)*TR0 195 : CONTR2=(TR3+TR4+TR5+TR6+TR7)*2.D0 196 : CONT11=(CONTR1+CONTR2)*SDETG 197 : **** 198 : **** 199 : TR10=0.D0 200 : DO K=1,IREF 201 : TR10=TR10+ 202 : $ (JAC(IDIM1,K,IPJA,ILJA)*IGHIG(K,IDIM2,IPJA,ILJA)) 203 : ENDDO 204 : TR10=TR10*(-2.D0) 205 : * 206 : TR0=0.D0 207 : DO K=1,IREF 208 : TR0=TR0+ 209 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,IDIM2,IPJA,ILJA)) 210 : ENDDO 211 : * 212 : TR11=0.D0 213 : DO J=1,IREF 214 : TR11=TR11+ 215 : $ (IG(IDIM4,J,IPJA,ILJA)*JAC(IDIM3,J,IPJA,ILJA)) 216 : ENDDO 217 : TR12=0.D0 218 : DO I=1,IREF 219 : TR12=TR12+ 220 : $ (IG(I,IDIM4,IPJA,ILJA)*JAC(IDIM3,I,IPJA,ILJA)) 221 : ENDDO 222 : 223 : * 224 : * 225 : TR13=0.D0 226 : DO K=1,IREF 227 : DO L=1,IREF 228 : TR13=TR13+ 229 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,L,IPJA,ILJA) 230 : $ *JAC(IDIM3,L,IPJA,ILJA)*IG(IDIM4,IDIM2,IPJA 231 : $ ,ILJA)) 232 : ENDDO 233 : ENDDO 234 : TR13=TR13*(-1.D0) 235 : TR14=0.D0 236 : DO K=1,IREF 237 : DO N=1,IREF 238 : TR14=TR14+ 239 : $ (JAC(IDIM1,K,IPJA,ILJA)*IG(K,IDIM4,IPJA,ILJA) 240 : $ *JAC(IDIM3,N,IPJA,ILJA)*IG(N,IDIM2,IPJA,ILJA)) 241 : ENDDO 242 : ENDDO 243 : TR14=TR14*(-1.D0) 244 : * 245 : TR15=0.D0 246 : IF (IDIM1.EQ.IDIM3) THEN 247 : TR15=TR15+IG(IDIM4,IDIM2,IPJA,ILJA) 248 : ENDIF 249 : * 250 : CONTR3=(TR11+TR12)*(TR10+(TR0*TRIGH))*SDETG 251 : CONTR4=(TR13+TR14+TR15)*2.D0*TRIGH*SDETG 252 : CONT12=(CONTR3+CONTR4) 253 : ***** 254 : ***** Finalement ! 255 : ***** 256 : * WRITE(IOIMP,*) 'CONT11=',CONT11 257 : * WRITE(IOIMP,*) 'CONT12=',CONT12 258 : CONTRI=(CONT11*(-0.5D0)) 259 : $ +(CONT12*0.25D0) 260 : * CONTRI=(CONT11*(-0.5D0)) 261 : C* 262 : C* Calcul du résidu (test) 263 : C* 264 : C ELSEIF (NOMLOI(1:4).EQ.'ADCS') THEN 264 : C CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 265 : C IF (IRET.NE.0) GOTO 9999 265 : C CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 266 : C IF (IRET.NE.0) GOTO 9999 267 : CC TR1=JAC(IDIM1,IDIM2) 268 : C TR1=0.D0 269 : C DO I=1,IREF 270 : C DO J=1,IREF 271 : C TR1=TR1+ 272 : C $ (IG(I,J)*JAC(IDIM1,J)*IGH(IDIM2,I)) 273 : C* $ (IG(I,J)*JAC(IDIM1,J)*IG(IDIM2,I)) 274 : C ENDDO 275 : C ENDDO 276 : C TR2=0.D0 277 : C DO I=1,IREF 278 : C DO M=1,IREF 279 : C TR2=TR2+ 280 : C $ (IG(I,IDIM2)*JAC(IDIM1,M)*IGH(M,I)) 281 : C ENDDO 282 : C ENDDO 283 : C CONTRI=(TR1+TR2)*SDETG 284 : C* CONTRI=(TR1+TR2) 285 : C* CONTRI=TR1 286 : C* 287 : C* Calcul du jacobien (test) 288 : C* 289 : C ELSEIF (NOMLOI(1:4).EQ.'ADCK') THEN 289 : C CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 290 : C IF (IRET.NE.0) GOTO 9999 290 : C CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 291 : C IF (IRET.NE.0) GOTO 9999 291 : C CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) 292 : C IF (IRET.NE.0) GOTO 9999 292 : C CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) 293 : C IF (IRET.NE.0) GOTO 9999 294 : CC IF ((IDIM1.EQ.IDIM3).AND.(IDIM2.EQ.IDIM4)) THEN 295 : CC CONTRI=1.D0 296 : CC ELSE 297 : CC CONTRI=0.D0 298 : CC ENDIF 299 : C TR0=0.D0 300 : C DO K=1,IREF 301 : C TR0=TR0+ 302 : C $ (JAC(IDIM1,K)*IG(K,IDIM2)) 303 : C ENDDO 304 : C* 305 : C TR1=0.D0 306 : C DO I=1,IREF 307 : C DO J=1,IREF 308 : C TR1=TR1+ 309 : C $ (IG(I,J)*JAC(IDIM3,J)*IGH(IDIM4,I)) 310 : C ENDDO 311 : C ENDDO 312 : C TR2=0.D0 313 : C DO I=1,IREF 314 : C DO M=1,IREF 315 : C TR2=TR2+ 316 : C $ (IG(I,IDIM4)*JAC(IDIM3,M)*IGH(M,I)) 317 : C ENDDO 318 : C ENDDO 319 : C* 320 : C* 321 : C TR3=0.D0 322 : C DO K=1,IREF 323 : C DO L=1,IREF 324 : C DO IO=1,IREF 325 : C TR3=TR3+ 326 : C $ (JAC(IDIM1,K)*IG(K,L)*JAC(IDIM3,L) 327 : C $ *IG(IDIM4,IO)*HIG(IO,IDIM2)) 328 : C* $ (JAC(IDIM1,K)*IG(K,L)*JAC(IDIM3,L) 329 : C* $ *IG(IDIM4,IO)*IG(IO,IDIM2)) 330 : C ENDDO 331 : C ENDDO 332 : C ENDDO 333 : C TR3=TR3*(-1.D0) 334 : C* 335 : C TR4=0.D0 336 : C DO K=1,IREF 337 : C DO L=1,IREF 338 : C DO IO=1,IREF 339 : C TR4=TR4+ 340 : C $ (JAC(IDIM1,K)*IGH(K,L)*IG(L,IDIM4) 341 : C $ *JAC(IDIM3,IO)*IG(IO,IDIM2)) 342 : C* $ (JAC(IDIM1,K)*IG(K,L)*IG(L,IDIM4) 343 : C* $ *JAC(IDIM3,IO)*IG(IO,IDIM2)) 344 : C ENDDO 345 : C ENDDO 346 : C ENDDO 347 : C TR4=TR4*(-1.D0) 348 : C* 349 : C TR5=0.D0 350 : C IF (IDIM1.EQ.IDIM3) THEN 351 : C DO L=1,IREF 352 : C TR5=TR5+ 353 : C $ (IGH(IDIM4,L)*IG(L,IDIM2)) 354 : C* $ (IG(IDIM4,L)*IG(L,IDIM2)) 355 : C ENDDO 356 : C ENDIF 357 : C* 358 : C TR6=0.D0 359 : C DO K=1,IREF 360 : C DO L=1,IREF 361 : C DO M=1,IREF 362 : C TR6=TR6+ 363 : C $ (JAC(IDIM1,K)*IGH(K,L)*IG(L,M) 364 : C $ *JAC(IDIM3,M)*IG(IDIM4,IDIM2)) 365 : C ENDDO 366 : C ENDDO 367 : C ENDDO 368 : C TR6=TR6*(-1.D0) 369 : C* 370 : C TR7=0.D0 371 : C DO K=1,IREF 372 : C DO N=1,IREF 373 : C DO IO=1,IREF 374 : C TR7=TR7+ 375 : C $ (JAC(IDIM1,K)*IG(K,IDIM4)*JAC(IDIM3,N) 376 : C $ *IG(N,IO)*HIG(IO,IDIM2)) 377 : C ENDDO 378 : C ENDDO 379 : C ENDDO 380 : C TR7=TR7*(-1.D0) 381 : C* 382 : C CONTR1=(TR1+TR2)*TR0 383 : C CONTR2=(TR3+TR4+TR5+TR6+TR7)*2.D0 384 : C CONTRI=(CONTR1+CONTR2)*SDETG 385 : C* CONTRI=(TR3+TR4+TR5+TR6+TR7)+(TR3B+TR4B+TR5B+TR6B+TR7B) 386 : C* CONTRI=(TR3+TR4+TR5+TR6+TR7) 387 : C* 388 : C* WRITE(IOIMP,*) 'CONTRI=',CONTRI 389 : C* 390 : C* Calcul du résidu (test) 391 : C* 392 : C ELSEIF (NOMLOI(1:4).EQ.'ADCS') THEN 392 : C CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 393 : C IF (IRET.NE.0) GOTO 9999 393 : C CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 394 : C IF (IRET.NE.0) GOTO 9999 395 : C TR11=0.D0 396 : C DO J=1,IREF 397 : C TR11=TR11+ 398 : C $ (IG(IDIM2,J)*JAC(IDIM1,J)) 399 : C ENDDO 400 : C TR12=0.D0 401 : C DO I=1,IREF 402 : C TR12=TR12+ 403 : C $ (IG(I,IDIM2)*JAC(IDIM1,I)) 404 : C ENDDO 405 : C CONTRI=(TR11+TR12)*TRIGH*SDETG 406 : C* 407 : C* Calcul du jacobien (test) 408 : C* 409 : C ELSEIF (NOMLOI(1:4).EQ.'ADCK') THEN 409 : C CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 410 : C IF (IRET.NE.0) GOTO 9999 410 : C CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 411 : C IF (IRET.NE.0) GOTO 9999 411 : C CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) 412 : C IF (IRET.NE.0) GOTO 9999 412 : C CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) 413 : C IF (IRET.NE.0) GOTO 9999 414 : C* 415 : C TR10=0.D0 416 : C DO K=1,IREF 417 : C DO L=1,IREF 418 : C TR10=TR10+ 419 : C $ (JAC(IDIM1,K)*IG(K,L)*HIG(L,IDIM2)) 420 : C ENDDO 421 : C ENDDO 422 : C TR10=TR10*(-2.D0) 423 : C* 424 : C TR0=0.D0 425 : C DO K=1,IREF 426 : C TR0=TR0+ 427 : C $ (JAC(IDIM1,K)*IG(K,IDIM2)) 428 : C ENDDO 429 : C* 430 : C TR11=0.D0 431 : C DO J=1,IREF 432 : C TR11=TR11+ 433 : C $ (IG(IDIM4,J)*JAC(IDIM3,J)) 434 : C ENDDO 435 : C TR12=0.D0 436 : C DO I=1,IREF 437 : C TR12=TR12+ 438 : C $ (IG(I,IDIM4)*JAC(IDIM3,I)) 439 : C ENDDO 440 : C 441 : C* 442 : C* 443 : C TR13=0.D0 444 : C DO K=1,IREF 445 : C DO L=1,IREF 446 : C TR13=TR13+ 447 : C $ (JAC(IDIM1,K)*IG(K,L)*JAC(IDIM3,L) 448 : C $ *IG(IDIM4,IDIM2)) 449 : C ENDDO 450 : C ENDDO 451 : C TR13=TR13*(-1.D0) 452 : C TR14=0.D0 453 : C DO K=1,IREF 454 : C DO N=1,IREF 455 : C TR14=TR14+ 456 : C $ (JAC(IDIM1,K)*IG(K,IDIM4) 457 : C $ *JAC(IDIM3,N)*IG(N,IDIM2)) 458 : C ENDDO 459 : C ENDDO 460 : C TR14=TR14*(-1.D0) 461 : C* 462 : C TR15=0.D0 463 : C IF (IDIM1.EQ.IDIM3) THEN 464 : C TR15=TR15+IG(IDIM4,IDIM2) 465 : C ENDIF 466 : C* 467 : CC CONTRI=(TR13+TR14+TR15)*2.D0 468 : C CONTR3=(TR11+TR12)*(TR10+(TR0*TRIGH))*SDETG 469 : C CONTR4=(TR13+TR14+TR15)*2.D0*TRIGH*SDETG 470 : CC CONTR1=(TR1+TR2)*TR0 471 : CC CONTR2=(TR3+TR4+TR5+TR6+TR7)*2.D0 472 : CC CONTRI=(CONTR1+CONTR2)*SDETG 473 : C CONTRI=(CONTR3+CONTR4) 474 : C* 475 : C WRITE(IOIMP,*) 'CONTRI=',CONTRI 476 : ELSE 477 : WRITE(IOIMP,*) 'Erreur grave' 478 : GOTO 9999 479 : ENDIF 480 : FC(IPFC,ILFC)= 481 : $ FC(IPFC,ILFC)+ 482 : $ CONTRI 483 : ENDDO 484 : ENDDO 485 : * 486 : * Normal termination 487 : * 488 : IRET=0 489 : RETURN 490 : * 491 : * Format handling 492 : * 493 : * 494 : * Error handling 495 : * 496 : 9999 CONTINUE 497 : IRET=1 498 : WRITE(IOIMP,*) 'An error was detected in subroutine cadxff' 499 : RETURN 500 : * 501 : * End of subroutine CADXFF 502 : * 503 : END 504 : 505 :
© Cast3M 2003 - All rights reserved.
Disclaimer