1 : C TAILCA SOURCE FANDEUR 10/08/30 21:15:47 6729 2 : C responsable : Mr MILLARD 3 : 4 : C======================================================================= 5 : C 6 : SUBROUTINE TAILCA(IPTRA,POI,IDIM,IFOUR,NBNO,MELE,IELE,MFR, IRET) 7 : C 8 : C======================================================================= 9 : C 10 : C CALCULE LE GRADIENT DE LA GEOMETRIE DEFORMEE 11 : C CALCULE LA TRANSFORMEE DES VECTEURS UNITAIRES DE LA 12 : C GEOMETRIE DE REFERENCE DANS LE REPERE GLOBAL 13 : C ENTREES 14 : C IPTRA=POINTEUR SUR UN SEGMENT CONTENANT LES INFOS : 15 : C -XEL(3,NBNO)=COORDONNEES LOCALES DE L ELEMENT 16 : C -SHP(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE 17 : C -SHPQSI(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (1 0 0) 18 : C _SHPETA(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (0 1 0) 19 : C -SHPDZE(6,NBNO)=DERIVEES PAR RAPPORT A LA GEOMETRIE DE REFERENCE de (0 0 1) 20 : C POI = POI POUR L INTEGRATION DU GRADIANT AU PT DE GAUSS 21 : C IDIM=DIMENSION 22 : C NBNO=NOMBRE DE NOEUDS 23 : C 24 : C SORTIES 25 : C VCOMP = COMPOSANTES DU PARAMETRE DE TAILLE DE L ELEMENT 26 : C ET DU POINT DE GAUSS COURANT 27 : C================================================================= 28 : C 29 : IMPLICIT INTEGER(I-N) 30 : IMPLICIT REAL*8 (A-H,O-Z) 31 : C 32 : PARAMETER (XZER=0.D0,UNDEMI=.5D0,UN=1.D0,DEUX=2.D0, 33 : & x1sur3=0.33333333333333333333333333333333333333333, 34 : & x1sra2=0.70710678118654752440084436210484904) 35 : C Note : x1sur3 = 1./3. et x1sra2 = 1./SQRT(2.) 36 : C 37 : SEGMENT MTRA1 38 : REAL*8 XEL(3,NBNN) 39 : REAL*8 VCOMP(NCOMP) 40 : REAL*8 SHP(6,NBNN),SHPZER(6,NBNN) 41 : REAL*8 SHPQSI(6,NBNN),SHPETA(6,NBNN),SHPDZE(6,NBNN) 42 : C* REAL*8 SHPGAU(6,NBNN) 43 : ENDSEGMENT 44 : 45 : NBNN=NBNO 46 : MTRA1=IPTRA 47 : 48 : C============================= 49 : C CAS 2 DIMENSIONS 50 : C OU CAS 3D COQUES MINCES 51 : C============================= 52 : IF (IDIM.EQ.3.AND.(MFR.EQ.1.OR.MFR.EQ.31)) GOTO 30 53 : C============================= 54 : 55 : POIG=ABS(POI) 56 : 57 : DXDQSI = XZER 58 : DXDETA = XZER 59 : DYDQSI = XZER 60 : DYDETA = XZER 61 : XORI = XZER 62 : YORI = XZER 63 : XQSI = XZER 64 : YQSI = XZER 65 : XETA = XZER 66 : YETA = XZER 67 : 68 : DO 100 I=1,NBNN 69 : DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) 70 : DXDETA=DXDETA+SHP(3,I)*XEL(1,I) 71 : DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) 72 : DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 73 : XORI=XORI+SHPZER(1,I)*XEL(1,I) 74 : YORI=YORI+SHPZER(1,I)*XEL(2,I) 75 : XQSI=XQSI+SHPQSI(1,I)*XEL(1,I) 76 : YQSI=YQSI+SHPQSI(1,I)*XEL(2,I) 77 : XETA=XETA+SHPETA(1,I)*XEL(1,I) 78 : YETA=YETA+SHPETA(1,I)*XEL(2,I) 79 : 100 CONTINUE 80 : 81 : UX=XQSI-XORI 82 : UY=YQSI-YORI 83 : VX=XETA-XORI 84 : VY=YETA-YORI 85 : WX=XQSI-XETA 86 : WY=YQSI-YETA 87 : 88 : RNU=SQRT(Ux*Ux+Uy*Uy) 89 : RNV=SQRT(Vx*Vx+Vy*Vy) 90 : RNW=SQRT(Wx*Wx+Wy*Wy) 91 : C En cas de vecteur(s) nul(s), on met sa(leur) norme(s) a un. 92 : C Le determinant "det" calcule plus loin sera nul d'où une erreur ! 93 : IF (RNU.EQ.XZER) RNU=UN 94 : IF (RNV.EQ.XZER) RNV=UN 95 : IF (RNW.EQ.XZER) RNW=UN 96 : C 97 : C------Calcule de la norme du gradient dans les directions associees----- 98 : C------------------------------a QSI ETA DZE----------------------------- 99 : C------------------et integration par la methode de Gauss---------------- 100 : C 101 : XLQSI=SQRT(POIG*(DXDQSI*DXDQSI+DYDQSI*DYDQSI)) 102 : YLETA=SQRT(POIG*(DXDETA*DXDETA+DYDETA*DYDETA)) 103 : ZLDZE=UN 104 : C 105 : C======================================================================== 106 : C--------------CALCUL DU TENSEUR DE PARAMETRE TAILLE--------------------- 107 : C======================================================================== 108 : C 109 : C----------CALCUL DE LA MATRICE DE PASSAGE ET DE SON INVERSE------------- 110 : C 111 : C========================================================== 112 : C---------------lorsqu il s agit de triangle,-------------- 113 : C---------- les deux + petit cotes sont pris en compte----- 114 : C========================================================== 115 : IF (IELE.EQ.4.OR.IELE.EQ.6) THEN 116 : IF (IELE.EQ.4) THEN 117 : rconst = x1sra2 118 : ELSE 119 : C ELSE IF (IELE.EQ.6) THEN 120 : rconst = UNDEMI 121 : ENDIF 122 : IF (RNW.LT.RNU.OR.RNW.LT.RNV) THEN 123 : IF (RNU.LT.RNV) THEN 124 : YLETA=SQRT(POIG)*RNW 125 : UUX=UX 126 : UUY=UY 127 : VVX=WX 128 : VVY=WY 129 : ELSE 130 : XLQSI=SQRT(POIG)*RNW 131 : UUX=WX 132 : UUY=WY 133 : VVX=VX 134 : VVY=VY 135 : ENDIF 136 : RNUU=SQRT(UUX*UUX+UUY*UUY) 137 : RNVV=SQRT(VVX*VVX+VVY*VVY) 138 : a=UUX/RNUU 139 : b=UUY/RNUU 140 : c=VVX/RNVV 141 : d=VVY/RNVV 142 : ELSE 143 : a=UX/RNU 144 : b=UY/RNU 145 : c=VX/RNV 146 : d=VY/RNV 147 : ENDIF 148 : XLQSI = XLQSI*rconst 149 : YLETA = YLETA*rconst 150 : C------------------------------------------ 151 : C pas de modification pour les quadrangles 152 : C------------------------------------------ 153 : ELSE 154 : a=UX/RNU 155 : b=UY/RNU 156 : c=VX/RNV 157 : d=VY/RNV 158 : ENDIF 159 : C 160 : C------------------------------------------- 161 : C calcul du produit scalaire 162 : C------------------------------------------- 163 : C 164 : XY=UNDEMI*(a*c+b*d)*SQRT(XLQSI*YLETA) 165 : C 166 : C----------------------------------------------------------------- 167 : C CALCUL DES COEFFICIENTS DANS LE REPERE GLOBAL 168 : C----------------------------------------------------------------- 169 : C write (*,*) 'XLQSI : ',XLQSI,' YLETA : ',YLETA 170 : C write (*,*) 171 : C write (*,*) 'RNU : ',RNU,' RNV : ',RNV,' RNW : ',RNW 172 : C write (*,*) 'a : ',a,' c : ',c 173 : C write (*,*) 'b : ',b,' d : ',d 174 : C write (*,*) '=================================================' 175 : C 176 : det=a*d-b*c 177 : IF (det.EQ.XZER) GOTO 999 178 : C 179 : XLXX=(d*d*XLQSI+b*b*YLETA-DEUX*b*d*XY)/(det*det) 180 : XLXY=(-c*d*XLQSI-b*a*YLETA+(b*c+a*d)*XY)/(det*det) 181 : C-------------------------------------------- 182 : C XLYX=XLXY 183 : C-------------------------------------------- 184 : YLYY=(c*c*XLQSI+a*a*YLETA-DEUX*a*c*XY)/(det*det) 185 : C------------------------------------------ 186 : C dans la direction normal au plan 187 : C il n y a pas de probleme d objectivite 188 : C a l utilisateur de rentrer epsi rupture 189 : C------------------------------------------ 190 : ZLZZ=ZLDZE 191 : PXX=(d*d+b*b)/(det*det) 192 : PXY=-(c*d+b*a)/(det*det) 193 : C-------------------------------------------- 194 : C PYX=PXY 195 : C-------------------------------------------- 196 : PYY=(a*a+c*c)/(det*det) 197 : 198 : VCOMP(1)=XLXX 199 : VCOMP(2)=YLYY 200 : VCOMP(3)=ZLZZ 201 : VCOMP(4)=XLXY 202 : VCOMP(5)=PXX 203 : VCOMP(6)=PYY 204 : VCOMP(7)=PXY 205 : 206 : GOTO 666 207 : 208 : C================================ 209 : C CAS TRIDIMENSIONEL 210 : C================================ 211 : 30 CONTINUE 212 : IF (IDIM.NE.3) GOTO 666 213 : IF (MFR.NE.1.AND.MFR.NE.31) GOTO 666 214 : 215 : POIG=ABS(POI)**x1sur3 216 : 217 : DXDQSI=XZER 218 : DXDETA=XZER 219 : DXDDZE=XZER 220 : DYDQSI=XZER 221 : DYDETA=XZER 222 : DYDDZE=XZER 223 : DZDQSI=XZER 224 : DZDETA=XZER 225 : DZDDZE=XZER 226 : 227 : XORI=XZER 228 : YORI=XZER 229 : ZORI=XZER 230 : XQSI=XZER 231 : YQSI=XZER 232 : ZQSI=XZER 233 : XETA=XZER 234 : YETA=XZER 235 : ZETA=XZER 236 : XDZE=XZER 237 : YDZE=XZER 238 : ZDZE=XZER 239 : 240 : DO 300 I=1,NBNO 241 : C 242 : DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) 243 : DXDETA=DXDETA+SHP(3,I)*XEL(1,I) 244 : DXDDZE=DXDDZE+SHP(4,I)*XEL(1,I) 245 : DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) 246 : DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 247 : DYDDZE=DYDDZE+SHP(4,I)*XEL(2,I) 248 : DZDQSI=DZDQSI+SHP(2,I)*XEL(3,I) 249 : DZDETA=DZDETA+SHP(3,I)*XEL(3,I) 250 : DZDDZE=DZDDZE+SHP(4,I)*XEL(3,I) 251 : 252 : XORI=XORI+SHPZER(1,I)*XEL(1,I) 253 : YORI=YORI+SHPZER(1,I)*XEL(2,I) 254 : ZORI=ZORI+SHPZER(1,I)*XEL(3,I) 255 : XQSI=XQSI+SHPQSI(1,I)*XEL(1,I) 256 : YQSI=YQSI+SHPQSI(1,I)*XEL(2,I) 257 : ZQSI=ZQSI+SHPQSI(1,I)*XEL(3,I) 258 : XETA=XETA+SHPETA(1,I)*XEL(1,I) 259 : YETA=YETA+SHPETA(1,I)*XEL(2,I) 260 : ZETA=ZETA+SHPETA(1,I)*XEL(3,I) 261 : XDZE=XDZE+SHPDZE(1,I)*XEL(1,I) 262 : YDZE=YDZE+SHPDZE(1,I)*XEL(2,I) 263 : ZDZE=ZDZE+SHPDZE(1,I)*XEL(3,I) 264 : 265 : 300 CONTINUE 266 : C 267 : C------Calcule de la norme du gradient dans les directions associees----- 268 : C------------------------------a QSI ETA DZE----------------------------- 269 : C------------------et integration par la methode de Gauss---------------- 270 : C 271 : XLQSI = POIG * SQRT(DXDQSI*DXDQSI+DYDQSI*DYDQSI+DZDQSI*DZDQSI) 272 : YLETA = POIG * SQRT(DXDETA*DXDETA+DYDETA*DYDETA+DZDETA*DZDETA) 273 : ZLDZE = POIG * SQRT(DXDDZE*DXDDZE+DYDDZE*DYDDZE+DZDDZE*DZDDZE) 274 : C 275 : C======================================================================== 276 : C--------------CALCUL DU TENSEUR DE PARAMETRE TAILLE--------------------- 277 : C======================================================================== 278 : UX=XQSI-XORI 279 : UY=YQSI-YORI 280 : UZ=ZQSI-ZORI 281 : VX=XETA-XORI 282 : VY=YETA-YORI 283 : VZ=ZETA-ZORI 284 : WX=XDZE-XORI 285 : WY=YDZE-YORI 286 : WZ=ZDZE-ZORI 287 : C 288 : VQSI=SQRT(UX*UX+UY*UY+UZ*UZ) 289 : VETA=SQRT(VX*VX+VY*VY+VZ*VZ) 290 : VDZE=SQRT(WX*WX+WY*WY+WZ*WZ) 291 : C En cas de vecteur(s) nul(s), on met sa(leur) norme(s) a un. 292 : C Le determinant "det" calcule plus loin sera nul d'où une erreur ! 293 : IF (VQSI.EQ.XZER) VQSI=UN 294 : IF (VETA.EQ.XZER) VETA=UN 295 : IF (VDZE.EQ.XZER) VDZE=UN 296 : 297 : C----------CALCUL DE LA MATRICE DE PASSAGE ET DE SON INVERSE------------- 298 : a=UX/VQSI 299 : b=UY/VQSI 300 : c=UZ/VQSI 301 : d=VX/VETA 302 : e=VY/VETA 303 : f=VZ/VETA 304 : o=WX/VDZE 305 : p=WY/VDZE 306 : q=WZ/VDZE 307 : C 308 : C----------------------------------------- 309 : C calcul des demi produits scalaires 310 : C----------------------------------------- 311 : C 312 : XY=UNDEMI*(a*d+b*e+c*f)*SQRT(XLQSI*YLETA) 313 : XZ=UNDEMI*(a*o+b*p+c*q)*SQRT(XLQSI*ZLDZE) 314 : YZ=UNDEMI*(d*o+e*p+f*q)*SQRT(YLETA*ZLDZE) 315 : C 316 : C--------------------------------------------------------------- 317 : C CALCUL DES COEFFICIENTS DANS LE REPERE GLOBAL 318 : C--------------------------------------------------------------- 319 : det=a*e*q+b*f*o+c*d*p-o*e*c-d*b*q-a*p*f 320 : C 321 : if (det.EQ.XZER) GOTO 999 322 : C 323 : C----calcul de l inverse de la matrice de passage---- 324 : C 325 : r=(e*q-f*p)/det 326 : s=(c*p-b*q)/det 327 : t=(b*f-c*e)/det 328 : u=(f*o-d*q)/det 329 : v=(a*q-c*o)/det 330 : w=(c*d-a*f)/det 331 : x=(d*p-o*e)/det 332 : y=(b*o-a*p)/det 333 : z=(a*e-b*d)/det 334 : C 335 : XLXX= XLQSI*r*r+YLETA*s*s+ZLDZE*t*t 336 : & +DEUX*(r*s*XY+s*t*YZ+r*t*XZ) 337 : C 338 : XLXY=XLQSI*r*u+YLETA*s*v+ZLDZE*t*w 339 : & +(s*u+r*v)*XY+(t*u+r*w)*XZ+(t*v+s*w)*YZ 340 : C 341 : XLXZ=XLQSI*r*x+YLETA*s*y+ZLDZE*t*z 342 : & +(s*x+r*y)*XY+(t*x+r*z)*XZ+(t*y+s*z)*YZ 343 : C 344 : YLYY=XLQSI*u*u+YLETA*v*v+ZLDZE*w*w 345 : & +DEUX*(u*v*XY+u*w*XZ+v*w*YZ) 346 : C 347 : YLYZ=XLQSI*x*u+YLETA*y*v+ZLDZE*z*w 348 : & +(x*v+y*u)*XY+(x*w+z*u)*XZ+(y*w+z*v)*YZ 349 : C 350 : ZLZZ=XLQSI*x*x+YLETA*y*y+ZLDZE*z*z 351 : & +DEUX*(x*y*XY+x*z*XZ+y*z*YZ) 352 : C 353 : PXX=r*r+s*s+t*t 354 : C 355 : PXY=r*u+s*v+t*w 356 : C 357 : PXZ=r*x+s*y+t*z 358 : C 359 : PYY=u*u+v*v+w*w 360 : C 361 : PYZ=x*u+y*v+z*w 362 : C 363 : PZZ=x*x+y*y+z*z 364 : C 365 : VCOMP( 1)=XLXX 366 : VCOMP( 2)=YLYY 367 : VCOMP( 3)=ZLZZ 368 : VCOMP( 4)=XLXY 369 : VCOMP( 5)=XLXZ 370 : VCOMP( 6)=YLYZ 371 : VCOMP( 7)=PXX 372 : VCOMP( 8)=PYY 373 : VCOMP( 9)=PZZ 374 : VCOMP(10)=PXY 375 : VCOMP(11)=PXZ 376 : VCOMP(12)=PYZ 377 : 378 : 666 CONTINUE 379 : IRET=1 380 : RETURN 381 : 382 : C------------------------- 383 : C DETERMINANT NUL 384 : C------------------------- 385 : 999 CONTINUE 385 : CALL ERREUR(664) 386 : IRET=0 387 : 388 : RETURN 389 : END 390 : 391 :
© Cast3M 2003 - All rights reserved.
Disclaimer