Télécharger epplum.eso

Retour à la liste

Numérotation des lignes :

  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