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    PV        18/11/20    21:15:04     1001           
  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 50 IB=1,NCOELE
  103.         DO 50 IA=1,NBPTEL
  104.           TAB(IA,IB)=0.D0
  105.   50  CONTINUE
  106. C
  107. C     BOUCLE SUR LES POINTS DE GAUSS
  108. C
  109.       DO 130 IC=1,NBPTEL
  110.  
  111.         tabphg(2)=0.D0
  112. *as xfem 2010_01_13
  113. *    Initialisation de SH1(Ni, Ni,x, Ni,y)
  114.         do i3=1,nbnn
  115.           xh1(i3)=0.D0
  116.           do i4=1,IDIM1
  117.             SH1(i4,i3)=SHPTOT(i4,i3,IC)
  118.           enddo
  119.         enddo
  120. *      Calcul de H et phi aux noeuds et au point de Gauss IC
  121.         do 131 i1=1,nbnn
  122.           mlree1=tlreel(i1)
  123.           if(mlree1.eq.0) goto 131
  124.           tabphg(1)=0.D0
  125.           do i2=1,nbnn
  126.             xphi1(i2)=mlree1.PROG(i2)
  127.             if (abs(xphi1(i2)).lt.1.d-7) then
  128.                xh1(i2)=0.D0
  129.             else
  130.                xh1(i2)=sign(1.d0,xphi1(i2))
  131.             endif
  132.             tabphg(1)=tabphg(1)+SH1(1,i2)*xphi1(i2)
  133.           enddo
  134.            if (abs(tabphg(1)).lt.1.d-7) then
  135.               tabphg(2)=0.D0
  136.             else
  137.               tabphg(2)=sign(1.d0,tabphg(1))
  138.             endif
  139.  
  140.  131    continue
  141. *    Calcul des H(x)-H(xi) :
  142.         do i3=1,nbnn
  143.          tabdh(i3)=tabphg(2)-xh1(i3)
  144.         enddo
  145. *   Calcul de SH1 :
  146.        call jacobix(XE1,TABA1,TABDH,SH1,IDIM,NBNN,DJac)
  147. C
  148. C       CALCUL DE LA MATRICE  F
  149. C
  150.         CALL ZERO(XJAC,3,3)
  151.         DO 140 ID=1,NBNN
  152.         DO 140 IE=1,IDIM
  153. *         r_z = XE2(IE,ID)
  154.           r_z = XE2(IE,ID)+(tabdh(ID)*TABA2(IE,ID))
  155.         DO 140 IF=1,IDIM
  156.           XJAC(IE,IF)=XJAC(IE,IF) + SH1(IF+1,ID)*r_z
  157.   140   CONTINUE
  158. *fin as xfem 2010_01_13
  159.         IF(IDIM.EQ.2) THEN
  160.            XJAC(3,3)=1.D0
  161.            IF(IFOU.EQ.0) THEN
  162. C
  163. CCCCCCCCCCCCC CAS AXISYMETRIQUE
  164. C
  165.              R1=0.
  166.              R2=0.
  167.              DO 150 ID=1,NBNN
  168.                R1=R1+SH1(1,ID)*XE1(1,ID)
  169.                R2=R2+SH1(1,ID)*XE2(1,ID)
  170.  150         CONTINUE
  171.              if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then
  172.                XJAC(3,3)=R2/R1
  173.              else
  174.                xjac(3,3)=xgrand*xzprec
  175.              endif
  176.            ENDIF
  177.         ENDIF
  178. C
  179. C
  180.         GO TO (500,600,700),KCAS
  181. C
  182. C         KCAS=1   CAS DES CONTRAINTES
  183. C         ----------------------------
  184. C
  185.   500   CONTINUE
  186. C
  187. C
  188. CCCCCCCCCCCC  CALCUL DE DETERMINANT DE F
  189. C
  190.         IF(IDIM.EQ.2) THEN
  191.           DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1)
  192.           DETF = DETF * XJAC (3,3)
  193.         ENDIF
  194.         IF(IDIM.EQ.3) THEN
  195.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  196.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  197.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  198.         ENDIF
  199.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  200.           DETF=1./(DETF)
  201.         else
  202.           DETF=XGRAND*xzprec
  203.         endif
  204. C
  205. C     CALCUL DES CONTRAINTES DE CAUCHY
  206. C
  207.         DO 160 ID=1,NCOELE
  208.          IND=IN(ID)
  209.          JND=JN(ID)
  210.          DO 170 IE=1,IDIM
  211.          DO 170 IF=1,IDIM
  212.            ICO=ITAB(IE,IF)
  213.            TAB(IC,ID)=TAB1(IC,ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  214.      1             +TAB(IC,ID)
  215.  170     CONTINUE
  216.  160    CONTINUE
  217. C
  218. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  219. C
  220.        IF(IDIM.EQ.2) THEN
  221.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF
  222.        ENDIF
  223.        GO TO 130
  224. C
  225. C         KCAS=2   CAS DES DEFORMATIONS
  226. C         -----------------------------
  227. C
  228.   600   CONTINUE
  229. C
  230. C
  231. CCCCCCCCCCCC  CALCUL DE L'INVERSE DE F
  232. C
  233.           CALL INVMA1(XJAC,3,3,KERRE)
  234.           IF(KERRE.NE.0) THEN
  235.              WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3)
  236. 77881        FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/))
  237.              RETURN
  238.           ENDIF
  239. C
  240. C     CALCUL DES DEFORMATIONS
  241. C
  242.         DO 260 ID=1,NCOELE
  243.          IND=IN(ID)
  244.          JND=JN(ID)
  245.          DO 270 IE=1,IDIM
  246.          DO 270 IF=1,IDIM
  247.            ICO=ITAB(IE,IF)
  248.            TAB(IC,ID)=TAB(IC,ID) +
  249.      1      FAC(ICO)*TAB1(IC,ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  250.  270     CONTINUE
  251.  260    CONTINUE
  252. C
  253. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  254. C
  255.        IF(IDIM.EQ.2) THEN
  256.          TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)
  257.        ENDIF
  258. C
  259.        GO TO 130
  260. C
  261. C         KCAS=3   CAS DE LA MATRICE DE HOOKE 
  262. C         -----------------------------------
  263. C
  264.   700   CONTINUE
  265. C
  266. CCCCCCCCCCCC  CALCUL DE L'INVERSE DU DETERMINANT DE F
  267. C
  268.         IF (IDIM.EQ.3) THEN
  269.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  270.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  271.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  272.         ELSE IF (IDIM.EQ.2) THEN
  273.           DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3)
  274.         ELSE IF (IDIM.EQ.1) THEN
  275.           DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3)
  276.         ENDIF
  277.         if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then
  278.           DETF=1./(DETF)
  279.         else
  280.           DETF=XGRAND*xzprec
  281.         endif
  282. C
  283.       IJ=1
  284.       DO 710 JJ=1,LHOOK
  285.         DO 710 II=1,LHOOK
  286.           DDHOOK(II,JJ)=PRODDI(IC,IJ)
  287.           IJ=IJ+1
  288.   710 CONTINUE
  289. *
  290.         CALL ZERO(DDHOMU,LHOOK,LHOOK)
  291.         KEY=2
  292. C
  293.        DO 720 LC=1,LHOOK
  294.  
  295.         CALL ZERO (VEC,LHOOK,1)
  296.         DO 760 ID=1,LHOOK
  297.          IND=IN(ID)
  298.          JND=JN(ID)
  299.          DO 770 IE=1,IDIM
  300.          DO 770 IF=1,IDIM
  301.            ICO=ITAB(IE,IF)
  302.            IF(ICO.EQ.LC) THEN
  303.              VEC(ID)=VEC(ID)+
  304.      &         FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID)
  305.            ENDIF
  306.  770     CONTINUE
  307.  760    CONTINUE
  308. C
  309. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  310. C
  311.        IF (IDIM.EQ.2) THEN
  312.          IF(LC.EQ.3)  VEC(3)=XJAC(3,3)*XJAC(3,3)
  313. *
  314.        ELSE IF (IDIM.EQ.1) THEN
  315.          IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2)
  316.          IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3)
  317.        ENDIF
  318. C
  319.        CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY)
  320. C
  321.         DO 761 ID=1,LHOOK
  322.          IND=IN(ID)
  323.          JND=JN(ID)
  324.          DO 771 IE=1,IDIM
  325.          DO 771 IF=1,IDIM
  326.            ICO=ITAB(IE,IF)
  327.            DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF
  328.      1             +DDHOMU(ID,LC)
  329.  771     CONTINUE
  330.  761    CONTINUE
  331. C
  332. C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  333. C
  334.        IF (IDIM.EQ.2) THEN
  335.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  336.        ELSE IF (IDIM.EQ.1) THEN
  337.          DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF
  338.          DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF
  339.        ENDIF
  340.  
  341.   720 CONTINUE
  342. C
  343.       IJ=1
  344.       DO 780 JJ=1,LHOOK
  345.         DO 780 II=1,LHOOK
  346.           PRODDO(IC,IJ)=DDHOMU(II,JJ) 
  347.           IJ=IJ+1
  348.   780 CONTINUE
  349. *
  350. *
  351.  130  CONTINUE
  352.       RETURN
  353.       END
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  

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