1 : C CCGTSU SOURCE GOUNAND 11/04/29 21:15:21 6947 2 : SUBROUTINE CCGTSU(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 : CCGTSU 9 : C DESCRIPTION : Calcul de la loi de comportement aux points de Gauss : 10 : C Tension de surface 11 : C 12 : C 13 : C LANGAGE : ESOPE 14 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 15 : C mél : gounand@semt2.smts.cea.fr 16 : C*********************************************************************** 17 : C APPELES : 18 : C APPELE PAR : 19 : C*********************************************************************** 20 : C ENTREES : 21 : C ENTREES/SORTIES : 22 : C SORTIES : - 23 : C TRAVAIL : 24 : C*********************************************************************** 25 : C VERSION : v1, 17/01/07, version initiale 26 : C HISTORIQUE : v1, 17/01/07, création 27 : C HISTORIQUE : 28 : C HISTORIQUE : 29 : C*********************************************************************** 30 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 31 : C en cas de modification de ce sous-programme afin de faciliter 32 : C la maintenance ! 33 : C*********************************************************************** 34 : -INC CCOPTIO 35 : -INC CCREEL 36 : CBEGININCLUDE SMCHAEL 37 : SEGMENT MCHAEL 38 : POINTEUR IMACHE(N1).MELEME 39 : POINTEUR ICHEVA(N1).MCHEVA 40 : ENDSEGMENT 41 : SEGMENT MCHEVA 42 : REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM) 43 : ENDSEGMENT 44 : SEGMENT LCHEVA 45 : POINTEUR LISCHE(NBCHE).MCHEVA 46 : ENDSEGMENT 47 : CENDINCLUDE SMCHAEL 48 : INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM,N1 49 : POINTEUR FC.MCHEVA 50 : POINTEUR LCOF.LCHEVA 51 : POINTEUR JMAJAC.MCHEVA 52 : POINTEUR JMIJAC.MCHEVA 53 : POINTEUR JDTJAC.MCHEVA 54 : POINTEUR JMAREG.MCHEVA 55 : POINTEUR JPC.MCHEVA 56 : POINTEUR JTS.MCHEVA 57 : CHARACTER*8 NOMLOI 58 : INTEGER ICOF 59 : * 60 : -INC TMXMAT 61 : POINTEUR JAC.MXMAT,JT.MXMAT 62 : POINTEUR G.MXMAT 63 : POINTEUR IG.MXMAT 64 : POINTEUR JPT.MXMAT 65 : POINTEUR JPTJT.MXMAT 66 : * 67 : LOGICAL LBID 68 : INTEGER LAXSP 69 : REAL*8 DEUPI,QUATPI,XR 70 : PARAMETER (DEUPI=2.D0*XPI,QUATPI=4.D0*XPI) 71 : * 72 : INTEGER IMPR,IRET 73 : * 74 : * Executable statements 75 : * 76 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ccgtsu' 77 : C IF (.NOT.(IDIM.EQ.1)) THEN 78 : C WRITE(IOIMP,*) 'IDIM=',IDIM,' ?' 79 : C GOTO 9999 80 : C ENDIF 81 : NLFC=FC.VELCHE(/6) 82 : NPFC=FC.VELCHE(/5) 83 : ICOF=0 84 : * 85 : * Récupération du coefficent de tension de surface 86 : * 87 : ICOF=ICOF+1 88 : JTS =LCOF.LISCHE(ICOF) 89 : NLTS=JTS.VELCHE(/6) 90 : NPTS=JTS.VELCHE(/5) 91 : * 92 : ICOF=ICOF+1 93 : JMAJAC=LCOF.LISCHE(ICOF) 94 : NLJA=JMAJAC.VELCHE(/6) 95 : NPJA=JMAJAC.VELCHE(/5) 96 : IREF=JMAJAC.VELCHE(/4) 97 : IREL=JMAJAC.VELCHE(/3) 98 : * 99 : IF (IREL.NE.IDIM) THEN 100 : WRITE(IOIMP,*) 'Erreur dims JMAJAC' 101 : GOTO 9999 102 : ENDIF 103 : * 104 : ICOF=ICOF+1 105 : ICOF=ICOF+1 106 : ICOF=ICOF+1 107 : * 108 : *! WRITE(IOIMP,*) 'IFOMOD=',IFOMOD 109 : IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN 110 : *! WRITE(IOIMP,*) 'MODE AXISYMETRIQUE !!!!!!' 111 : LAXSP=1 112 : ELSEIF (IFOMOD.EQ.5) THEN 113 : *! WRITE(IOIMP,*) 'MODE SPHERIQUE !!!!!!' 114 : LAXSP=2 115 : ELSE 116 : *! WRITE(IOIMP,*) 'MODE PLAN !!!!!!' 117 : LAXSP=0 118 : ENDIF 119 : *! WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 120 : *! WRITE(IOIMP,*) 'LAXSP=',LAXSP 121 : * 122 : IF (LAXSP.GT.0) THEN 123 : ICOF=ICOF+1 124 : ICOF=ICOF+1 125 : JPC=LCOF.LISCHE(ICOF) 126 : NLPC=JPC.VELCHE(/6) 127 : NPPC=JPC.VELCHE(/5) 128 : IF (((NLPC.NE.1).AND.(NLPC.NE.NLFC)).OR. 129 : $ ((NPPC.NE.1).AND.(NPPC.NE.NPFC))) THEN 130 : WRITE(IOIMP,*) 'Erreur dims JPC' 131 : GOTO 9999 132 : ENDIF 133 : ELSE 134 : JPC=0 135 : NLPC=0 136 : NPPC=0 137 : ENDIF 138 : * 139 : * Objets temporaires 140 : * 141 : LDIM1=IREL 142 : LDIM2=IREF 143 : SEGINI,JAC 144 : SEGINI,JPT 145 : LDIM1=IREF 146 : LDIM2=IREL 147 : SEGINI,JT 148 : LDIM1=IREF 149 : LDIM2=IREF 150 : SEGINI,G 151 : SEGINI,IG 152 : LDIM1=IREL 153 : LDIM2=IREL 154 : SEGINI,JPTJT 155 : * 156 : DO ILFC=1,NLFC 157 : IF (NLTS.EQ.1) THEN 158 : ILTS=1 159 : ELSE 160 : ILTS=ILFC 161 : ENDIF 162 : IF (NLJA.EQ.1) THEN 163 : ILJA=1 164 : ELSE 165 : ILJA=ILFC 166 : ENDIF 167 : IF (NLPC.EQ.1) THEN 168 : ILPC=1 169 : ELSE 170 : ILPC=ILFC 171 : ENDIF 172 : 173 : 174 : DO IPFC=1,NPFC 175 : IF (NPTS.EQ.1) THEN 176 : IPTS=1 177 : ELSE 178 : IPTS=IPFC 179 : ENDIF 180 : IF (NPJA.EQ.1) THEN 181 : IPJA=1 182 : ELSE 183 : IPJA=IPFC 184 : ENDIF 185 : IF (NPPC.EQ.1) THEN 186 : IPPC=1 187 : ELSE 188 : IPPC=IPFC 189 : ENDIF 190 : * 191 : * Valeur du coefficient 192 : * 193 : IF (JTS.NE.0) THEN 194 : XTS=JTS.VELCHE(1,1,1,1,IPTS,ILTS) 195 : ENDIF 196 : * 197 : * Première coordonnée 198 : * 199 : IF (JPC.NE.0) THEN 200 : XR=JPC.VELCHE(1,1,1,1,IPPC,ILPC) 201 : ENDIF 202 : * 203 : * Copie du jacobien 204 : * 204 : CALL MAMA(JMAJAC.VELCHE(1,1,1,1,IPJA,ILJA),IREL,IREF, 205 : $ 'COPIE ', 206 : $ JAC.XMAT,IREL,IREF, 207 : $ IMPR,IRET) 208 : IF (IRET.NE.0) GOTO 9999 209 : C SEGPRT,JAC 210 : * 211 : * Calcul de la métrique G 212 : * 213 : * Calcul de Jt 213 : CALL MAMA(JAC.XMAT,IREL,IREF, 214 : $ 'TRANSPOS',JT.XMAT,IREF,IREL, 215 : $ IMPR,IRET) 216 : IF (IRET.NE.0) GOTO 9999 217 : * Calcul de G=JtJ 217 : CALL MAMAMA(JT.XMAT,IREF,IREL,JAC.XMAT,IREL,IREF, 218 : $ 'FOIS ',G.XMAT,IREF,IREF, 219 : $ IMPR,IRET) 220 : IF (IRET.NE.0) GOTO 9999 221 : * 222 : * Calcul de l'inverse, du déterminant et trace de l'inverse 223 : * 224 : LBID=.FALSE. 224 : CALL GEOLI2(IREF,1,1,G.XMAT,IG.XMAT,DETG,LBID,IMPR,IRET) 225 : IF (IRET.NE.0) GOTO 9999 226 : SDETG=SQRT(DETG) 227 : * Calcul de JPT=J (JtJ)-1 227 : CALL MAMAMA(JAC.XMAT,IREL,IREF,IG.XMAT,IREF,IREF, 228 : $ 'FOIS ',JPT.XMAT,IREL,IREF, 229 : $ IMPR,IRET) 230 : IF (IRET.NE.0) GOTO 9999 231 : * Calcul de JPT=J (JtJ)-1 231 : CALL MAMAMA(JPT.XMAT,IREL,IREF,JT.XMAT,IREF,IREL, 232 : $ 'FOIS ',JPTJT.XMAT,IREL,IREL, 233 : $ IMPR,IRET) 234 : IF (IRET.NE.0) GOTO 9999 235 : * 236 : * Calcul de la fonctionnelle 237 : * 238 : IF (NOMLOI.EQ.'TSUF ') THEN 239 : CONTRI=XTS*SDETG 240 : IF (LAXSP.EQ.1) THEN 241 : CONTRI=CONTRI*2.D0*XPI*XR 242 : ELSEIF (LAXSP.EQ.2) THEN 243 : CONTRI=CONTRI*4.D0*XPI*XR*XR 244 : ENDIF 245 : * 246 : * Calcul du résidu 247 : * 248 : ELSEIF (NOMLOI(1:4).EQ.'TSUR') THEN 248 : CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 249 : IF (IRET.NE.0) GOTO 9999 249 : CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 250 : IF (IRET.NE.0) GOTO 9999 251 : IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN 252 : IF (LAXSP.EQ.1) THEN 253 : CONTRI=XTS*2.D0*XPI*SDETG 254 : ELSEIF (LAXSP.EQ.2) THEN 255 : CONTRI=XTS*8.D0*XPI*XR*SDETG 256 : ELSE 257 : WRITE(IOIMP,*) 'LAXSP=',LAXSP 258 : GOTO 9999 259 : ENDIF 260 : ELSE 261 : CONTRI=XTS*SDETG*JPT.XMAT(IDIM1,IDIM2) 262 : IF (LAXSP.EQ.1) THEN 263 : CONTRI=CONTRI*DEUPI*XR 264 : ELSEIF (LAXSP.EQ.2) THEN 265 : CONTRI=CONTRI*QUATPI*XR*XR 266 : ENDIF 267 : ENDIF 268 : * 269 : * Calcul du jacobien 270 : * 271 : ELSEIF (NOMLOI(1:4).EQ.'TSUJ') THEN 271 : CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 272 : IF (IRET.NE.0) GOTO 9999 272 : CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 273 : IF (IRET.NE.0) GOTO 9999 273 : CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) 274 : IF (IRET.NE.0) GOTO 9999 274 : CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) 275 : IF (IRET.NE.0) GOTO 9999 276 : IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0) 277 : $ .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN 278 : IF (LAXSP.EQ.2) THEN 279 : CONTRI=XTS*8.D0*XPI*SDETG 280 : ELSE 281 : WRITE(IOIMP,*) 'LAXSP=',LAXSP 282 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 283 : GOTO 9999 284 : ENDIF 285 : ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN 286 : IF (LAXSP.GE.1) THEN 287 : CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM3,IDIM4) 288 : IF (LAXSP.EQ.2) THEN 289 : CONTRI=CONTRI*4.D0*XR 290 : ENDIF 291 : ELSE 292 : WRITE(IOIMP,*) 'LAXSP=',LAXSP 293 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 294 : GOTO 9999 295 : ENDIF 296 : ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN 297 : IF (LAXSP.GE.1) THEN 298 : CONTRI=XTS*DEUPI*SDETG*JPT.XMAT(IDIM1,IDIM2) 299 : IF (LAXSP.EQ.2) THEN 300 : CONTRI=CONTRI*4.D0*XR 301 : ENDIF 302 : ELSE 303 : WRITE(IOIMP,*) 'LAXSP=',LAXSP 304 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 305 : GOTO 9999 306 : ENDIF 307 : ELSE 308 : CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4) 309 : IF (IDIM1.EQ.IDIM3) THEN 310 : CONTRI=CONTRI+IG.XMAT(IDIM2,IDIM4) 311 : ENDIF 312 : CONTRI=CONTRI- 313 : $ (JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4)) 314 : CONTRI=CONTRI- 315 : $ (JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4)) 316 : CONTRI=CONTRI*XTS*SDETG 317 : IF (LAXSP.EQ.1) THEN 318 : CONTRI=CONTRI*DEUPI*XR 319 : ELSEIF (LAXSP.EQ.2) THEN 320 : CONTRI=CONTRI*QUATPI*XR*XR 321 : ENDIF 322 : ENDIF 323 : * 324 : * Calcul d'une partie du jacobien 325 : * 326 : ELSEIF (NOMLOI(1:3).EQ.'TSU') THEN 326 : CALL CH2INT(NOMLOI(4:4),IDIM0,IMPR,IRET) 327 : IF (IRET.NE.0) GOTO 9999 327 : CALL CH2INT(NOMLOI(5:5),IDIM1,IMPR,IRET) 328 : IF (IRET.NE.0) GOTO 9999 328 : CALL CH2INT(NOMLOI(6:6),IDIM2,IMPR,IRET) 329 : IF (IRET.NE.0) GOTO 9999 329 : CALL CH2INT(NOMLOI(7:7),IDIM3,IMPR,IRET) 330 : IF (IRET.NE.0) GOTO 9999 330 : CALL CH2INT(NOMLOI(8:8),IDIM4,IMPR,IRET) 331 : IF (IRET.NE.0) GOTO 9999 332 : IF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0) 333 : $ .AND.(IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN 334 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 335 : GOTO 9999 336 : ELSEIF ((IDIM1.EQ.1).AND.(IDIM2.EQ.0)) THEN 337 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 338 : GOTO 9999 339 : ELSEIF ((IDIM3.EQ.1).AND.(IDIM4.EQ.0)) THEN 340 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 341 : GOTO 9999 342 : ELSE 343 : IF (IDIM0.EQ.1) THEN 344 : CONTRI=JPT.XMAT(IDIM1,IDIM2)*JPT.XMAT(IDIM3,IDIM4) 345 : ELSEIF (IDIM0.EQ.2) THEN 346 : IF (IDIM1.EQ.IDIM3) THEN 347 : CONTRI=IG.XMAT(IDIM2,IDIM4) 348 : ELSE 349 : CONTRI=0.D0 350 : ENDIF 351 : ELSEIF (IDIM0.EQ.3) THEN 352 : CONTRI= 353 : $ -(JPTJT.XMAT(IDIM1,IDIM3)*IG.XMAT(IDIM2,IDIM4) 354 : $ ) 355 : ELSEIF (IDIM0.EQ.4) THEN 356 : CONTRI=-(JPT.XMAT(IDIM3,IDIM2)*JPT.XMAT(IDIM1,IDIM4 357 : $ )) 358 : ELSE 359 : WRITE(IOIMP,*) 'NOMLOI=',NOMLOI 360 : GOTO 9999 361 : ENDIF 362 : CONTRI=CONTRI*XTS*SDETG 363 : IF (LAXSP.EQ.1) THEN 364 : CONTRI=CONTRI*DEUPI*XR 365 : ELSEIF (LAXSP.EQ.2) THEN 366 : CONTRI=CONTRI*QUATPI*XR*XR 367 : ENDIF 368 : ENDIF 369 : *! WRITE(IOIMP,*) 'CONTRI=',CONTRI 370 : ELSE 371 : WRITE(IOIMP,*) 'Erreur grave' 372 : GOTO 9999 373 : ENDIF 374 : FC.VELCHE(1,1,1,1,IPFC,ILFC)= 375 : $ FC.VELCHE(1,1,1,1,IPFC,ILFC)+ 376 : $ CONTRI 377 : ENDDO 378 : ENDDO 379 : SEGSUP,JPTJT 380 : SEGSUP,JPT 381 : SEGSUP,IG 382 : SEGSUP,G 383 : SEGSUP,JT 384 : SEGSUP,JAC 385 : * 386 : * Normal termination 387 : * 388 : IRET=0 389 : RETURN 390 : * 391 : * Format handling 392 : * 393 : * 394 : * Error handling 395 : * 396 : 9999 CONTINUE 397 : IRET=1 398 : WRITE(IOIMP,*) 'An error was detected in subroutine ccgtsu' 399 : RETURN 400 : * 401 : * End of subroutine CCGTSU 402 : * 403 : END 404 : 405 : 406 : 407 : 408 : 409 :
© Cast3M 2003 - All rights reserved.
Disclaimer