Télécharger epplum.eso

Retour à la liste

Numérotation des lignes :

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

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