Esope source files

   1 : C CALJBR    SOURCE    MAGN      08/10/10    21:15:00     6185           
   2 :       SUBROUTINE CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
   3 :      *NES,ND,NP,NPG,IAXI,AIRE,AJ,SGN)
   4 : C***********************************************************************
   5 : C     Ce Sp calcule le gradient des fonctions de formes pour un element
   6 : C     courant a partir des gradients de l'element de reference et des
   7 : C     coordonnees des points de l'element
   8 : C     Il calcule le produit PgSq poids d'integration x element d'aire et
   9 : C     le Jacobien
  10 : C     Pour les elements coques ou filaires les normales et tangentes
  11 : C     sont donnees a la place du Jacobien
  12 : C     les fonctions de forme FN sont donnees en entree pour calculer
  13 : C     les elements d'aire en axisymetrique
  14 : C
  15 : C     Entree :
  16 : C     NES         Dimension d'espace de l'element de reference
  17 : C                 (different de ND pour les elements coques ou filaires)
  18 : C     ND          Dimension d'espace
  19 : C     NP          Nombre de noeuds de l'element
  20 : C     NPG         Nombre de points d'integration
  21 : C     IAXI        =0 en 3D et en 2D PLAN
  22 : C     IAXI        NE 0 en axi symetrique   axe de symetrie  x
  23 : C     FN(NP,NPG) Valeurs des fonctions de forme aux points d'integration
  24 : C     GR(NES,NP,NPG) Valeurs des gradients dans l'element de reference
  25 : C     PG(NPG)        Poids d'integration
  26 : C     XYZ(ND,NP)     Coordonnees des noeuds de l'element
  27 : C
  28 : C     Sortie
  29 : C     HR(NES,NP,NPG) Valeurs des gradients dans l'element courant
  30 : C     PGSQ(NPG)      poids d'integration x element d'aire
  31 : C  ATTENTION   en axi : poids d'integration x element d'aire x 2 pi R
  32 : C     RPG(NPG)       Rayon des points d'integration
  33 : C     AIRE           Aire de l'element
  34 : C     AJ(ND,ND,NPG)  Valeurs du jacobien pour chaque point d'integration
  35 : C     SGN Pour les elements massif =1 si meme orientation element de Ref
  36 : C     SGN Pour les elements massif =-1 si orientation opposee
  37 : C     SGN 0 pour les autres types d'elements
  38 : C
  39 : C
  40 : C     CE SP TRAITE LES ELEMENTS VOLUMIQUES, SURFACIQUES et FILAIRES
  41 : C     DANS LES CAS 2D (PLAN et AXI) ET 3D
  42 : C
  43 : C     Dans AJ on stoke l'inverse du Jacobien (AJ=1/J) si l'element est
  44 : C     volumique
  45 : C     Sinon pour les elements coques on stoke tangentes et normales
  46 : C                          AJ=|tx ty|     |tx ty tz|
  47 : C                             |nx ny|  ou |ux uy uz|
  48 : C                                         |nx ny nz|
  49 : C
  50 : C
  51 : C     CALCUL INTERMEDIAIRE DE L'ELEMENT D'AIRE   SQ=DET(J)
  52 : C
  53 : C***********************************************************************
  54 :       IMPLICIT INTEGER(I-N)
  55 :       IMPLICIT REAL*8 (A-H,O-Z)
  56 : C
  57 :       REAL*8 FN(NP,NPG),GR(NES,NP,NPG),HR(ND,NP,NPG)
  58 :       REAL*8 PG(NPG),XYZ(ND,NP),PGSQ(NPG),RPG(NPG),AJ(ND,ND,NPG)
  59 :       REAL*8 AG(3,3,25)
  60 : -INC CCREEL
  61 : C
  62 : C***
  63 :       SGN=0.D0
  63 :       CALL INITD(RPG,NPG,1.D0)
  64 : 
  65 :       IF(NES.EQ.1.AND.ND.EQ.2)THEN
  66 : 
  66 :       IF(NPG.GT.25)CALL ARRET(0)
  67 :       AIRE=0.D0
  68 :       DO 110 L=1,NPG
  69 :       AJX=0.D0
  70 :       AJY=0.D0
  71 :       DO 111 I=1,NP
  72 :       AJX=AJX+GR(1,I,L)*XYZ(1,I)
  73 :       AJY=AJY+GR(1,I,L)*XYZ(2,I)
  74 :   111 CONTINUE
  75 :       AJN=(AJX*AJX+AJY*AJY)**0.5D0
  76 : C     write(6,*)' AJX,AJY=',AJX,AJY,AJN
  77 :       DET=AJX*AJX+AJY*AJY
  78 : 
  79 :       AJ(1,1,L)=AJX/AJN
  80 :       AJ(2,1,L)=AJY/AJN
  81 :       AJ(1,2,L)=-AJY/AJN
  82 :       AJ(2,2,L)=AJX/AJN
  83 : 
  84 :       AG(1,1,L)=AJX/DET
  85 :       AG(2,1,L)=-AJY/DET
  86 :       AG(1,2,L)=AJY/DET
  87 :       AG(2,2,L)=AJX/DET
  88 : 
  89 :       PGSQ(L)=PG(L)*AJN
  90 :       AIRE=AIRE+PGSQ(L)
  91 :   110 CONTINUE
  92 : C     write(6,*)' AJ ds CALJBR',AIRE
  93 : C     write(6,1002)AJ
  94 :       DO 131 L=1,NPG
  95 :       DO 131 I=1,NP
  96 :       DO 131 N=1,ND
  97 :       U=0.D0
  98 :       DO 132 M=1,NES
  99 :   132 U=U+AG(M,N,L)*GR(M,I,L)
 100 :   131 HR(N,I,L)=U
 101 : C     write(6,*)'GR'
 102 : C     write(6,1002)gr
 103 : C     write(6,*)'HR'
 104 : C     write(6,1002)hr
 105 : 
 106 :       ELSEIF(NES.EQ.1.AND.ND.EQ.3)THEN
 107 : 
 107 :       IF(NPG.GT.25)CALL ARRET(0)
 108 :       AIRE=0.D0
 109 :       DO 310 L=1,NPG
 110 :       AJX=0.D0
 111 :       AJY=0.D0
 112 :       AJZ=0.D0
 113 :       BJX=0.D0
 114 :       BJY=0.D0
 115 :       BJZ=0.D0
 116 :       DO 311 I=1,NP
 117 :       AJX=AJX+GR(1,I,L)*XYZ(1,I)
 118 :       AJY=AJY+GR(1,I,L)*XYZ(2,I)
 119 :       AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
 120 : C     BJX=BJX+GR(2,I,L)*XYZ(1,I)
 121 : C     BJY=BJY+GR(2,I,L)*XYZ(2,I)
 122 : C     BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
 123 :   311 CONTINUE
 124 : C     write(6,*)ajx,ajy,ajz
 125 :       IF(ABS(AJX).NE.0.D0.AND.ABS(AJY).NE.0.D0)THEN
 126 :       BJX=AJY
 127 :       BJY=-AJX
 128 :       BJZ=AJZ
 129 :       ELSE
 130 :       BJX=1.D0
 131 :       BJY=1.D0
 132 :       BJZ=0.D0
 133 :       ENDIF
 134 : 
 135 :       XB=AJY*BJZ-AJZ*BJY
 136 :       YB=AJZ*BJX-AJX*BJZ
 137 :       ZB=AJX*BJY-AJY*BJX
 138 : 
 139 :       AJN=(XB*XB+YB*YB*ZB*ZB)**0.5D0+1.D-5
 140 : C     write(6,*)' AJN=',AJN
 141 :       PGSQ(L)=PG(L)*AJN
 142 :       AIRE=AIRE+PGSQ(L)
 143 : 
 144 :       AJ(1,1,L)=AJX/AJN
 145 :       AJ(2,1,L)=AJY/AJN
 146 :       AJ(3,1,L)=AJZ/AJN
 147 :       AJ(1,2,L)=BJX/AJN
 148 :       AJ(2,2,L)=BJY/AJN
 149 :       AJ(3,2,L)=BJZ/AJN
 150 :       AJ(1,3,L)=XB/AJN
 151 :       AJ(2,3,L)=YB/AJN
 152 :       AJ(3,3,L)=ZB/AJN
 153 : 
 154 :       AG(1,1,L)=AJX
 155 :       AG(2,1,L)=AJY
 156 :       AG(3,1,L)=AJZ
 157 :       AG(1,2,L)=BJX
 158 :       AG(2,2,L)=BJY
 159 :       AG(3,2,L)=BJZ
 160 :       AG(1,3,L)=XB
 161 :       AG(2,3,L)=YB
 162 :       AG(3,3,L)=ZB
 163 : 
 164 :       D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
 165 :       D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
 166 :       D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
 167 :       D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
 168 :       D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
 169 :       D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
 170 :       D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
 171 :       D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
 172 :       D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
 173 :       VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
 174 :       RVINT=1.D0/VINT
 175 :       AG(1,1,L)= RVINT*D11
 176 :       AG(1,2,L)=-RVINT*D12
 177 :       AG(1,3,L)= RVINT*D13
 178 :       AG(2,1,L)=-RVINT*D21
 179 :       AG(2,2,L)= RVINT*D22
 180 :       AG(2,3,L)=-RVINT*D23
 181 :       AG(3,1,L)= RVINT*D31
 182 :       AG(3,2,L)=-RVINT*D32
 183 :       AG(3,3,L)= RVINT*D33
 184 : 
 185 :   310 CONTINUE
 186 : C     write(6,*)' AJ ds CALJBR'
 187 : C     write(6,1002)AJ
 188 :       DO 331 L=1,NPG
 189 :       DO 331 I=1,NP
 190 :       DO 331 N=1,ND
 191 :       U=0.D0
 192 :       DO 332 M=1,NES
 193 :   332 U=U+AG(M,N,L)*GR(M,I,L)
 194 :   331 HR(N,I,L)=U
 195 : C     write(6,*)'GR'
 196 : C     write(6,1002)gr
 197 : C     write(6,*)'HR'
 198 : C     write(6,1002)hr
 199 : 
 200 :       ELSEIF(NES.EQ.2.AND.ND.EQ.3)THEN
 201 : C     write(6,*)' Je passe ICI'
 202 :                                                                                 
 203 :       AIRE=0.D0
 204 :       DO 210 L=1,NPG
 205 :       AJX=0.D0
 206 :       AJY=0.D0
 207 :       AJZ=0.D0
 208 :       BJX=0.D0
 209 :       BJY=0.D0
 210 :       BJZ=0.D0
 211 :       DO 211 I=1,NP
 212 :       AJX=AJX+GR(1,I,L)*XYZ(1,I)
 213 :       AJY=AJY+GR(1,I,L)*XYZ(2,I)
 214 :       AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
 215 :       BJX=BJX+GR(2,I,L)*XYZ(1,I)
 216 :       BJY=BJY+GR(2,I,L)*XYZ(2,I)
 217 :       BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
 218 :   211 CONTINUE
 219 : 
 220 :       XB=AJY*BJZ-AJZ*BJY
 221 :       YB=AJZ*BJX-AJX*BJZ
 222 :       ZB=AJX*BJY-AJY*BJX
 223 : 
 224 :       AJN=(XB*XB+YB*YB+ZB*ZB)**0.5D0
 225 : 
 226 :       PGSQ(L)=PG(L)*AJN
 227 :       AIRE=AIRE+PGSQ(L)
 228 : 
 229 :       AJ(1,1,L)=AJX/AJN
 230 :       AJ(2,1,L)=AJY/AJN
 231 :       AJ(3,1,L)=AJZ/AJN
 232 :       AJ(1,2,L)=BJX/AJN
 233 :       AJ(2,2,L)=BJY/AJN
 234 :       AJ(3,2,L)=BJZ/AJN
 235 :       AJ(1,3,L)=XB/AJN
 236 :       AJ(2,3,L)=YB/AJN
 237 :       AJ(3,3,L)=ZB/AJN
 238 : 
 239 :       AG(1,1,L)=AJX
 240 :       AG(2,1,L)=AJY
 241 :       AG(3,1,L)=AJZ
 242 :       AG(1,2,L)=BJX
 243 :       AG(2,2,L)=BJY
 244 :       AG(3,2,L)=BJZ
 245 :       AG(1,3,L)=XB
 246 :       AG(2,3,L)=YB
 247 :       AG(3,3,L)=ZB
 248 : 
 249 :       D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
 250 :       D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
 251 :       D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
 252 :       D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
 253 :       D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
 254 :       D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
 255 :       D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
 256 :       D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
 257 :       D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
 258 :       VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
 259 :       RVINT=1.D0/VINT
 260 :       AG(1,1,L)= RVINT*D11
 261 :       AG(1,2,L)=-RVINT*D12
 262 :       AG(1,3,L)= RVINT*D13
 263 :       AG(2,1,L)=-RVINT*D21
 264 :       AG(2,2,L)= RVINT*D22
 265 :       AG(2,3,L)=-RVINT*D23
 266 :       AG(3,1,L)= RVINT*D31
 267 :       AG(3,2,L)=-RVINT*D32
 268 :       AG(3,3,L)= RVINT*D33
 269 : 
 270 :  210  CONTINUE
 271 : C                                                                               
 272 :       DO 231 L=1,NPG
 273 :       DO 231 I=1,NP
 274 :       DO 231 N=1,ND
 275 :       U=0.D0
 276 :       DO 232 M=1,NES
 277 :   232 U=U+AG(M,N,L)*GR(M,I,L)
 278 :   231 HR(N,I,L)=U
 279 : C     write(6,*)'GR'
 280 : C     write(6,1002)gr
 281 : C     write(6,*)'HR'
 282 : C     write(6,1002)hr
 283 :                                                                                 
 284 :       ELSE
 285 : 
 286 :       DO 10 L=1,NPG
 287 :       DO 10 M=1,ND
 288 :       DO 10 N=1,ND
 289 :       AJT=0.D0
 290 :       DO 11 I=1,NP
 291 :       AJT=AJT+GR(M,I,L)*XYZ(N,I)
 292 :    11 CONTINUE
 293 :       AJ(N,M,L)=AJT
 294 :    10 CONTINUE
 295 : 
 296 : C
 297 :       DO 20 L=1,NPG
 298 :       IF(ND.EQ.1)THEN
 299 :       VINT=AJ(1,1,L)
 300 : C     VINT=ABS(VINT)
 301 :       AJ(1,1,L)=1.D0/VINT
 302 :       ELSEIF(ND.EQ.2)THEN
 303 :       VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
 304 : C     VINT=ABS(VINT)
 305 :       RVINT=1.D0/VINT
 306 :       D11=AJ(2,2,L)
 307 :       D12=AJ(1,2,L)
 308 :       D21=AJ(2,1,L)
 309 :       D22=AJ(1,1,L)
 310 :       AJ(1,1,L)= RVINT*D11
 311 :       AJ(1,2,L)=-RVINT*D12
 312 :       AJ(2,1,L)=-RVINT*D21
 313 :       AJ(2,2,L)= RVINT*D22
 314 :       ELSEIF(ND.EQ.3)THEN
 315 :       D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L)
 316 :       D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L)
 317 :       D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L)
 318 :       D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L)
 319 :       D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L)
 320 :       D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L)
 321 :       D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L)
 322 :       D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L)
 323 :       D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L)
 324 :       VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31
 325 : C     VINT=ABS(VINT)
 326 :       RVINT=1.D0/VINT
 327 :       AJ(1,1,L)= RVINT*D11
 328 :       AJ(1,2,L)=-RVINT*D12
 329 :       AJ(1,3,L)= RVINT*D13
 330 :       AJ(2,1,L)=-RVINT*D21
 331 :       AJ(2,2,L)= RVINT*D22
 332 :       AJ(2,3,L)=-RVINT*D23
 333 :       AJ(3,1,L)= RVINT*D31
 334 :       AJ(3,2,L)=-RVINT*D32
 335 :       AJ(3,3,L)= RVINT*D33
 336 :       ELSE
 337 :       VINT=0.D0
 338 :       ENDIF
 339 : 
 340 :       PGSQ(L)=VINT
 340 : C***  CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER)
 341 : C
 342 :    20 CONTINUE
 343 : C
 344 :       DO 31 L=1,NPG
 345 :       DO 31 I=1,NP
 346 :       DO 31 N=1,ND
 347 :       U=0.D0
 348 :       DO 32 M=1,ND
 349 :    32 U=U+AJ(M,N,L)*GR(M,I,L)
 350 :    31 HR(N,I,L)=U
 351 :       U=0.D0
 352 :       V=0.D0
 353 :       DO 4 L=1,NPG
 354 :       W=ABS(PGSQ(L))
 355 :       U=U+PG(L)*W
 356 :       V=V+PG(L)*PGSQ(L)
 357 :     4 PGSQ(L)=PG(L)*W
 358 :       AIRE=U
 359 :       SGN=SIGN(1.D0,V)
 360 : 
 361 :       ENDIF
 362 : 
 363 :       IF(IAXI.EQ.0)RETURN
 364 : C
 365 :       DO 45 L=1,NPG
 366 :       RPGT  =0.D0
 367 :       DO 46 I=1,NP
 368 :       RPGT=RPGT+XYZ(1,I)*FN(I,L)
 369 :  46   CONTINUE
 370 :       RPG(L)=RPGT
 371 :  45   CONTINUE
 372 : C
 373 :       DEUPI=2.D0*XPI
 374 :       U=0.D0
 375 :       DO 47 L=1,NPG
 376 : C?    PGSQ(L)=PGSQ(L)*DEUPI*RPG(L)
 377 :       U=U+PGSQ(L)*DEUPI*RPG(L)
 378 :  47   CONTINUE
 379 :       AIRE=U
 380 : C
 381 :       RETURN
 382 :  1002 FORMAT(10(1X,1PE11.4))
 383 :  1001 FORMAT(20(1X,I5))
 384 :       END
 385 : 
 386 :  
 387 :  
 388 :  
 389 :  
 390 :  
 391 :  
 392 :  
 393 :  
 394 :  
 395 :  
 396 :