Télécharger vect.eso

Retour à la liste

Numérotation des lignes :

  1. C VECT SOURCE CHAT 05/01/13 04:05:16 5004
  2. SUBROUTINE VECT(TENS,VAL,VP,P,S,TAMP)
  3. C======================================================================
  4. C
  5. C SOUS PROGRAMME DE CALCUL
  6. C POUR
  7. C TENSEUR DE DEFORMATION
  8. C
  9. C VERSION 1.0
  10. C -----------
  11. C
  12. C
  13. C CALCUL DE :
  14. C
  15. C 1- Valeurs et vecteurs propres
  16. C 2- Tenseur Q
  17. C
  18. C======================================================================
  19. C
  20. C CREATION : F.CORMERY
  21. C E.N.S.M.A - LMPM
  22. C JUILLET 1993
  23. C
  24. C======================================================================
  25. IMPLICIT INTEGER(I-N)
  26. IMPLICIT REAL*8 (A-H,O-Z)
  27. C**********************************************************************
  28. C DIMENSIONS ET DATA
  29. C**********************************************************************
  30. C
  31. CC DIMENSION TENS(6),VP(3)
  32. C N9 N15
  33. CC * ,P(3,3),S(6),
  34. C N24
  35. CC * VAL(3,3),TAMP(3,3)
  36. C
  37. DIMENSION TENS(*),VP(*)
  38. * ,P(3,*),S(*),
  39. * VAL(3,*),TAMP(3,*)
  40. DATA ZERO/0.D0/,UN/1.D0/,
  41. * PRECIS/1.D-08/,DPRECS/1.D-08/
  42. C----------------------------------------------------------------------
  43. AMAX1(X,Y,Z,U,V,W)=DMAX1(X,Y,Z,U,V,W)
  44. C------------------
  45. MT=10
  46. C**********************************************************************
  47. C INITIALISATION
  48. C**********************************************************************
  49. * DO 5 J=1,6
  50. * DO 6 K=1,6
  51. * P(J,K)=ZERO
  52. * Q(J,K)=ZERO
  53. * QPLUS(J,K)=ZERO
  54. *6 CONTINUE
  55. 5 CONTINUE
  56. DO 55 J=1,3
  57. DO 55 K=1,3
  58. 55 P(J,K)=ZERO
  59. C----------------------------------------------------------------------
  60. C**********************************************************************
  61. C NORMALISATION DU TENSEUR A
  62. C**********************************************************************
  63. C
  64. C----------------------------------------------------------------------
  65. C Trouver la valeur max de TENS(I)
  66. C----------------------------------------------------------------------
  67. DO 3 I=1,6
  68. S(I)=ABS(TENS(I))
  69. 3 CONTINUE
  70. C--------------
  71. TMAX=AMAX1(S(1),S(2),S(3),S(4),S(5),S(6))
  72. IF(TMAX.EQ.0.D0)TMAX=UN
  73. C----------------------------------------------------------------------
  74. C Normaliser a un la composante de TENS(I) la plus grande
  75. C----------------------------------------------------------------------
  76. DO 4 I=1,6
  77. TENS(I)=TENS(I)/TMAX
  78. IF(ABS(TENS(I)).LE.1E-15) TENS(I)=0.D0
  79. 4 CONTINUE
  80. C------------------------ cas axes principaux -------------------------
  81. NN=0
  82. DO 234 IV=4,6
  83. IF(ABS(TENS(IV)).LE.1E-15) NN=NN+1
  84. 234 CONTINUE
  85. IF(NN.EQ.3)THEN
  86. VP(1)=TENS(1)
  87. VP(2)=TENS(2)
  88. VP(3)=TENS(3)
  89. DO 235 I=1,3
  90. P(I,1)=0
  91. P(I,2)=0
  92. P(I,3)=0
  93. P(I,I)=1
  94. 235 CONTINUE
  95. goto 98
  96. ENDIF
  97. C***********************************************************************
  98. C CALCUL DES VALEURS PROPRES
  99. C***********************************************************************
  100. CALL VALPRP(TENS(1),TENS(2),TENS(3),TENS(6),TENS(4),TENS(5),
  101. * VP(1),VP(2),VP(3))
  102. C***********************************************************************
  103. C CALCUL DES VECTEURS PROPRES
  104. C***********************************************************************
  105. IM=2
  106. IF(ABS(VP(1)-VP(2)).LT.1E-8)THEN
  107. VP(1)=VP(3)
  108. VP(2)=VP(2)
  109. VP(3)=VP(2)
  110. IM=2
  111. ENDIF
  112. IF(ABS(VP(1)-VP(3)).LT.1E-8)THEN
  113. VP(1)=VP(2)
  114. VP(2)=VP(3)
  115. IM=2
  116. ENDIF
  117. IF(ABS(VP(2)-VP(3)).LT.1E-8)THEN
  118. IM=2
  119. ENDIF
  120. C----------------------------------------------------------------------
  121. IMM=0
  122. C----------------------------------------------------------------------
  123. DO 10 I=1,IM
  124. C-------------------
  125. SDET1=(TENS(2)-VP(I))*(TENS(3)-VP(I))-TENS(4)**2
  126. SDET2=(TENS(1)-VP(I))*(TENS(3)-VP(I))-TENS(5)**2
  127. SDET3=(TENS(1)-VP(I))*(TENS(2)-VP(I))-TENS(6)**2
  128. SDET4=TENS(6)*(TENS(3)-VP(I))-TENS(4)*TENS(5)
  129. SDET5=-TENS(5)*(TENS(2)-VP(I))+TENS(6)*TENS(4)
  130. SDET6=TENS(4)*(TENS(1)-VP(I))-TENS(6)*TENS(5)
  131. C----------------------------------------------------------------------
  132. C WRITE(10,*)'MINEURS :'
  133. C WRITE(10,*)SDET1,SDET2,SDET3,SDET4,SDET5,SDET6
  134. C----------------------------------------------------------------------
  135.  
  136. IF (ABS(SDET1).GT.DPRECS) THEN
  137. P(I,1)=UN
  138. P(I,2)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET1
  139. P(I,3)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET1
  140. GOTO 96
  141. C-------------------
  142. ENDIF
  143. IF (ABS(SDET2).GT.DPRECS) THEN
  144. P(I,2)=UN
  145. P(I,1)=((-TENS(6)*(TENS(3)-VP(I)))+TENS(4)*TENS(5))/SDET2
  146. P(I,3)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET2
  147. GOTO 96
  148. C-------------------
  149. ENDIF
  150. IF (ABS(SDET3).GT.DPRECS) THEN
  151. P(I,3)=UN
  152. P(I,1)=((-TENS(5)*(TENS(2)-VP(I)))+TENS(4)*TENS(6))/SDET3
  153. P(I,2)=((-TENS(4)*(TENS(1)-VP(I)))+TENS(5)*TENS(6))/SDET3
  154. GOTO 96
  155. C--------------------
  156. ENDIF
  157. IF (ABS(SDET4).GT.DPRECS) THEN
  158. P(I,1)=UN
  159. P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET4
  160. P(I,3)=((TENS(4)*(TENS(1)-VP(I)))-TENS(5)*TENS(6))/SDET4
  161. GOTO 96
  162. C--------------------
  163. ENDIF
  164. IF (ABS(SDET5).GT.DPRECS) THEN
  165. P(I,1)=(-(tens(4)**2)+(tens(2)-vp(i)))/sdet5
  166. P(I,2)=(-tens(6)*(tens(3)-vp(i))+tens(5)*tens(4))/sdet5
  167. P(I,3)=1
  168. GOTO 96
  169. C--------------------
  170. ENDIF
  171. IF (ABS(SDET6).GT.DPRECS) THEN
  172. P(I,3)=UN
  173. P(I,1)=((-TENS(5)*TENS(4))+(TENS(3)-vp(i))*TENS(6))/SDET6
  174. P(I,2)=((-(TENS(3)-vp(i))*(TENS(1)-VP(I)))+TENS(5)**2)/SDET6
  175. ENDIF
  176. C-----------------------------------------------------------------------
  177. SSDET1=TENS(1)-VP(I)
  178. SSDET2=TENS(2)-VP(I)
  179. SSDET3=TENS(3)-VP(I)
  180. C-------------------CAS PARTICULIERS------------------------------------
  181. IF (ABS(SSDET1).LE.PRECIS) THEN
  182. P(I,1)=1
  183. P(I,2)=0
  184. P(I,3)=0
  185. GOTO 96
  186. ENDIF
  187. IF (ABS(SSDET2).LE.PRECIS) THEN
  188. P(I,1)=0
  189. P(I,2)=1
  190. P(I,3)=0
  191. GOTO 96
  192. ENDIF
  193. IF (ABS(SSDET3).LE.PRECIS) THEN
  194. P(I,1)=0
  195. P(I,2)=0
  196. P(I,3)=1
  197. GOTO 96
  198. ENDIF
  199. IF (ABS(SSDET1).GT.PRECIS) THEN
  200. P(I,1)=-(TENS(6)+TENS(5))/SSDET1
  201. P(I,2)=1
  202. P(I,3)=1
  203. C-------------------
  204. GOTO 96
  205. ENDIF
  206. IF (ABS(SSDET2).GT.PRECIS) THEN
  207. P(I,1)=1
  208. P(I,2)=-(TENS(6)+TENS(4))/SSDET2
  209. P(I,3)=1
  210. GOTO 96
  211. C-------------------
  212. ENDIF
  213. IF (ABS(SSDET3).GT.PRECIS) THEN
  214. P(I,1)=1
  215. P(I,2)=1
  216. P(I,3)=-(TENS(5)+TENS(4))/SSDET3
  217. GOTO 96
  218. ENDIF
  219. C-----------------------------------------------------------------------
  220. WRITE(MT,*)'ERREUR DANS VPROP.FOR'
  221. WRITE(MT,1010)
  222. 1010 FORMAT(1X,'TENSEUR A SYM. D ORDRE 2 :',
  223. * /1X,'--------------------------'/)
  224. WRITE(MT,1001)TENS(1),TENS(6),TENS(5)
  225. 1001 FORMAT(15X,'* ',3e20.7,' *')
  226. WRITE(MT,1002)TENS(6),TENS(2),TENS(4)
  227. 1002 FORMAT(15X,'* ',3e20.7,' *')
  228. WRITE(MT,1003)TENS(5),TENS(4),TENS(3)
  229. 1003 FORMAT(15X,'* ',3e20.7,' *'/)
  230. STOP
  231. C-----------------------------------------------------------------------
  232. 96 CONTINUE
  233. C-----------------------------------------------------------------------
  234. 10 CONTINUE
  235. 98 CONTINUE
  236. IF(IM.EQ.2)THEN
  237. P(3,1)=P(1,2)*P(2,3)-P(1,3)*P(2,2)
  238. P(3,2)=P(1,3)*P(2,1)-P(1,1)*P(2,3)
  239. P(3,3)=P(1,1)*P(2,2)-P(1,2)*P(2,1)
  240. ENDIF
  241. C------------------------------------------------------------------------
  242. C On verifie que la base formee est bien directe
  243. C------------------------------------------------------------------------
  244. DIR=P(3,1)*(P(1,2)*P(2,3)-P(1,3)*P(2,2))+
  245. * P(3,2)*(P(1,3)*P(2,1)-P(1,1)*P(2,3))+
  246. * P(3,3)*(P(1,1)*P(2,2)-P(1,2)*P(2,1))
  247. IF(DIR.LT.ZERO)THEN
  248. DO 12 J=1,3
  249. TAMP1=VP(2)
  250. VP(2)=VP(1)
  251. VP(1)=TAMP1
  252. TAMP(1,J)=P(1,J)
  253. P(1,J)=P(2,J)
  254. P(2,J)=TAMP(1,J)
  255. 12 CONTINUE
  256. ENDIF
  257. C------------------------------------------------------------------------
  258. C Normalisation des vecteurs propres
  259. C------------------------------------------------------------------------
  260. DO 11 J=1,5
  261. DO 11 I=1,3
  262. IF(ABS(P(I,1)).LT.1.E-15)P(I,1)=0.D0
  263. IF(ABS(P(I,2)).LT.1.E-15)P(I,2)=0.D0
  264. IF(ABS(P(I,3)).LT.1.E-15)P(I,3)=0.D0
  265. RAC=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  266. P(I,1)=P(I,1)/RAC
  267. P(I,2)=P(I,2)/RAC
  268. P(I,3)=P(I,3)/RAC
  269. 11 CONTINUE
  270. C-----------------------------------------------------------------------
  271. C Multiplication par le facteur de normalisation
  272. C-----------------------------------------------------------------------
  273. DO 110 I=1,6
  274. TENS(I)=TENS(I)*TMAX
  275. 110 CONTINUE
  276. C-----------------------------------------------------------------------
  277. DO 341 I=1,3
  278. DO 341 J=1,3
  279. VAL(I,J)=P(I,J)
  280. 341 CONTINUE
  281. DO 342 I=1,3
  282. VP(I)=VP(I)*tmax
  283. 342 CONTINUE
  284. C-------------------------------------------------------------------------
  285. RETURN
  286. END
  287. C-----------------------------------------------------------------------
  288.  
  289.  
  290.  

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