Télécharger la plaquette Cast3M


fr en

Télécharger epplum.eso

Retour à la liste

Numérotation des lignes : 

epplum

  1. C EPPLUM    SOURCE    GF238795  18/02/05    21:15:14     9726           
  2.       SUBROUTINE EPPLUM(TENS,PPLUS,IG,VAL1,VP1,QPLUS,Q,VP,
  3.      .                  P,S,TAMP)
  4. C======================================================================
  5. C
  6. C                     SOUS PROGRAMME DE CALCUL
  7. C                               POUR
  8. C                     TENSEUR DE DEFORMATION
  9. C
  10. C                            VERSION 1.0
  11. C                            -----------
  12. C
  13. C
  14. C                            CALCUL DE :
  15. C
  16. C                 1- Valeurs et vecteurs propres
  17. C                 2- Tenseur Q
  18. C                 3- Tenseur Q+
  19. C                 4- Tenseur d ordre 4  P+
  20. C
  21. C======================================================================
  22. C
  23. C     CREATION : F.CORMERY
  24. C                E.N.S.M.A  -  LMPM
  25. C                DEC 1992
  26. C
  27. C======================================================================
  28.       IMPLICIT INTEGER(I-N)
  29.       IMPLICIT REAL*8 (A-H,O-Z)
  30. C**********************************************************************
  31. C                      DIMENSIONS ET DATA
  32. C**********************************************************************
  33. C                                        N36       N72    N75
  34. CC      DIMENSION PPLUS(3,3,3,3),TENS(6),QPLUS(6,6),Q(6,6),VP(3)
  35. C                  N84   N90
  36. CC     *          ,P(3,3),S(6),
  37. C                                  N99
  38. CC     *           VP1(3),VAL1(3,3),TAMP(3,3)
  39. C
  40.       DIMENSION PPLUS(3,3,3,*),TENS(*),QPLUS(6,*),Q(6,*),VP(*)
  41.      *          ,P(3,*),S(*),
  42.      *           VP1(*),VAL1(3,*),TAMP(3,*)
  43.       DATA ZERO/0.D0/,UN/1.D0/,
  44.      *    PRECIS/1.D-08/,DPRECS/1.D-08/
  45.       INTEGER IM
  46. C----------------------------------------------------------------------
  47.       AMAX1(X,Y,Z,U,V,W)= MAX(X,Y,Z,U,V,W)
  48. C------------------
  49.       IM=0
  50.       MT=10
  51. C**********************************************************************
  52. C                          INITIALISATION
  53. C**********************************************************************
  54.       DO 5 J=1,6
  55.          DO 6 K=1,6
  56. *            P(J,K)=ZERO
  57.             Q(J,K)=ZERO
  58.             QPLUS(J,K)=ZERO
  59. 6        CONTINUE
  60. 5     CONTINUE
  61.       DO 55 J=1,3
  62.       DO 55 K=1,3
  63.  55   P(J,K)=ZERO
  64. C----------------------------------------------------------------------
  65.       TENS(4)=TENS(4)/2
  66.       TENS(5)=TENS(5)/2
  67.       TENS(6)=TENS(6)/2
  68. C**********************************************************************
  69. C                     NORMALISATION DU TENSEUR A
  70. C**********************************************************************
  71. C
  72. C----------------------------------------------------------------------
  73. C                  Trouver la valeur max de TENS(I)
  74. C----------------------------------------------------------------------
  75.       DO 3 I=1,6
  76.          S(I)=ABS(TENS(I))
  77. 3     CONTINUE
  78. C--------------
  79.       TMAX=AMAX1(S(1),S(2),S(3),S(4),S(5),S(6))
  80.       IF(TMAX.EQ.0.D0)TMAX=UN
  81. C----------------------------------------------------------------------
  82. C       Normaliser a un la composante de TENS(I) la plus grande
  83. C----------------------------------------------------------------------
  84.       DO 4 I=1,6
  85.          TENS(I)=TENS(I)/TMAX
  86.          IF(ABS(TENS(I)).LE.1E-15) TENS(I)=0.D0
  87. 4     CONTINUE
  88. C------------------------ cas axes principaux -------------------------
  89.       NN=0
  90.       DO 234 IV=4,6
  91.          IF(ABS(TENS(IV)).LE.1E-15) NN=NN+1
  92. 234   CONTINUE
  93.       IF(NN.EQ.3)THEN
  94.       VP(1)=TENS(1)
  95.       VP(2)=TENS(2)
  96.       VP(3)=TENS(3)
  97.       DO 235 I=1,3
  98.                P(I,1)=0
  99.                P(I,2)=0
  100.                P(I,3)=0
  101.                P(I,I)=1
  102. 235   CONTINUE
  103.       goto 98
  104.       ENDIF
  105. C***********************************************************************
  106. C                   CALCUL DES VALEURS  PROPRES
  107. C***********************************************************************
  108.       CALL VALPRP(TENS(1),TENS(2),TENS(3),TENS(6),TENS(4),TENS(5),
  109.      *            VP(1),VP(2),VP(3))
  110. C***********************************************************************
  111. C                   CALCUL DES VECTEURS PROPRES
  112. C***********************************************************************
  113.          IM=2
  114.          IF(ABS(VP(1)-VP(2)).LT.1E-08)THEN
  115.          VP(1)=VP(3)
  116.          VP(2)=VP(2)
  117.          VP(3)=VP(2)
  118.          IM=2
  119.          ENDIF
  120.          IF(ABS(VP(1)-VP(3)).LT.1E-08)THEN
  121.          VP(1)=VP(2)
  122.          VP(2)=VP(3)
  123.          IM=2
  124.          ENDIF
  125.          IF(ABS(VP(2)-VP(3)).LT.1E-08)THEN
  126.          IM=2
  127.          ENDIF
  128. C----------------------------------------------------------------------
  129.          IMM=0
  130. C----------------------------------------------------------------------
  131.        DO 10 I=1,IM
  132. C-------------------
  133.          SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2
  134.          SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2
  135.          SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2
  136.          SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5)
  137.          SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4)
  138.          SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5)
  139. C----------------------------------------------------------------------
  140. C       WRITE(10,*)'MINEURS :'
  141. C       WRITE(10,*)SDET1,SDET2,SDET3,SDET4,SDET5,SDET6
  142. C----------------------------------------------------------------------
  143.  
  144.              IF (ABS(SDET1).GT.DPRECS) THEN
  145.                P(I,1)=UN
  146.                P(I,2)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET1
  147.                P(I,3)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET1
  148.                GOTO 96
  149. C-------------------
  150.              ENDIF
  151.              IF (ABS(SDET2).GT.DPRECS) THEN
  152.                P(I,2)=UN
  153.                P(I,1)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET2
  154.                P(I,3)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET2
  155.                GOTO 96
  156. C-------------------
  157.              ENDIF
  158.              IF (ABS(SDET3).GT.DPRECS) THEN
  159.                P(I,3)=UN
  160.                P(I,1)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET3
  161.                P(I,2)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET3
  162.                GOTO 96
  163. C--------------------
  164.              ENDIF
  165.              IF (ABS(SDET4).GT.DPRECS) THEN
  166.                P(I,1)=UN
  167.            P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET4
  168.                P(I,3)=((TENS(4)*(TENS(1)-VP(I)))-TENS(5)*TENS(6))/SDET4
  169.                GOTO 96
  170. C--------------------
  171.              ENDIF
  172.              IF (ABS(SDET5).GT.DPRECS) THEN
  173.                P(I,1)=(-(tens(4)**2)+(tens(2)-vp(i)))/sdet5
  174.                P(I,2)=(-tens(6)*(tens(3)-vp(i))+tens(5)*tens(4))/sdet5
  175.                P(I,3)=1
  176.                GOTO 96
  177. C--------------------
  178.              ENDIF
  179.              IF (ABS(SDET6).GT.DPRECS) THEN
  180.                P(I,3)=UN
  181.                P(I,1)=((-TENS(5)*TENS(4))+(TENS(3)-vp(i))*TENS(6))/SDET6
  182.             P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET6
  183.              ENDIF
  184. C-----------------------------------------------------------------------
  185.             SSDET1=TENS(1)-VP(I)
  186.             SSDET2=TENS(2)-VP(I)
  187.             SSDET3=TENS(3)-VP(I)
  188. C-------------------CAS PARTICULIERS------------------------------------
  189.             IF (ABS(SSDET1).LE.PRECIS) THEN
  190.                                      P(I,1)=1
  191.                                      P(I,2)=0
  192.                                      P(I,3)=0
  193.             GOTO 96
  194.             ENDIF
  195.             IF (ABS(SSDET2).LE.PRECIS) THEN
  196.                                      P(I,1)=0
  197.                                      P(I,2)=1
  198.                                      P(I,3)=0
  199.             GOTO 96
  200.             ENDIF
  201.             IF (ABS(SSDET3).LE.PRECIS) THEN
  202.                                      P(I,1)=0
  203.                                      P(I,2)=0
  204.                                      P(I,3)=1
  205.             GOTO 96
  206.             ENDIF
  207.             IF (ABS(SSDET1).GT.PRECIS) THEN
  208.                                      P(I,1)=-(TENS(6)+TENS(5))/SSDET1
  209.                                      P(I,2)=1
  210.                                      P(I,3)=1
  211. C-------------------
  212.             GOTO 96
  213.             ENDIF
  214.             IF (ABS(SSDET2).GT.PRECIS) THEN
  215.                                      P(I,1)=1
  216.                                      P(I,2)=-(TENS(6)+TENS(4))/SSDET2
  217.                                      P(I,3)=1
  218.             GOTO 96
  219. C-------------------
  220.             ENDIF
  221.             IF (ABS(SSDET3).GT.PRECIS) THEN
  222.                                      P(I,1)=1
  223.                                      P(I,2)=1
  224.                                      P(I,3)=-(TENS(5)+TENS(4))/SSDET3
  225.            GOTO 96
  226.            ENDIF
  227. C-----------------------------------------------------------------------
  228.          WRITE(MT,*)'ERREUR DANS VPROP.FOR'
  229.       WRITE(MT,1010)
  230. 1010  FORMAT(1X,'TENSEUR A  SYM. D ORDRE 2 :',
  231.      *       /1X,'--------------------------'/)
  232.       WRITE(MT,1001)TENS(1),TENS(6),TENS(5)
  233. 1001  FORMAT(15X,'*  ',3e20.7,'  *')
  234.       WRITE(MT,1002)TENS(6),TENS(2),TENS(4)
  235. 1002  FORMAT(15X,'*  ',3e20.7,'  *')
  236.       WRITE(MT,1003)TENS(5),TENS(4),TENS(3)
  237. 1003  FORMAT(15X,'*  ',3e20.7,'  *'/)
  238.          STOP
  239. C-----------------------------------------------------------------------
  240. 96       CONTINUE
  241. C-----------------------------------------------------------------------
  242. 10    CONTINUE
  243. 98       CONTINUE
  244.                IF(IM.EQ.2)THEN
  245.                  P(3,1)=P(1,2)*P(2,3)-P(1,3)*P(2,2)
  246.                  P(3,2)=P(1,3)*P(2,1)-P(1,1)*P(2,3)
  247.                  P(3,3)=P(1,1)*P(2,2)-P(1,2)*P(2,1)
  248.                ENDIF
  249. C------------------------------------------------------------------------
  250. C       On verifie que la base formee est bien directe
  251. C------------------------------------------------------------------------
  252.                  DIR=P(3,1)*(P(1,2)*P(2,3)-P(1,3)*P(2,2))+
  253.      *               P(3,2)*(P(1,3)*P(2,1)-P(1,1)*P(2,3))+
  254.      *               P(3,3)*(P(1,1)*P(2,2)-P(1,2)*P(2,1))
  255.                  IF(DIR.LT.ZERO)THEN
  256.                    DO 12 J=1,3
  257.                      TAMP1=VP(2)
  258.                      VP(2)=VP(1)
  259.                      VP(1)=TAMP1
  260.                      TAMP(1,J)=P(1,J)
  261.                      P(1,J)=P(2,J)
  262.                      P(2,J)=TAMP(1,J)
  263. 12                 CONTINUE
  264.                  ENDIF
  265. C------------------------------------------------------------------------
  266. C              Normalisation des vecteurs propres
  267. C------------------------------------------------------------------------
  268.       DO 11 J=1,5
  269.              DO 11 I=1,3
  270.                IF(ABS(P(I,1)).LT.1.E-15)P(I,1)=0.D0
  271.                IF(ABS(P(I,2)).LT.1.E-15)P(I,2)=0.D0
  272.                IF(ABS(P(I,3)).LT.1.E-15)P(I,3)=0.D0
  273.                RAC=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  274.                P(I,1)=P(I,1)/RAC
  275.                P(I,2)=P(I,2)/RAC
  276.                P(I,3)=P(I,3)/RAC
  277. 11           CONTINUE
  278. C------------------------------------------------------------------------
  279. C            Factorisation par le facteur de normalisation
  280. C------------------------------------------------------------------------
  281. C               VP1(I)=VP(I)*TMAX
  282. C************************************************************************
  283. C                     Calcul de Q et Q+
  284. C************************************************************************
  285.          DO 223 I=1,3
  286.              DO 50 J=1,3
  287.                 DO 40 K=1,3
  288.                    Q(J,K)=Q(J,K)+P(I,J)*P(I,K)
  289.                    IF (VP(I).GT.PRECIS) THEN
  290.                                  QPLUS(J,K)=QPLUS(J,K)+P(I,J)*P(I,K)
  291.                    ENDIF
  292. 40              CONTINUE
  293. 50           CONTINUE
  294. C-------------------------------------------------------------------------
  295. C               Fin de la boucle sur les valeur propres
  296. C-------------------------------------------------------------------------
  297. 223     CONTINUE
  298. C-----------------------------------------------------------------------
  299. C         Verification si la base orthogonale est bien construite
  300. C-----------------------------------------------------------------------
  301.              DO 51 J=1,3
  302.              DO 51 K=1,3
  303.                    IF (ABS(Q(J,K)).LE.1E-15)Q(J,K)=ZERO
  304.                    IF (ABS(QPLUS(J,K)).LE.1E-15)QPLUS(J,K)=ZERO
  305. C-----------------------
  306.       IF (J.NE.K)THEN
  307.       IF(ABS(Q(J,K)).GT.1E-3)THEN
  308.                           WRITE(MT,*)'ERREUR DANS VPROP.FOR:'
  309.                           WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.'
  310. C----------------------
  311.       DO 118 I=1,6
  312.          TENS(I)=TENS(I)
  313. 118   CONTINUE
  314. C-----------------------
  315.       WRITE(MT,1110)
  316. 1110  FORMAT(/1X,'TENSEUR A  SYM. D ORDRE 2 :',
  317.      *       /1X,'--------------------------'/)
  318.       WRITE(MT,1001)TENS(1),TENS(6),TENS(5)
  319. 1101  FORMAT(15X,'*  ',3e20.7,'  *')
  320.       WRITE(MT,1002)TENS(6),TENS(2),TENS(4)
  321. 1102  FORMAT(15X,'*  ',3e20.7,'  *')
  322.       WRITE(MT,1003)TENS(5),TENS(4),TENS(3)
  323. 1103  FORMAT(15X,'*  ',3e20.7,'  *')
  324.       DO 31 I=1,3
  325.          SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2
  326.          SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2
  327.          SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2
  328.          SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5)
  329.          SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4)
  330.          SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5)
  331.       write(mt,*)sdet1,sdet2,sdet3,sdet4,sdet5,sdet6
  332.                WRITE(MT,5000)VP(I),(P(I,L),L=1,3)
  333. 5000   FORMAT(3X,'VALEUR PRO. :',E12.5,' VECTEUR PROP.:',3E12.5)
  334. 31    CONTINUE
  335. C----------------
  336.       WRITE(MT,*)'********** EPPLUSM********************'
  337.       WRITE(MT,*)'  TENSEURS Q :'
  338.       WRITE(MT,*)
  339. C----------------
  340.       DO 131 I=1,3
  341.       WRITE(MT,7000)(Q(I,L),L=1,3)
  342. 7000  FORMAT(5X,'*  ',3E12.5,'  *')
  343. 131   CONTINUE
  344. C----------------
  345.       WRITE(MT,*)
  346.       WRITE(MT,*)'  TENSEURS Q+ :'
  347.       WRITE(MT,*)
  348. C----------------
  349.       DO 132 I=1,3
  350.       WRITE(MT,7001)(QPLUS(I,L),L=1,3)
  351. 7001  FORMAT(5X,'*  ',3E12.5,'  *')
  352. 132   CONTINUE
  353. C----------------
  354.       STOP
  355. C----------------
  356.       ENDIF
  357.       ENDIF
  358. C-----------------------
  359.       IF (J.EQ.K)THEN
  360.       IF(ABS(Q(J,K)-UN).GT.1E-3)THEN
  361.                           WRITE(MT,*)'ERREUR DANS VPROP.FOR:'
  362.                           WRITE(MT,*)'- LA BASE N EST PAS ORTHOGONALE.'
  363.       ENDIF
  364.       ENDIF
  365. 51    CONTINUE
  366. 56    CONTINUE
  367. C**********************************************************************
  368. C                    CALCUL DE L OPERATEUR P+
  369. C**********************************************************************
  370.       DO 100 I=1,3
  371.          DO 100 J=1,3
  372.             DO 100 K=1,3
  373.                DO 100 L=1,3
  374.                PPLUS(I,J,K,L)=ZERO
  375.                do 114 M=1,3
  376.                do 114 N=1,3
  377.                         PPLUS(I,J,K,L)=PPLUS(I,J,K,L)+QPLUS(I,M)
  378.      *       *QPLUS(J,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N))/4.D0+
  379.      *        QPLUS(J,M)*QPLUS(I,N)*(Q(K,M)*Q(L,N)+Q(L,M)*Q(K,N))
  380.      *       /4.D0
  381. 114    continue
  382.               IF(ABS(PPLUS(I,J,K,L)).LE.1E-15)PPLUS(I,J,L,K)=0.D0
  383. 100    CONTINUE
  384. C-----------------------
  385. c      if (ig.eq.1)then
  386. c      WRITE(MT,1117)
  387. c1117  FORMAT(1X,'TENSEUR A  SYM. D ORDRE 2 :')
  388. c      WRITE(MT,1001)TENS(1),TENS(6),TENS(5)
  389. c      WRITE(MT,1002)TENS(6),TENS(2),TENS(4)
  390. c      WRITE(MT,1003)TENS(5),TENS(4),TENS(3)
  391. c      DO 234 I=1,3
  392. c               WRITE(MT,5000)VP(I),(P(I,L),L=1,3)
  393. c234    CONTINUE
  394. C----------------
  395. c      WRITE(MT,*)'  TENSEURS Q :'
  396. C----------------
  397. c      DO 231 I=1,3
  398. c      WRITE(MT,7000)(Q(I,L),L=1,3)
  399. c231   CONTINUE
  400. c      endif
  401. C-----------------------------------------------------------------------
  402. C             Multiplication par le facteur de normalisation
  403. C-----------------------------------------------------------------------
  404.       DO 110 I=1,6
  405.          TENS(I)=TENS(I)*TMAX
  406. 110   CONTINUE
  407. C-----------------------------------------------------------------------
  408.       TENS(4)=2*TENS(4)
  409.       TENS(5)=2*TENS(5)
  410.       TENS(6)=2*TENS(6)
  411. C----------------------------------------------------------------------
  412.       DO 341 I=1,3
  413.       DO 341 J=1,3
  414.       VAL1(I,J)=P(I,J)
  415. 341   CONTINUE
  416.       DO 342 I=1,3
  417.       VP1(I)=VP(I)*tmax
  418. 342   CONTINUE
  419. C-------------------------------------------------------------------------
  420.       RETURN
  421.       END
  422. C-----------------------------------------------------------------------
  423.  
  424.  
  425.  
  426.  
  427.  

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