Télécharger la plaquette Cast3M


fr en

Télécharger piocaf.eso

Retour à la liste

Numérotation des lignes : 

piocaf

  1. C PIOCAF    SOURCE    PV090527  24/03/14    21:15:02     11870          
  2.       SUBROUTINE PIOCAF(NBNN,nbsh,IDIM,TAB1,NCOELE,NBPTEL,SHP,XE1,
  3.      1                   XE2,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   NBSH             = NOMBRE DE fonctions de forme
  13. C   IDIM             = DIMENSION DE L ESPACE SUPPORT
  14. C
  15. C   TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF
  16. C
  17. C   NCOELE           = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES
  18. C
  19. C   NBPTEL           = NOMBRE DE POINTS DE GAUSS
  20. C   SHP(6,NBNN,NBPTEL)= FONCTIONS DE FORME
  21. C
  22. C   KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS
  23. C
  24. C TABLEAUX DE TRAVAIL
  25. C--------------------
  26. C   XE1(3,NBNN)    = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART
  27. C
  28. C   XE2(3,NBNN)    = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL
  29. C
  30. C   SH1(6,NBNN)    = FONCTIONS DE FORME EN UN POINT DE GAUSS
  31. C
  32. C  SORTIES
  33. C---------
  34. C   TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY
  35. C
  36. C
  37. C     AOUT 85
  38. C     MODIF PEGON FEV 90  CAS BIDIM
  39. C     PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91
  40. C
  41. C=======================================================================
  42.       IMPLICIT INTEGER(I-N)
  43.       IMPLICIT REAL*8(A-H,O-Z)
  44. -INC CCREEL
  45. *
  46.       SEGMENT MWRK6
  47.        INTEGER ITRES1(NBPTEL)
  48.        REAL*8 PRODDI(NBPTEL,LHOO2),PRODDO(NBPTEL,LHOO2)
  49.        REAL*8 DDHOOK(LHOOK,LHOOK),DDHOMU(LHOOK,LHOOK)
  50.        REAL*8 VEC(LHOOK),VEC2(LHOOK)
  51.       ENDSEGMENT
  52. *
  53. C
  54.       dimension r_z(3)
  55.       DIMENSION TAB1(NBPTEL,*)
  56.       DIMENSION TAB(NBPTEL,*)
  57. *as xfem 2010_01_13
  58.       DIMENSION SHP(6,NBSH,*)
  59. *      DIMENSION SHP(6,NBNN,*)
  60. *fin as xfem 2010_01_13
  61.       DIMENSION XE1(3,*),XE2(3,*),SH1(6,*)
  62. C
  63. C     TABLEAUX DE TRAVAIL DIMENSIONNES ICI
  64. C
  65.       DIMENSION XJAC(3,3),FAC(6),FACI(6)
  66. C
  67. C     TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO
  68. C     DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS
  69. C
  70.       DIMENSION IN(6),JN(6),ITAB(3,3)
  71. C
  72.       DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/
  73.       DATA FACI/1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/
  74.       DATA IN/1,2,3,1,1,2/
  75.       DATA JN/1,2,3,2,3,3/
  76. C
  77.       DATA ITAB(1,1),ITAB(2,1),ITAB(3,1)/1,4,5/
  78.       DATA ITAB(1,2),ITAB(2,2),ITAB(3,2)/4,2,6/
  79.       DATA ITAB(1,3),ITAB(2,3),ITAB(3,3)/5,6,3/
  80. C
  81.       KERRE=0
  82. C
  83. C     MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS
  84. C
  85. **    DO 50 IA=1,NBPTEL
  86. **    DO 50 IB=1,NCOELE
  87. **    TAB(IA,IB)=0.D0
  88. **50  CONTINUE
  89.       CALL ZERO(TAB,NBPTEL,NCOELE)
  90. C
  91. C     BOUCLE SUR LES POINTS DE GAUSS
  92. C
  93.       DO 130 IC=1,NBPTEL
  94. C
  95.         if (nbsh.ne.nbnn) then
  96.           CALL HPRIMEX(XE1,NBNN,IDIM,NBSH,SHP,IC,SH1,DJAC)
  97.         else
  98.           CALL HPRIME(XE1,NBNN,IDIM,SHP,IC,SH1,DJAC)
  99.         endif
  100. C
  101. C       CALCUL DE LA MATRICE  F (=XJAC)
  102. C
  103.         CALL ZERO(XJAC,3,3)
  104.         DO  IF=1,IDIM
  105.           JF = IF + 1
  106.           r_z1 = 0.D0
  107.           r_z2 = 0.D0
  108.           r_z3 = 0.D0
  109.           do  id=1,nbnn
  110.                              r_z1    = r_z1    + SH1(JF,ID)*XE2(1,ID)
  111.               if (idim.ge.2) r_z2    = r_z2    + SH1(JF,ID)*XE2(2,ID)
  112.               if (idim.ge.3) r_z3    = r_z3    + SH1(JF,ID)*XE2(3,ID)
  113.           enddo
  114.           if (idim.ge.1) XJAC(1,IF) = r_z1
  115.           if (idim.ge.2) XJAC(2,IF) = r_z2
  116.           if (idim.ge.3) XJAC(3,IF) = r_z3
  117.         enddo
  118. C
  119. CC - MODES DE CALCUL EN DEFORMATIONS "PLANES" GENERALISEES
  120. CC CAS 3D :
  121.         IF (IDIM.EQ.3) THEN
  122. C         RIEN FAIRE !
  123. CC CAS 2D :
  124.         ELSE IF(IDIM.EQ.2) THEN
  125. CCCC CAS 2D AXISYMETRIQUE
  126.            IF (IFOU.EQ.0) THEN
  127.              R1=0.D0
  128.              R2=0.D0
  129.              DO 150 ID=1,NBNN
  130.                R1=R1+SH1(1,ID)*XE1(1,ID)
  131.                R2=R2+SH1(1,ID)*XE2(1,ID)
  132.  150         CONTINUE
  133.              if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  134.                XJAC(3,3)=R2/R1
  135.              else
  136.                xjac(3,3)=xgrand*xzprec
  137.              endif
  138.            ELSE
  139. CCCC CAS 2D PLAN
  140.              XJAC(3,3)=1.D0
  141. CCCC CAS 2D PLAN DEFO GENE
  142. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  143.              IF (IFOU.EQ.-3) THEN
  144.                XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  145.              ENDIF
  146.            ENDIF
  147. CC CAS 1D :
  148.         ELSE IF (IDIM.EQ.1) THEN
  149. CCCC CAS 1D PLAN
  150.            IF (IFOU.GE.3 .AND. IFOU.LE.11) THEN
  151.              XJAC(2,2) = 1.D0
  152.              XJAC(3,3) = 1.D0
  153. * Rq : "Deplacement" UY du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  154.              IF (IFOU.EQ.7 .OR. IFOU.EQ.8 .OR. IFOU.EQ.11) THEN
  155.                XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  156.              ENDIF
  157. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP)
  158. c*           IF (IFOU.EQ.9 .OR. IFOU.EQ.10 .OR. IFOU.EQ.11) THEN
  159.              IF (IFOU.GE.9) THEN
  160.                XJAC(3,3) = XJAC(3,3) + XE2(3,1)
  161.              ENDIF
  162. CCCC CAS 1D AXIS et SPHE
  163.            ELSE
  164.              CALL DISTRR(XE1,SH1,NBNN,R1)
  165.              CALL DISTRR(XE2,SH1,NBNN,R2)
  166.              if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  167.                 FR2R1 = R2 /R1
  168.              else
  169.                 FR2R1=XGRAND*XZPREC
  170.              endif
  171.              XJAC(3,3) = FR2R1
  172. CCCC CAS 1D SPHE
  173.              IF (IFOU.EQ.15) THEN
  174.                XJAC(2,2) = FR2R1
  175.              ELSE
  176. CCCC CAS 1D AXIS
  177.                XJAC(2,2) = 1.D0
  178. * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP)
  179.                IF (IFOU.EQ.14) THEN
  180.                  XJAC(2,2) = XJAC(2,2) + XE2(2,1)
  181.                ENDIF
  182.              ENDIF
  183.            ENDIF
  184.         ENDIF
  185. C
  186.         GO TO (500,600,700),KCAS
  187. C
  188. C         KCAS=1   CAS DES CONTRAINTES
  189. C         ----------------------------
  190. C
  191.   500   CONTINUE
  192. C
  193. CCCCCCCCCCCC  CALCUL DE L'INVERSE DU DETERMINANT DE F
  194. C
  195.         IF (IDIM.EQ.3) THEN
  196.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  197.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  198.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  199.         ELSE IF (IDIM.EQ.2) THEN
  200.           DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  201.         ELSE IF (IDIM.EQ.1) THEN
  202.           DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  203.         ENDIF
  204.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  205.           DETF=1./(DETF)
  206.         else
  207.           DETF=XGRAND*xzprec
  208.         endif
  209. C
  210. C     CALCUL DES CONTRAINTES DE CAUCHY
  211. C
  212.         DO 160 ID=1,NCOELE
  213.          IND=IN(ID)
  214.          JND=JN(ID)
  215.          TABCD = TAB(IC,ID)
  216.          DO 170 IE=1,IDIM
  217.          TABAU1=0.D0
  218.          DO 171 IF=1,IDIM
  219.            ICO=ITAB(IF,IE)
  220.            TABAU1=TAB1(IC,ICO)*XJAC(JND,IF)
  221.      1             +TABAU1
  222.  171     CONTINUE
  223.          TABAU1 = TABAU1 * XJAC(IND,IE)
  224.          TABCD = TABCD + TABAU1
  225.  170     CONTINUE
  226.          TAB(IC,ID)=TABCD*DETF  + TAB(IC,ID)
  227.  160    CONTINUE
  228. C
  229. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  230. C
  231.        IF (IDIM.EQ.2) THEN
  232.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  233.        ELSE IF (IDIM.EQ.1) THEN
  234.          TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  235.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  236.        ENDIF
  237.        GO TO 130
  238. C
  239. C         KCAS=2   CAS DES DEFORMATIONS
  240. C         -----------------------------
  241. C
  242.   600   CONTINUE
  243. C
  244. C
  245. CCCCCCCCCCCC  CALCUL DE L'INVERSE DE F
  246. C
  247.           CALL INVMA1(XJAC,3,3,KERRE)
  248.           IF(KERRE.NE.0) THEN
  249.              WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  250. 77881        FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  251.              RETURN
  252.           ENDIF
  253. C
  254. C     CALCUL DES DEFORMATIONS
  255. C
  256.         DO 260 ID=1,NCOELE
  257.          IND=IN(ID)
  258.          JND=JN(ID)
  259.          TABCD=0.D0
  260.          DO 271 IE=1,IDIM
  261. *  unrolling de la boucle 270 pour mettre xjac en registre
  262. **       TABAU1=0.D0
  263. **       DO 270 IF=1,IDIM
  264. **         ICO=ITAB(IF,IE)
  265. **         TABAU1=TABAU1+
  266. **   1      FAC(ICO)*TAB1(IC,ICO)*XJAC(IF,JND)
  267. *270     CONTINUE
  268.          ICO1=ITAB(1,IE)
  269.          tabau1=          FAC(ICO1)*TAB1(IC,ICO1)*XJAC(1,JND)
  270.          ICO2=ITAB(2,IE)
  271.          tabau1= tabau1 + FAC(ICO2)*TAB1(IC,ICO2)*XJAC(2,JND)
  272.          if(idim.eq.3) then
  273.          ICO3=ITAB(3,IE)
  274.          tabau1= tabau1 + FAC(ICO3)*TAB1(IC,ICO3)*XJAC(3,JND)
  275.          endif
  276.          TABCD=TABCD+TABAU1*XJAC(IE,IND)
  277.  271     CONTINUE
  278.          TAB(IC,ID)=TAB(IC,ID)+TABCD*FACI(ID)
  279.  260    CONTINUE
  280. C
  281. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  282. C
  283.        IF(IDIM.EQ.2) THEN
  284.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  285.        ELSE IF(IDIM.EQ.1) THEN
  286.          TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF
  287.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  288.        ENDIF
  289. C
  290.        GO TO 130
  291. C
  292. C         KCAS=3   CAS DE LA MATRICE DE HOOKE 
  293. C         -----------------------------------
  294. C
  295.   700   CONTINUE
  296. C
  297. CCCCCCCCCCCC  CALCUL DE L'INVERSE DU DETERMINANT DE F
  298. C
  299.         IF (IDIM.EQ.3) THEN
  300.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  301.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  302.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  303.         ELSE IF (IDIM.EQ.2) THEN
  304.           DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  305.         ELSE IF (IDIM.EQ.1) THEN
  306.           DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  307.         ENDIF
  308.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  309.           DETF=1./(DETF)
  310.         else
  311.           DETF=XGRAND*xzprec
  312.         endif
  313. C
  314.       IJ=1
  315.       DO  JJ=1,LHOOK
  316.         DO  II=1,LHOOK
  317.           DDHOOK(II,JJ)=PRODDI(IC,IJ)
  318.           IJ=IJ+1
  319.         enddo
  320.       enddo
  321. *
  322.         CALL ZERO(DDHOMU,LHOOK,LHOOK)
  323.         KEY=2
  324. C
  325.        DO 720 LC=1,LHOOK
  326.  
  327.         CALL ZERO (VEC,LHOOK,1)
  328.         DO 760 ID=1,LHOOK
  329.          IND=IN(ID)
  330.          JND=JN(ID)
  331.          DO  IE=1,IDIM
  332.          DO  IF=1,IDIM
  333.            ICO=ITAB(IF,IE)
  334.            IF(ICO.EQ.LC) THEN
  335.              VEC(ID)=VEC(ID)+
  336.      &         FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)
  337.            ENDIF
  338.          enddo
  339.          enddo
  340.          VEC(ID)=VEC(ID)*FACI(ID)
  341.  760    CONTINUE
  342. C
  343. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  344. C
  345.        IF (IDIM.EQ.2) THEN
  346.          IF(LC.EQ.3)  VEC(3)=XJAC(3,3)*XJAC(3,3)
  347. *
  348.        ELSE IF (IDIM.EQ.1) THEN
  349.          IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  350.          IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  351.        ENDIF
  352. C
  353.        CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  354. C
  355.         DO 761 ID=1,LHOOK
  356.          IND=IN(ID)
  357.          JND=JN(ID)
  358.          DO  IE=1,IDIM
  359.          DO  IF=1,IDIM
  360.            ICO=ITAB(IF,IE)
  361.            DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  362.      1             +DDHOMU(ID,LC)
  363.          enddo
  364.          enddo
  365.  761    CONTINUE
  366. C
  367. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  368. C
  369.        IF (IDIM.EQ.2) THEN
  370.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  371.        ELSE IF (IDIM.EQ.1) THEN
  372.          DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  373.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  374.        ENDIF
  375.  
  376.   720 CONTINUE
  377. C
  378.       IJ=1
  379.       DO  JJ=1,LHOOK
  380.         DO  II=1,LHOOK
  381.           PRODDO(IC,IJ)=DDHOMU(II,JJ) 
  382.           IJ=IJ+1
  383.         enddo
  384.       enddo
  385. *
  386. *
  387.  130  CONTINUE
  388.       RETURN
  389.       END
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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