Télécharger la plaquette Cast3M


fr en

Télécharger caljbr.eso

Retour à la liste

Numérotation des lignes : 

caljbr

  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
  64.       CALL INITD(RPG,NPG,1.D0)
  65.  
  66.       IF(NES.EQ.1.AND.ND.EQ.2)THEN
  67.  
  68.       IF(NPG.GT.25)CALL ARRET(0)
  69.       AIRE=0.D0
  70.       DO 110 L=1,NPG
  71.       AJX=0.D0
  72.       AJY=0.D0
  73.       DO 111 I=1,NP
  74.       AJX=AJX+GR(1,I,L)*XYZ(1,I)
  75.       AJY=AJY+GR(1,I,L)*XYZ(2,I)
  76.   111 CONTINUE
  77.       AJN=(AJX*AJX+AJY*AJY)**0.5D0
  78. C     write(6,*)' AJX,AJY=',AJX,AJY,AJN
  79.       DET=AJX*AJX+AJY*AJY
  80.  
  81.       AJ(1,1,L)=AJX/AJN
  82.       AJ(2,1,L)=AJY/AJN
  83.       AJ(1,2,L)=-AJY/AJN
  84.       AJ(2,2,L)=AJX/AJN
  85.  
  86.       AG(1,1,L)=AJX/DET
  87.       AG(2,1,L)=-AJY/DET
  88.       AG(1,2,L)=AJY/DET
  89.       AG(2,2,L)=AJX/DET
  90.  
  91.       PGSQ(L)=PG(L)*AJN
  92.       AIRE=AIRE+PGSQ(L)
  93.   110 CONTINUE
  94. C     write(6,*)' AJ ds CALJBR',AIRE
  95. C     write(6,1002)AJ
  96.       DO 131 L=1,NPG
  97.       DO 131 I=1,NP
  98.       DO 131 N=1,ND
  99.       U=0.D0
  100.       DO 132 M=1,NES
  101.   132 U=U+AG(M,N,L)*GR(M,I,L)
  102.   131 HR(N,I,L)=U
  103. C     write(6,*)'GR'
  104. C     write(6,1002)gr
  105. C     write(6,*)'HR'
  106. C     write(6,1002)hr
  107.  
  108.       ELSEIF(NES.EQ.1.AND.ND.EQ.3)THEN
  109.  
  110.       IF(NPG.GT.25)CALL ARRET(0)
  111.       AIRE=0.D0
  112.       DO 310 L=1,NPG
  113.       AJX=0.D0
  114.       AJY=0.D0
  115.       AJZ=0.D0
  116.       BJX=0.D0
  117.       BJY=0.D0
  118.       BJZ=0.D0
  119.       DO 311 I=1,NP
  120.       AJX=AJX+GR(1,I,L)*XYZ(1,I)
  121.       AJY=AJY+GR(1,I,L)*XYZ(2,I)
  122.       AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
  123. C     BJX=BJX+GR(2,I,L)*XYZ(1,I)
  124. C     BJY=BJY+GR(2,I,L)*XYZ(2,I)
  125. C     BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
  126.   311 CONTINUE
  127. C     write(6,*)ajx,ajy,ajz
  128.       IF(ABS(AJX).NE.0.D0.AND.ABS(AJY).NE.0.D0)THEN
  129.       BJX=AJY
  130.       BJY=-AJX
  131.       BJZ=AJZ
  132.       ELSE
  133.       BJX=1.D0
  134.       BJY=1.D0
  135.       BJZ=0.D0
  136.       ENDIF
  137.  
  138.       XB=AJY*BJZ-AJZ*BJY
  139.       YB=AJZ*BJX-AJX*BJZ
  140.       ZB=AJX*BJY-AJY*BJX
  141.  
  142.       AJN=(XB*XB+YB*YB*ZB*ZB)**0.5D0+1.D-5
  143. C     write(6,*)' AJN=',AJN
  144.       PGSQ(L)=PG(L)*AJN
  145.       AIRE=AIRE+PGSQ(L)
  146.  
  147.       AJ(1,1,L)=AJX/AJN
  148.       AJ(2,1,L)=AJY/AJN
  149.       AJ(3,1,L)=AJZ/AJN
  150.       AJ(1,2,L)=BJX/AJN
  151.       AJ(2,2,L)=BJY/AJN
  152.       AJ(3,2,L)=BJZ/AJN
  153.       AJ(1,3,L)=XB/AJN
  154.       AJ(2,3,L)=YB/AJN
  155.       AJ(3,3,L)=ZB/AJN
  156.  
  157.       AG(1,1,L)=AJX
  158.       AG(2,1,L)=AJY
  159.       AG(3,1,L)=AJZ
  160.       AG(1,2,L)=BJX
  161.       AG(2,2,L)=BJY
  162.       AG(3,2,L)=BJZ
  163.       AG(1,3,L)=XB
  164.       AG(2,3,L)=YB
  165.       AG(3,3,L)=ZB
  166.  
  167.       D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
  168.       D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
  169.       D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
  170.       D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
  171.       D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
  172.       D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
  173.       D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
  174.       D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
  175.       D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
  176.       VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
  177.       RVINT=1.D0/VINT
  178.       AG(1,1,L)= RVINT*D11
  179.       AG(1,2,L)=-RVINT*D12
  180.       AG(1,3,L)= RVINT*D13
  181.       AG(2,1,L)=-RVINT*D21
  182.       AG(2,2,L)= RVINT*D22
  183.       AG(2,3,L)=-RVINT*D23
  184.       AG(3,1,L)= RVINT*D31
  185.       AG(3,2,L)=-RVINT*D32
  186.       AG(3,3,L)= RVINT*D33
  187.  
  188.   310 CONTINUE
  189. C     write(6,*)' AJ ds CALJBR'
  190. C     write(6,1002)AJ
  191.       DO 331 L=1,NPG
  192.       DO 331 I=1,NP
  193.       DO 331 N=1,ND
  194.       U=0.D0
  195.       DO 332 M=1,NES
  196.   332 U=U+AG(M,N,L)*GR(M,I,L)
  197.   331 HR(N,I,L)=U
  198. C     write(6,*)'GR'
  199. C     write(6,1002)gr
  200. C     write(6,*)'HR'
  201. C     write(6,1002)hr
  202.  
  203.       ELSEIF(NES.EQ.2.AND.ND.EQ.3)THEN
  204. C     write(6,*)' Je passe ICI'
  205.  
  206.       AIRE=0.D0
  207.       DO 210 L=1,NPG
  208.       AJX=0.D0
  209.       AJY=0.D0
  210.       AJZ=0.D0
  211.       BJX=0.D0
  212.       BJY=0.D0
  213.       BJZ=0.D0
  214.       DO 211 I=1,NP
  215.       AJX=AJX+GR(1,I,L)*XYZ(1,I)
  216.       AJY=AJY+GR(1,I,L)*XYZ(2,I)
  217.       AJZ=AJZ+GR(1,I,L)*XYZ(3,I)
  218.       BJX=BJX+GR(2,I,L)*XYZ(1,I)
  219.       BJY=BJY+GR(2,I,L)*XYZ(2,I)
  220.       BJZ=BJZ+GR(2,I,L)*XYZ(3,I)
  221.   211 CONTINUE
  222.  
  223.       XB=AJY*BJZ-AJZ*BJY
  224.       YB=AJZ*BJX-AJX*BJZ
  225.       ZB=AJX*BJY-AJY*BJX
  226.  
  227.       AJN=(XB*XB+YB*YB+ZB*ZB)**0.5D0
  228.  
  229.       PGSQ(L)=PG(L)*AJN
  230.       AIRE=AIRE+PGSQ(L)
  231.  
  232.       AJ(1,1,L)=AJX/AJN
  233.       AJ(2,1,L)=AJY/AJN
  234.       AJ(3,1,L)=AJZ/AJN
  235.       AJ(1,2,L)=BJX/AJN
  236.       AJ(2,2,L)=BJY/AJN
  237.       AJ(3,2,L)=BJZ/AJN
  238.       AJ(1,3,L)=XB/AJN
  239.       AJ(2,3,L)=YB/AJN
  240.       AJ(3,3,L)=ZB/AJN
  241.  
  242.       AG(1,1,L)=AJX
  243.       AG(2,1,L)=AJY
  244.       AG(3,1,L)=AJZ
  245.       AG(1,2,L)=BJX
  246.       AG(2,2,L)=BJY
  247.       AG(3,2,L)=BJZ
  248.       AG(1,3,L)=XB
  249.       AG(2,3,L)=YB
  250.       AG(3,3,L)=ZB
  251.  
  252.       D11=AG(2,2,L)*AG(3,3,L)-AG(3,2,L)*AG(2,3,L)
  253.       D12=AG(1,2,L)*AG(3,3,L)-AG(3,2,L)*AG(1,3,L)
  254.       D13=AG(1,2,L)*AG(2,3,L)-AG(2,2,L)*AG(1,3,L)
  255.       D21=AG(2,1,L)*AG(3,3,L)-AG(3,1,L)*AG(2,3,L)
  256.       D22=AG(1,1,L)*AG(3,3,L)-AG(3,1,L)*AG(1,3,L)
  257.       D23=AG(1,1,L)*AG(2,3,L)-AG(2,1,L)*AG(1,3,L)
  258.       D31=AG(2,1,L)*AG(3,2,L)-AG(3,1,L)*AG(2,2,L)
  259.       D32=AG(1,1,L)*AG(3,2,L)-AG(3,1,L)*AG(1,2,L)
  260.       D33=AG(1,1,L)*AG(2,2,L)-AG(2,1,L)*AG(1,2,L)
  261.       VINT=AG(1,1,L)*D11-AG(1,2,L)*D21+AG(1,3,L)*D31
  262.       RVINT=1.D0/VINT
  263.       AG(1,1,L)= RVINT*D11
  264.       AG(1,2,L)=-RVINT*D12
  265.       AG(1,3,L)= RVINT*D13
  266.       AG(2,1,L)=-RVINT*D21
  267.       AG(2,2,L)= RVINT*D22
  268.       AG(2,3,L)=-RVINT*D23
  269.       AG(3,1,L)= RVINT*D31
  270.       AG(3,2,L)=-RVINT*D32
  271.       AG(3,3,L)= RVINT*D33
  272.  
  273.  210  CONTINUE
  274. C
  275.       DO 231 L=1,NPG
  276.       DO 231 I=1,NP
  277.       DO 231 N=1,ND
  278.       U=0.D0
  279.       DO 232 M=1,NES
  280.   232 U=U+AG(M,N,L)*GR(M,I,L)
  281.   231 HR(N,I,L)=U
  282. C     write(6,*)'GR'
  283. C     write(6,1002)gr
  284. C     write(6,*)'HR'
  285. C     write(6,1002)hr
  286.  
  287.       ELSE
  288.  
  289.       DO 10 L=1,NPG
  290.       DO 10 M=1,ND
  291.       DO 10 N=1,ND
  292.       AJT=0.D0
  293.       DO 11 I=1,NP
  294.       AJT=AJT+GR(M,I,L)*XYZ(N,I)
  295.    11 CONTINUE
  296.       AJ(N,M,L)=AJT
  297.    10 CONTINUE
  298.  
  299. C
  300.       DO 20 L=1,NPG
  301.       IF(ND.EQ.1)THEN
  302.       VINT=AJ(1,1,L)
  303. C     VINT=ABS(VINT)
  304.       AJ(1,1,L)=1.D0/VINT
  305.       ELSEIF(ND.EQ.2)THEN
  306.       VINT=AJ(1,1,L)*AJ(2,2,L)-AJ(1,2,L)*AJ(2,1,L)
  307. C     VINT=ABS(VINT)
  308.       RVINT=1.D0/VINT
  309.       D11=AJ(2,2,L)
  310.       D12=AJ(1,2,L)
  311.       D21=AJ(2,1,L)
  312.       D22=AJ(1,1,L)
  313.       AJ(1,1,L)= RVINT*D11
  314.       AJ(1,2,L)=-RVINT*D12
  315.       AJ(2,1,L)=-RVINT*D21
  316.       AJ(2,2,L)= RVINT*D22
  317.       ELSEIF(ND.EQ.3)THEN
  318.       D11=AJ(2,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(2,3,L)
  319.       D12=AJ(1,2,L)*AJ(3,3,L)-AJ(3,2,L)*AJ(1,3,L)
  320.       D13=AJ(1,2,L)*AJ(2,3,L)-AJ(2,2,L)*AJ(1,3,L)
  321.       D21=AJ(2,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(2,3,L)
  322.       D22=AJ(1,1,L)*AJ(3,3,L)-AJ(3,1,L)*AJ(1,3,L)
  323.       D23=AJ(1,1,L)*AJ(2,3,L)-AJ(2,1,L)*AJ(1,3,L)
  324.       D31=AJ(2,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(2,2,L)
  325.       D32=AJ(1,1,L)*AJ(3,2,L)-AJ(3,1,L)*AJ(1,2,L)
  326.       D33=AJ(1,1,L)*AJ(2,2,L)-AJ(2,1,L)*AJ(1,2,L)
  327.       VINT=AJ(1,1,L)*D11-AJ(1,2,L)*D21+AJ(1,3,L)*D31
  328. C     VINT=ABS(VINT)
  329.       RVINT=1.D0/VINT
  330.       AJ(1,1,L)= RVINT*D11
  331.       AJ(1,2,L)=-RVINT*D12
  332.       AJ(1,3,L)= RVINT*D13
  333.       AJ(2,1,L)=-RVINT*D21
  334.       AJ(2,2,L)= RVINT*D22
  335.       AJ(2,3,L)=-RVINT*D23
  336.       AJ(3,1,L)= RVINT*D31
  337.       AJ(3,2,L)=-RVINT*D32
  338.       AJ(3,3,L)= RVINT*D33
  339.       ELSE
  340.       VINT=0.D0
  341.       ENDIF
  342.  
  343.       PGSQ(L)=VINT
  344. C***  CALL DLAIN(ND,ND,AJ(1,1,L),KAUX,B,IER)
  345. C
  346.    20 CONTINUE
  347. C
  348.       DO 31 L=1,NPG
  349.       DO 31 I=1,NP
  350.       DO 31 N=1,ND
  351.       U=0.D0
  352.       DO 32 M=1,ND
  353.    32 U=U+AJ(M,N,L)*GR(M,I,L)
  354.    31 HR(N,I,L)=U
  355.       U=0.D0
  356.       V=0.D0
  357.       DO 4 L=1,NPG
  358.       W=ABS(PGSQ(L))
  359.       U=U+PG(L)*W
  360.       V=V+PG(L)*PGSQ(L)
  361.     4 PGSQ(L)=PG(L)*W
  362.       AIRE=U
  363.       SGN=SIGN(1.D0,V)
  364.  
  365.       ENDIF
  366.  
  367.       IF(IAXI.EQ.0)RETURN
  368. C
  369.       DO 45 L=1,NPG
  370.       RPGT  =0.D0
  371.       DO 46 I=1,NP
  372.       RPGT=RPGT+XYZ(1,I)*FN(I,L)
  373.  46   CONTINUE
  374.       RPG(L)=RPGT
  375.  45   CONTINUE
  376. C
  377.       DEUPI=2.D0*XPI
  378.       U=0.D0
  379.       DO 47 L=1,NPG
  380. C?    PGSQ(L)=PGSQ(L)*DEUPI*RPG(L)
  381.       U=U+PGSQ(L)*DEUPI*RPG(L)
  382.  47   CONTINUE
  383.       AIRE=U
  384. C
  385.       RETURN
  386.  1002 FORMAT(10(1X,1PE11.4))
  387.  1001 FORMAT(20(1X,I5))
  388.       END
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales