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

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