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

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