Télécharger la plaquette Cast3M


fr en

Télécharger piocax.eso

Retour à la liste

Numérotation des lignes : 

piocax

  1. C PIOCAX    SOURCE    CB215821  26/08/24    21:17:40     12622          
  2.       SUBROUTINE PIOCAX(NBNN,IDIM,TAB1,NCOELE,NBPTEL,IPMINT,XE1,XE2,
  3.      1           TABA,MRACC,SH1,TAB,MWRK6,LHOOK,
  4.      2           IFOU,KCAS,KERRE)
  5. C=======================================================================
  6. C
  7. C    TRANSFORME LES CONTRAINTES DE PIOLA KIRCHHOFF EN CONTRAINTES DE
  8. C                             CAUCHY
  9. C  ENTREE
  10. C -------
  11. C   NBNN             = NOMBRE DE POINTS PAR ELEMENTS
  12. C   IDIM             = DIMENSION DE L ESPACE SUPPORT
  13. C
  14. C   TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF
  15. C
  16. C   NCOELE           = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES
  17. C
  18. C   NBPTEL           = NOMBRE DE POINTS DE GAUSS
  19. C   IPMINT           = POINTEUR DES FONCTIONS DE FORME
  20. C   TABA             = pointeur tableau avec ddl de saut
  21. c   MRACC            = pointeur tableau de raccourci pour les
  22. C                      enrichissements elementaires
  23. C
  24. C   KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS
  25. C
  26. C TABLEAUX DE TRAVAIL
  27. C--------------------
  28. C   XE1(3,NBNN)    = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART
  29. C
  30. C   XE2(3,NBNN)    = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL
  31. C
  32. C   SH1(6,NBNN)    = FONCTIONS DE FORME EN UN POINT DE GAUSS
  33. C
  34. C  SORTIES
  35. C---------
  36. C   TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY
  37. C
  38. C
  39. C     AOUT 85
  40. C     MODIF PEGON FEV 90  CAS BIDIM
  41. C     PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91
  42. C
  43. C=======================================================================
  44.       IMPLICIT INTEGER(I-N)
  45.       IMPLICIT REAL*8(A-H,O-Z)
  46. C
  47. -INC CCREEL
  48. -INC SMLREEL
  49. -INC SMINTE
  50. *
  51.       SEGMENT MWRK6
  52.        INTEGER ITRES1(NBPTEL)
  53.        REAL*8 PRODDI(NBPTEL,LHOO2),PRODDO(NBPTEL,LHOO2)
  54.        REAL*8 DDHOOK(LHOOK,LHOOK),DDHOMU(LHOOK,LHOOK)
  55.        REAL*8 VEC(LHOOK),VEC2(LHOOK)
  56.       ENDSEGMENT
  57. *
  58.  
  59.       DIMENSION TAB1(NBPTEL,*)
  60.       DIMENSION TAB(NBPTEL,*)
  61.       DIMENSION XE1(3,*),XE2(3,*),SH1(6,*)
  62. *as xfem 2010_01_13
  63.       PARAMETER (NBNNMAX=20)
  64.       SEGMENT MRACC
  65.         INTEGER TLREEL(NBNN)
  66.       ENDSEGMENT
  67. *     ddl de saut enrichissement
  68.       SEGMENT TABA
  69.         REAL*8 TABA1(IDIM,NBNN),TABA2(IDIM,NBNN)
  70.       ENDSEGMENT
  71.  
  72. *     phi et H aux nnoeuds ; phi et H au point de Gauss courant
  73.       DIMENSION xphi1(NBNNMAX),xh1(NBNNMAX),TABPHG(2),tabdh(nbnnMAX)
  74. *fin as xfem 2010_01_13
  75. C
  76. C     TABLEAUX DE TRAVAIL DIMENSIONNES ICI
  77. C
  78.       DIMENSION XJAC(3,3),FAC(6)
  79. C
  80. C     TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO
  81. C     DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS
  82. C
  83.       DIMENSION IN(6),JN(6),ITAB(3,3)
  84. C
  85.       DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/
  86.       DATA IN/1,2,3,1,1,2/
  87.       DATA JN/1,2,3,2,3,3/
  88. C
  89.       DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/
  90.       DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/
  91.       DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/
  92. C
  93.       KERRE=0
  94.  
  95.       MINTE = IPMINT
  96.       NBSH = SHPTOT(/2)
  97.  
  98.       IDIM1=IDIM+1
  99. C
  100. C     MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS
  101. C
  102.       DO 77882 IB=1,NCOELE
  103.         DO 50 IA=1,NBPTEL
  104.           TAB(IA,IB)=0.D0
  105.   50  CONTINUE
  106. 77882 CONTINUE
  107. C
  108. C     BOUCLE SUR LES POINTS DE GAUSS
  109. C
  110.       DO 130 IC=1,NBPTEL
  111.  
  112.         tabphg(2)=0.D0
  113. *as xfem 2010_01_13
  114. *    Initialisation de SH1(Ni, Ni,x, Ni,y)
  115.         do i3=1,nbnn
  116.           xh1(i3)=0.D0
  117.           do i4=1,IDIM1
  118.             SH1(i4,i3)=SHPTOT(i4,i3,IC)
  119.           enddo
  120.         enddo
  121. *      Calcul de H et phi aux noeuds et au point de Gauss IC
  122.         do 131 i1=1,nbnn
  123.           mlree1=tlreel(i1)
  124.           if(mlree1.eq.0) goto 131
  125.           tabphg(1)=0.D0
  126.           do i2=1,nbnn
  127.             xphi1(i2)=mlree1.PROG(i2)
  128.             if (abs(xphi1(i2)).lt.1.d-7) then
  129.                xh1(i2)=0.D0
  130.             else
  131.                xh1(i2)=sign(1.d0,xphi1(i2))
  132.             endif
  133.             tabphg(1)=tabphg(1)+SH1(1,i2)*xphi1(i2)
  134.           enddo
  135.            if (abs(tabphg(1)).lt.1.d-7) then
  136.               tabphg(2)=0.D0
  137.             else
  138.               tabphg(2)=sign(1.d0,tabphg(1))
  139.             endif
  140.  
  141.  131    continue
  142. *    Calcul des H(x)-H(xi) :
  143.         do i3=1,nbnn
  144.          tabdh(i3)=tabphg(2)-xh1(i3)
  145.         enddo
  146. *   Calcul de SH1 :
  147.        call jacobix(XE1,TABA1,TABDH,SH1,IDIM,NBNN,DJac)
  148. C
  149. C       CALCUL DE LA MATRICE  F
  150. C
  151.         CALL ZERO(XJAC,3,3)
  152.         DO 77884 ID=1,NBNN
  153.         DO 77883 IE=1,IDIM
  154. *         r_z = XE2(IE,ID)
  155.           r_z = XE2(IE,ID)+(tabdh(ID)*TABA2(IE,ID))
  156.         DO 140 IF=1,IDIM
  157.           XJAC(IE,IF)=XJAC(IE,IF) + SH1(IF+1,ID)*r_z
  158.   140   CONTINUE
  159. 77883   CONTINUE
  160. 77884   CONTINUE
  161. *fin as xfem 2010_01_13
  162.         IF(IDIM.EQ.2) THEN
  163.            XJAC(3,3)=1.D0
  164.            IF(IFOU.EQ.0) THEN
  165. C
  166. CCCCCCCCCCCCC CAS AXISYMETRIQUE
  167. C
  168.              R1=0.
  169.              R2=0.
  170.              DO 150 ID=1,NBNN
  171.                R1=R1+SH1(1,ID)*XE1(1,ID)
  172.                R2=R2+SH1(1,ID)*XE2(1,ID)
  173.  150         CONTINUE
  174.              if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  175.                XJAC(3,3)=R2/R1
  176.              else
  177.                xjac(3,3)=xgrand*xzprec
  178.              endif
  179.            ENDIF
  180.         ENDIF
  181. C
  182. C
  183.         GO TO (500,600,700),KCAS
  184. C
  185. C         KCAS=1   CAS DES CONTRAINTES
  186. C         ----------------------------
  187. C
  188.   500   CONTINUE
  189. C
  190. C
  191. CCCCCCCCCCCC  CALCUL DE DETERMINANT DE F
  192. C
  193.         IF(IDIM.EQ.2) THEN
  194.           DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1)
  195.           DETF = DETF * XJAC (3,3)
  196.         ENDIF
  197.         IF(IDIM.EQ.3) THEN
  198.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  199.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  200.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  201.         ENDIF
  202.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  203.           DETF=1./(DETF)
  204.         else
  205.           DETF=XGRAND*xzprec
  206.         endif
  207. C
  208. C     CALCUL DES CONTRAINTES DE CAUCHY
  209. C
  210.         DO 160 ID=1,NCOELE
  211.          IND=IN(ID)
  212.          JND=JN(ID)
  213.          DO 77885 IE=1,IDIM
  214.          DO 170 IF=1,IDIM
  215.            ICO=ITAB(IE,IF)
  216.            TAB(IC,ID)=TAB1(IC,ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  217.      1             +TAB(IC,ID)
  218.  170     CONTINUE
  219. 77885    CONTINUE
  220.  160    CONTINUE
  221. C
  222. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  223. C
  224.        IF(IDIM.EQ.2) THEN
  225.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  226.        ENDIF
  227.        GO TO 130
  228. C
  229. C         KCAS=2   CAS DES DEFORMATIONS
  230. C         -----------------------------
  231. C
  232.   600   CONTINUE
  233. C
  234. C
  235. CCCCCCCCCCCC  CALCUL DE L'INVERSE DE F
  236. C
  237.           CALL INVMA1(XJAC,3,3,KERRE)
  238.           IF(KERRE.NE.0) THEN
  239.              WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  240. 77881        FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  241.              RETURN
  242.           ENDIF
  243. C
  244. C     CALCUL DES DEFORMATIONS
  245. C
  246.         DO 260 ID=1,NCOELE
  247.          IND=IN(ID)
  248.          JND=JN(ID)
  249.          DO 77886 IE=1,IDIM
  250.          DO 270 IF=1,IDIM
  251.            ICO=ITAB(IE,IF)
  252.            TAB(IC,ID)=TAB(IC,ID) +
  253.      1      FAC(ICO)*TAB1(IC,ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  254.  270     CONTINUE
  255. 77886    CONTINUE
  256.  260    CONTINUE
  257. C
  258. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  259. C
  260.        IF(IDIM.EQ.2) THEN
  261.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  262.        ENDIF
  263. C
  264.        GO TO 130
  265. C
  266. C         KCAS=3   CAS DE LA MATRICE DE HOOKE 
  267. C         -----------------------------------
  268. C
  269.   700   CONTINUE
  270. C
  271. CCCCCCCCCCCC  CALCUL DE L'INVERSE DU DETERMINANT DE F
  272. C
  273.         IF (IDIM.EQ.3) THEN
  274.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  275.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  276.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  277.         ELSE IF (IDIM.EQ.2) THEN
  278.           DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  279.         ELSE IF (IDIM.EQ.1) THEN
  280.           DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  281.         ENDIF
  282.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  283.           DETF=1./(DETF)
  284.         else
  285.           DETF=XGRAND*xzprec
  286.         endif
  287. C
  288.       IJ=1
  289.       DO 77887 JJ=1,LHOOK
  290.         DO 710 II=1,LHOOK
  291.           DDHOOK(II,JJ)=PRODDI(IC,IJ)
  292.           IJ=IJ+1
  293.   710 CONTINUE
  294. 77887 CONTINUE
  295. *
  296.         CALL ZERO(DDHOMU,LHOOK,LHOOK)
  297.         KEY=2
  298. C
  299.        DO 720 LC=1,LHOOK
  300.  
  301.         CALL ZERO (VEC,LHOOK,1)
  302.         DO 760 ID=1,LHOOK
  303.          IND=IN(ID)
  304.          JND=JN(ID)
  305.          DO 77888 IE=1,IDIM
  306.          DO 770 IF=1,IDIM
  307.            ICO=ITAB(IE,IF)
  308.            IF(ICO.EQ.LC) THEN
  309.              VEC(ID)=VEC(ID)+
  310.      &         FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  311.            ENDIF
  312.  770     CONTINUE
  313. 77888    CONTINUE
  314.  760    CONTINUE
  315. C
  316. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  317. C
  318.        IF (IDIM.EQ.2) THEN
  319.          IF(LC.EQ.3)  VEC(3)=XJAC(3,3)*XJAC(3,3)
  320. *
  321.        ELSE IF (IDIM.EQ.1) THEN
  322.          IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  323.          IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  324.        ENDIF
  325. C
  326.        CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  327. C
  328.         DO 761 ID=1,LHOOK
  329.          IND=IN(ID)
  330.          JND=JN(ID)
  331.          DO 77889 IE=1,IDIM
  332.          DO 771 IF=1,IDIM
  333.            ICO=ITAB(IE,IF)
  334.            DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  335.      1             +DDHOMU(ID,LC)
  336.  771     CONTINUE
  337. 77889    CONTINUE
  338.  761    CONTINUE
  339. C
  340. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  341. C
  342.        IF (IDIM.EQ.2) THEN
  343.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  344.        ELSE IF (IDIM.EQ.1) THEN
  345.          DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  346.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  347.        ENDIF
  348.  
  349.   720 CONTINUE
  350. C
  351.       IJ=1
  352.       DO 77890 JJ=1,LHOOK
  353.         DO 780 II=1,LHOOK
  354.           PRODDO(IC,IJ)=DDHOMU(II,JJ) 
  355.           IJ=IJ+1
  356.   780 CONTINUE
  357. 77890 CONTINUE
  358. *
  359. *
  360.  130  CONTINUE
  361.       RETURN
  362.       END
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  

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