1 : C RTENS1 SOURCE PV 09/03/12 21:34:02 6325 2 : SUBROUTINE RTENS1(IPCHE1,IFOMEM,IMOT,IPTV2,IELEME,IVAVEC,IVACOM, 3 : & IVARES,IDEFO,IINTE,MELE,NPINT,NVEC,V1,V2,W2,W3, 4 : & CENTR1,CENTR2,AXEI1,IER1) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8(A-H,O-Z) 7 : *-----------------------------------------------------------------------* 8 : * Operateur RTENS : cas de la formulation massive * 9 : * * 10 : * IPCHE1 (e) pointeur sur un MCHAML de caracteristiques * 11 : * = 0 si isotropie * 12 : * IFOMEM (e) = IFOUR de CCOPTIO * 13 : * IMOT (e) indique le type de repere desire (cf RTENS) * 14 : * IPTV2 (e) pointeur sur le 2nd point repere * 15 : * IELEME (e) pointeur sur le segment MELEME (actif) * 16 : * IVAVEC (e/s) pointeur sur un segment MPTVAL (actif) * 17 : * IVACOM (e/s) pointeur sur un segment MPTVAL (actif) * 18 : * IVARES (e/s) pointeur sur un segment MPTVAL (actif) * 19 : * IDEFO (e) =1 : tenseur de deformations (contraintes sinon) * 20 : * IINTE (e) pointeur sur le segment MINTE (actif) * 21 : * MELE (e) numero de l'element-fini dans NOMTP * 22 : * NPINT (e) nombre de points d'integration (coques) * 23 : * NVEC (e) nombre de composantes du futur MCHAML * 24 : * V1 (e) coordonnees et norme du 1er vecteur * 25 : * V2 (e) coordonnees et norme du 2nd vecteur * 26 : * W2 (e) coordonnees d'un 1er vecteur de travail * 27 : * W3 (e) coordonnees d'un 2nd vecteur de travail * 28 : * CENTR1 (e) coordonnees du 1er point repere * 29 : * CENTR2 (e) coordonnees du 2nd point repere * 30 : * AXEI1 (e) coordonnees du vecteur de l'axe de symetrie * 31 : * IER1 (s) code d'erreur pour desactivation dans RTENS * 32 : * D.R.-M. le 17/3/94 * 33 : *-----------------------------------------------------------------------* 34 : -INC CCOPTIO 35 : -INC CCHAMP 36 : -INC SMCHAML 37 : -INC SMINTE 38 : -INC SMCOORD 39 : -INC SMELEME 40 : * 41 : * MWRK1,3,4 initialises dans RTENS1 42 : * 43 : SEGMENT MWRK1 44 : REAL*8 XEL(3,NBNN),XEL2(3,NBNN) 45 : ENDSEGMENT 46 : * 47 : SEGMENT MWRK3 48 : REAL*8 A(NDIM,NDIM),R(NDIM,NDIM),RT(NDIM,NDIM),TRAV(NDIM,NDIM) 49 : ENDSEGMENT 50 : * 51 : SEGMENT MWRK4 52 : REAL*8 XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) 53 : REAL*8 TXR1(IDIM,IDIM),VALVEC(NVEC) 54 : ENDSEGMENT 55 : * 56 : * Les MPTVAL recueillent les donnees pour le MCHAML resultat 57 : * IVAL contient les pointeurs des MELVAL du nouveau MCHAML 58 : * 59 : SEGMENT MPTVAL 60 : INTEGER IPOS(NS) , NSOF(NS) 61 : INTEGER IVAL(NCOSOU) 62 : CHARACTER*16 TYVAL(NCOSOU) 63 : ENDSEGMENT 64 : * 65 : DIMENSION VECWRK(3),V1(4),V2(4),W2(3),W3(3) 66 : DIMENSION CENTR1(3),CENTR2(3),AXEI1(3),VECX(3),VECY(3) 67 : DIMENSION UR(3),UTHETA(3),UPHI(3),UN(3),UT(3),XIGAU(3) 68 : * 69 : IER1 = 0 70 : BIDON = 0.D0 71 : ERRAXI = 0 72 : MELEME = IELEME 73 : NBNN = NUM(/1) 74 : NBELEM = NUM(/2) 75 : MINTE = IINTE 76 : NBPGAU = POIGAU(/1) 77 : * 78 : NDIM=IDIM 79 : IF (IFOMEM.EQ.1) NDIM=IDIM+1 80 : SEGINI MWRK3 81 : * 82 : IF (IPCHE1.EQ.0.AND.IMOT.EQ.0) THEN 83 : * 84 : * Repere cartesien : on veut le tenseur dans un repere defini 85 : * par un ou deux vecteurs. On construit la matrice de rotation 86 : * qui fait passer du repere general au repere defini par V1 87 : * (et V2 en 3D) 88 : * 89 : IF (IDIM.EQ.2) THEN 90 : R(1,1)=V1(1)/V1(4) 91 : R(2,1)=V1(2)/V1(4) 92 : R(1,2)= -1.D0 * V1(2)/V1(4) 93 : R(2,2)=V1(1)/V1(4) 94 : IF (IFOMEM.EQ.1) THEN 95 : R(3,3) = 1.D0 96 : ENDIF 97 : ELSE IF (IDIM.EQ.3) THEN 98 : IF (IPTV2.EQ.0) THEN 98 : CALL ERREUR(338) 99 : SEGSUP MWRK3 100 : IER1 = 1 101 : RETURN 102 : ENDIF 103 : R(1,1)=V1(1)/V1(4) 104 : R(2,1)=V1(2)/V1(4) 105 : R(3,1)=V1(3)/V1(4) 106 : R(1,2)=W2(1) 107 : R(2,2)=W2(2) 108 : R(3,2)=W2(3) 109 : R(1,3)=W3(1) 110 : R(2,3)=W3(2) 111 : R(3,3)=W3(3) 112 : ENDIF 112 : CALL TRSPOD(R,NDIM,NDIM,RT) 113 : ELSE 114 : * 115 : * Le repere choisi n'est pas cartesien. On recupere les fonctions 116 : * de forme et leurs derivees au centre de l'element pour calculer 117 : * les axes locaux 118 : * 119 : NLG=NUMGEO(MELE) 119 : CALL RESHPT(1,NBNN,NLG,MELE,NPINT,IPT1,IRT1) 120 : MINTE2=IPT1 121 : SEGACT MINTE2 122 : SEGINI MWRK4,MWRK1 123 : ENDIF 124 : * 125 : * Boucle sur les elements 126 : * 127 : DO 1010 IB=1,NBELEM 128 : * 129 : * Recherche des coordonnees des noeuds de l'element IB 130 : * 131 : IF (IMOT.NE.0.OR.IPCHE1.NE.0) 131 : $ CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XEL) 132 : * 133 : IF (IPCHE1.NE.0) THEN 134 : * 135 : * >>> Repere d'Orthotropie <<< 136 : * 137 : * Calcul des axes locaux pour les materiaux orthotropes, 138 : * anisotropes et unidirectionnels 139 : * 140 : NBSH=MINTE2.SHPTOT(/2) 140 : CALL RLOCAL(XEL,MINTE2.SHPTOT,NBSH,NBNN,TXR1) 141 : if (nbsh.eq.-1) then 141 : call erreur(525) 142 : return 143 : endif 144 : IF (IERR.NE.0) THEN 145 : SEGSUP MWRK1,MWRK3,MWRK4 146 : IER1 = 1 147 : RETURN 148 : ENDIF 149 : * 150 : * RECHERCHE DE NBGMAX 151 : * 152 : NBGMAX=1 153 : MPTVAL=IVAVEC 154 : DO 1311 IV=1,NVEC 155 : IF (IVAL(IV).NE.0) THEN 156 : MELVAL=IVAL(IV) 157 : NBGMAX=MAX(NBGMAX,VELCHE(/1)) 158 : ENDIF 159 : 1311 CONTINUE 160 : 161 : ENDIF 162 : 163 : 164 : * 165 : * Boucle sur les points de Gauss 166 : * 167 : 168 : DO 1010 IGAU=1,NBPGAU 169 : 170 : 171 : 172 : * 173 : * >>> Repere d'Orthotropie <<< 174 : * 175 : *------------------------------------------------------------ 176 : * MLR 13/8/99 ON MET LE CALCUL DES AXES DANS LA BOUCLE 177 : * SUR LES POINTS DE GAUSS ET NON PAS EN DEHORS 178 : *------------------------------------------------------------ 179 : 180 : IF (IPCHE1.NE.0) THEN 181 : IF(IGAU.EQ.1.OR.NBGMAX.GT.1) THEN 182 : 183 : MPTVAL=IVAVEC 184 : 185 : DO 1011 IV=1,NVEC 186 : IF (IVAL(IV).NE.0) THEN 187 : MELVAL=IVAL(IV) 188 : IBMN=MIN(IB,VELCHE(/2)) 189 : IGMN=MIN(IGAU,VELCHE(/1)) 190 : VALVEC(IV)=VELCHE(IGMN,IBMN) 191 : ELSE 192 : VALVEC(IV)=0.D0 193 : ENDIF 194 : 1011 CONTINUE 194 : CALL RGLOB(VALVEC,IDIM,TXR1,XLOC,XGLOB) 195 : IF (IERR.NE.0) THEN 196 : SEGSUP MWRK1,MWRK3,MWRK4 197 : IER1 = 1 198 : RETURN 199 : ENDIF 200 : DO 1012 IC=1,IDIM 201 : DO 1012 IL=1,IDIM 202 : R(IL,IC)=XGLOB(IL,IC) 203 : 1012 CONTINUE 204 : IF (IDIM.EQ.2.AND.IFOMEM.EQ.1) R(3,3)=1.D0 204 : CALL TRSPOD(R,NDIM,NDIM,RT) 205 : 206 : ENDIF 207 : ENDIF 208 : 209 : 210 : *------------------------------------------------------------ 211 : 212 : 213 : * 214 : * Sous-zones du MCHAML avant rotation 215 : * 216 : MPTVAL=IVACOM 217 : * 218 : * Initialisations 219 : * 220 : IF (IMOT.NE.0) THEN 221 : DO 1013 IC = 1,IDIM 222 : XIGAU(IC) = 0.D0 223 : DO 1013 IL=1,XEL(/2) 224 : XIGAU(IC)=XIGAU(IC)+(SHPTOT(1,IL,IGAU)*XEL(IC,IL)) 225 : 1013 CONTINUE 226 : * 227 : SCAL=0.D0 228 : DO 1014 IL=1,IDIM 229 : UTHETA(IL)=0.D0 230 : UPHI(IL)=0.D0 231 : UT(IL)=0.D0 232 : VECX(IL)=0.D0 233 : VECY(IL)=0.D0 234 : UR(IL)=XIGAU(IL)-CENTR1(IL) 235 : UN(IL)=UR(IL) 236 : 1014 CONTINUE 237 : DO 1015 IL=1,IDIM 238 : SCAL=SCAL+UR(IL)*V1(IL) 239 : 1015 CONTINUE 240 : DO 1016 IL=1,IDIM 241 : UR(IL)=UR(IL)-SCAL*V1(IL) 242 : 1016 CONTINUE 243 : * 244 : SCAL=0.D0 245 : DO 1017 IL=1,IDIM 246 : SCAL=SCAL+UR(IL)**2 247 : 1017 CONTINUE 248 : SCAL = SQRT(SCAL) 249 : IF (SCAL.EQ.0.D0) THEN 249 : CALL ERREUR(642) 250 : IER1 = 1 251 : GOTO 1010 252 : ENDIF 253 : IF (IDIM.EQ.3) THEN 253 : CALL NORMER(UR) 253 : CALL NORMER(UN) 253 : CALL PROVEC(V1,UR,UTHETA) 254 : ELSE 255 : * 256 : * Dimension 2 : 'POLA' 257 : * 258 : UR(1)=UR(1)/SCAL 259 : UR(2)=UR(2)/SCAL 260 : R(1,1)=UR(1) 261 : R(1,2)=-UR(2) 262 : R(2,2)=UR(1) 263 : R(2,1)=UR(2) 264 : IF (IFOMEM.EQ.1) THEN 265 : R(3,3)=1D0 266 : ENDIF 266 : CALL TRSPOD (R,NDIM,NDIM,RT) 267 : ENDIF 268 : * 269 : * Debut des calculs de R pour les autres reperes 270 : * 271 : * -- Cas CYLINDRIQUE -- 272 : * 273 : IF (IMOT.EQ.2) THEN 274 : DO 1019 IL=1,IDIM 275 : R(IL,1)=UR(IL) 276 : R(IL,2)=UTHETA(IL) 277 : R(IL,3)=V1(IL) 278 : 1019 CONTINUE 278 : CALL TRSPOD (R,NDIM,NDIM,RT) 279 : ELSE 280 : 281 : ENDIF 282 : * 283 : * -- Cas SPHERIQUE -- 284 : * 285 : IF (IMOT.EQ.3) THEN 286 : UR(1)=UN(1) 287 : UR(2)=UN(2) 288 : UR(3)=UN(3) 289 : UPHI(1)=UTHETA(1) 290 : UPHI(2)=UTHETA(2) 291 : UPHI(3)=UTHETA(3) 291 : CALL PROVEC (UPHI,UR,UTHETA) 292 : DO 1021 IL=1,IDIM 293 : R(IL,1)=UR(IL) 294 : R(IL,2)=UTHETA(IL) 295 : R(IL,3)=UPHI(IL) 296 : 1021 CONTINUE 296 : CALL TRSPOD (R,NDIM,NDIM,RT) 297 : ENDIF 298 : * 299 : * -- Cas TORIQUE CIRCULAIRE -- 300 : * 301 : IF (IMOT.EQ.4) THEN 302 : VECWRK(1)=CENTR2(1)-CENTR1(1) 303 : VECWRK(2)=CENTR2(2)-CENTR1(2) 304 : VECWRK(3)=CENTR2(3)-CENTR1(3) 304 : CALL NORME(VECWRK,SCAL) 305 : UN(1)=UN(1)-SCAL*UR(1) 306 : UN(2)=UN(2)-SCAL*UR(2) 307 : UN(3)=UN(3)-SCAL*UR(3) 307 : CALL NORMER(UN) 307 : CALL PROVEC(UN,UTHETA,UT) 308 : DO 1022 IL=1,IDIM 309 : R(IL,1)=UTHETA(IL) 310 : R(IL,2)=UT(IL) 311 : R(IL,3)=UN(IL) 312 : 1022 CONTINUE 312 : CALL TRSPOD (R,NDIM,NDIM,RT) 313 : ENDIF 314 : * 315 : * -- Cas TORIQUE CARTESIEN -- 316 : * 317 : IF (IMOT.EQ.5) THEN 318 : DO 1023 IL=1,IDIM 319 : R(IL,1)=UR(IL) 320 : R(IL,2)=UTHETA(IL) 321 : R(IL,3)=V1(IL) 322 : 1023 CONTINUE 322 : CALL TRSPOD (R,NDIM,NDIM,RT) 323 : ENDIF 324 : ENDIF 325 : * 326 : * Tenseur avant changement de repere 327 : * 328 : MELVAL=IVAL(1) 329 : IGMN = MIN(IGAU,VELCHE(/1)) 330 : IBMN = MIN(IB, VELCHE(/2)) 331 : A(1,1) = VELCHE(IGMN,IBMN) 332 : * 333 : MELVAL=IVAL(2) 334 : IGMN = MIN(IGAU,VELCHE(/1)) 335 : IBMN = MIN(IB, VELCHE(/2)) 336 : A(2,2) = VELCHE(IGMN,IBMN) 337 : * 338 : MELVAL=IVAL(4) 339 : IGMN = MIN(IGAU,VELCHE(/1)) 340 : IBMN = MIN(IB, VELCHE(/2)) 341 : A(1,2) = VELCHE(IGMN,IBMN) 342 : * 343 : IF (IDEFO.EQ.1) A(1,2)=A(1,2)/2.D0 344 : A(2,1)=A(1,2) 345 : * 346 : IF (IFOMEM.LT.1) GOTO 6610 347 : * 348 : MELVAL=IVAL(3) 349 : IGMN = MIN(IGAU,VELCHE(/1)) 350 : IBMN = MIN(IB, VELCHE(/2)) 351 : A(3,3) = VELCHE(IGMN,IBMN) 352 : * 353 : MELVAL=IVAL(5) 354 : IGMN = MIN(IGAU,VELCHE(/1)) 355 : IBMN = MIN(IB, VELCHE(/2)) 356 : A(3,1) = VELCHE(IGMN,IBMN) 357 : * 358 : MELVAL=IVAL(6) 359 : IGMN = MIN(IGAU,VELCHE(/1)) 360 : IBMN = MIN(IB, VELCHE(/2)) 361 : A(3,2) = VELCHE(IGMN,IBMN) 362 : * 363 : IF (IDEFO.EQ.1) A(3,1)=A(3,1)/2.D0 364 : IF (IDEFO.EQ.1) A(3,2)=A(3,2)/2.D0 365 : A(1,3)=A(3,1) 366 : A(2,3)=A(3,2) 367 : * 368 : MELVAL=IVAL(3) 369 : IGMN = MIN(IGAU,VELCHE(/1)) 370 : IBMN = MIN(IB, VELCHE(/2)) 371 : A(3,3) = VELCHE(IGMN,IBMN) 372 : * 373 : 6610 CONTINUE 374 : * 375 : MELVAL=IVAL(3) 376 : IGMN = MIN(IGAU,VELCHE(/1)) 377 : IBMN = MIN(IB, VELCHE(/2)) 378 : AUX = VELCHE(IGMN,IBMN) 379 : * t 380 : * >>> Rotation du tenseur : A = R A R <<< 381 : * 381 : CALL MULMAT(TRAV,A,R,NDIM,NDIM,NDIM) 381 : CALL MULMAT(A,RT,TRAV,NDIM,NDIM,NDIM) 382 : * 383 : * Tenseur apres changement de repere 384 : * Sous-zones du MCHAML resultat 385 : * 386 : MPTVAL=IVARES 387 : * 388 : MELVAL=IVAL(1) 389 : VELCHE(IGAU,IB) = A(1,1) 390 : * 391 : MELVAL=IVAL(2) 392 : VELCHE(IGAU,IB) = A(2,2) 393 : * 394 : IF (IDEFO.EQ.1) A(1,2)=A(1,2)*2.D0 395 : * 396 : MELVAL=IVAL(4) 397 : VELCHE(IGAU,IB) = A(1,2) 398 : * 399 : IF (IFOMEM.LT.1) THEN 400 : * 401 : MELVAL=IVAL(3) 402 : VELCHE(IGAU,IB)= AUX 403 : * 404 : ELSE 405 : * 406 : MELVAL=IVAL(3) 407 : VELCHE(IGAU,IB)=A(3,3) 408 : * 409 : IF (IDEFO.EQ.1) A(3,1)=A(3,1)*2.D0 410 : IF (IDEFO.EQ.1) A(3,2)=A(3,2)*2.D0 411 : * 412 : MELVAL=IVAL(5) 413 : VELCHE(IGAU,IB)= A(3,1) 414 : * 415 : MELVAL=IVAL(6) 416 : VELCHE(IGAU,IB)=A(3,2) 417 : * 418 : ENDIF 419 : * 420 : 1010 CONTINUE 421 : SEGSUP MWRK3 422 : IF (IPCHE1.NE.0.OR.IMOT.NE.0) THEN 423 : SEGSUP MWRK1,MWRK4 424 : SEGDES MINTE2 425 : ENDIF 426 : * 427 : RETURN 428 : END 429 : 430 : 431 : 432 : 433 : 434 : 435 : 436 : 437 :
© Cast3M 2003 - All rights reserved.
Disclaimer