Télécharger dcov.eso

Retour à la liste

Numérotation des lignes :

  1. C DCOV SOURCE PV 13/04/12 21:15:27 7756
  2. SUBROUTINE DCOV
  3. C
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8(A-H,O-Z)
  6. C
  7. -INC CCOPTIO
  8. -INC SMELEME
  9. -INC SMCOORD
  10. -INC SMRIGID
  11. C
  12. DIMENSION ZLAM(3),ZINVL(3),XCOV(3,3)
  13. CHARACTER*8 MOTS(1),MOTLA(4)
  14. CHARACTER*4 MOTLO(2)
  15. CHARACTER*9 MOTD(1)
  16. C
  17. DATA MOTS/'SIGMA '/
  18. DATA MOTLA/'LAMBDA ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/
  19. DATA MOTLO/'EXPO','GAUS'/
  20. DATA MOTD/'DIRECTION'/
  21. C
  22. C epsilon servant à différents tests
  23. C
  24. EPS = 1.D-12
  25. C
  26. IDIR=0
  27. C
  28. C Lecture du maillage
  29. C
  30. CALL LIROBJ('MAILLAGE',MELEME,1,IRET)
  31. IF (IERR.NE.0) RETURN
  32. C
  33. C Lecture de la loi utilisée
  34. C
  35. CALL LIRMOT(MOTLO,2,ILOI,1)
  36. IF (ILOI.EQ.0) RETURN
  37. C
  38. C Lecture du mot-clé 'SIGMA'
  39. C
  40. CALL LIRMOT(MOTS,1,IRET,1)
  41. IF (IRET.EQ.0) RETURN
  42. C
  43. C Lecture de la valeur de sigma (strictement supérieure à 0.D0)
  44. C
  45. CALL LIRREE(ZSIG,1,IRET)
  46. IF (IERR.NE.0) RETURN
  47. IF (ZSIG.LE.0.D0) THEN
  48. REAERR(1) = REAL(ZSIG)
  49. REAERR(2) = REAL(0.D0)
  50. MOTERR(1:8) = ' SIGMA '
  51. CALL ERREUR(41)
  52. RETURN
  53. ENDIF
  54. C
  55. C Lecture du mot-clé 'LAMBDA' ou 'LAMBDA1'
  56. C
  57. IF (IDIM.EQ.1) THEN
  58. CALL LIRMOT(MOTLA,1,IANI,1)
  59. ELSE
  60. CALL LIRMOT(MOTLA,2,IANI,1)
  61. ENDIF
  62. IF (IANI.EQ.0) RETURN
  63. C
  64. C Lecture de la valeur de lambda (ou lambda1)
  65. C (strictement supérieure à 0.D0)
  66. C
  67. CALL LIRREE(ZLAM(1),1,IRET)
  68. IF (IERR.NE.0) RETURN
  69. IF (ZLAM(1).LE.0.D0) THEN
  70. REAERR(1) = REAL(ZLAM(1))
  71. REAERR(2) = REAL(0.D0)
  72. IF (IANI.EQ.1) MOTERR(1:8) = ' LAMBDA '
  73. IF (IANI.EQ.2) MOTERR(1:8) = ' LAMBDA1'
  74. CALL ERREUR(41)
  75. RETURN
  76. ENDIF
  77. C
  78. IF (IANI.EQ.2) THEN
  79. C
  80. C Lecture éventuelle du mot-clé 'LAMBDA2' (cas anisotrope)
  81. C
  82. CALL LIRMOT(MOTLA(3),1,IRET,1)
  83. IF (IRET.EQ.0) RETURN
  84. C
  85. C Lecture éventuelle de la valeur de lambda2 (cas anisotrope)
  86. C (strictement supérieure à 0.D0)
  87. C
  88. CALL LIRREE(ZLAM(2),1,IRET)
  89. IF (IERR.NE.0) RETURN
  90. IF (ZLAM(2).LE.0.D0) THEN
  91. REAERR(1) = REAL(ZLAM(2))
  92. REAERR(2) = REAL(0.D0)
  93. MOTERR(1:8) = ' LAMBDA2'
  94. CALL ERREUR(41)
  95. RETURN
  96. ENDIF
  97. C
  98. IF (IDIM.EQ.3) THEN
  99. C
  100. C Lecture éventuelle du mot-clé 'LAMBDA3' (cas anisotrope en 3d)
  101. C
  102. CALL LIRMOT(MOTLA(4),1,IRET,1)
  103. IF (IRET.EQ.0) RETURN
  104. C
  105. C Lecture éventuelle de la valeur de lambda3 (cas anisotrope en 3d)
  106. C (strictement supérieure à 0.D0)
  107. C
  108. CALL LIRREE(ZLAM(3),1,IRET)
  109. IF (IERR.NE.0) RETURN
  110. IF (ZLAM(3).LE.0.D0) THEN
  111. REAERR(1) = REAL(ZLAM(3))
  112. REAERR(2) = REAL(0.D0)
  113. MOTERR(1:8) = ' LAMBDA3'
  114. CALL ERREUR(41)
  115. RETURN
  116. ENDIF
  117. C
  118. ENDIF
  119. C
  120. C Lecture éventuelle et optionnelle du mot-clé 'DIRECTION'
  121. C (cas anisotrope, directions d'anisotropie non parallèles aux axes)
  122. C
  123. CALL LIRMOT(MOTD,1,IDIR,0)
  124. IF (IERR.NE.0) RETURN
  125. C
  126. IF (IDIR.NE.0) THEN
  127. C
  128. C Lecture éventuelle de la valeur de VEC1
  129. C (cas anisotrope, directions d'anisotropie non parallèles aux axes)
  130. C
  131. CALL LIROBJ('POINT',IPTV1,1,IRET)
  132. IF (IERR.NE.0) RETURN
  133. C
  134. IF (IDIM.EQ.3) THEN
  135. C
  136. C Lecture éventuelle de la valeur de VEC2
  137. C (cas anisotrope en 3d,
  138. C directions d'anisotropie non parallèles aux axes)
  139. C
  140. CALL LIROBJ('POINT',IPTV2,1,IRET)
  141. IF (IERR.NE.0) RETURN
  142. C
  143. ENDIF
  144. C
  145. ENDIF
  146. C
  147. ENDIF
  148. C
  149. C Dans le cas anisotrope,
  150. C directions d'anisotropie non parallèles aux axes,
  151. C on vérifie que le(s) vecteur(s) introduit(s) sous forme de point ne
  152. C sont pas de longueur nulle
  153. C si tel est le cas, on le(s) norme
  154. C
  155. IF (IDIR.NE.0) THEN
  156. C
  157. IREF1 = (IPTV1-1)*(IDIM+1)
  158. SCOV1 = 0.D0
  159. DO 1 I=1,IDIM
  160. XCOV(1,I) = XCOOR(IREF1+I)
  161. SCOV1 = SCOV1+(XCOV(1,I)*XCOV(1,I))
  162. 1 CONTINUE
  163. C
  164. IF (SCOV1.LT.EPS) THEN
  165. CALL ERREUR(239)
  166. RETURN
  167. ELSE
  168. DO 2 I=1,IDIM
  169. XCOV(1,I) = XCOV(1,I)/SCOV1
  170. 2 CONTINUE
  171. ENDIF
  172. C
  173. IF (IDIM.EQ.3) THEN
  174. C
  175. IREF2 = (IPTV2-1)*(IDIM+1)
  176. SCOV2 = 0.D0
  177. DO 3 I=1,IDIM
  178. XCOV(2,I) = XCOOR(IREF2+I)
  179. SCOV2 = SCOV2+(XCOV(2,I)*XCOV(2,I))
  180. 3 CONTINUE
  181. C
  182. IF (SCOV2.LT.EPS) THEN
  183. CALL ERREUR(239)
  184. RETURN
  185. ELSE
  186. DO 4 I=1,IDIM
  187. XCOV(2,I) = XCOV(2,I)/SCOV2
  188. 4 CONTINUE
  189. ENDIF
  190. C
  191. C Dans le cas anisotrope 3d,
  192. C directions d'anisotropie non parallèles aux axes,
  193. C on vérifie que les vecteurs introduits sous forme de point
  194. C sont bien orthogonaux
  195. C
  196. SCOV12 = 0.D0
  197. DO 5 I=1,IDIM
  198. SCOV12 = SCOV12+(XCOV(1,I)*XCOV(2,I))
  199. 5 CONTINUE
  200. IF (SCOV12.GT.EPS) THEN
  201. CALL ERREUR(161)
  202. RETURN
  203. ENDIF
  204. C
  205. ENDIF
  206. C
  207. ENDIF
  208. C
  209. IF (IDIR.NE.0) THEN
  210. C
  211. C Dans le cas anisotrope 2d,
  212. C directions d'anisotropie non parallèles aux axes,
  213. C on complète le repère par rotation de +90 degrés
  214. C
  215. IF (IDIM.EQ.2) THEN
  216. XCOV(2,1) = -1.D0*XCOV(1,2)
  217. XCOV(2,2) = XCOV(1,1)
  218. ENDIF
  219. C
  220. C Dans le cas anisotrope 3d,
  221. C directions d'anisotropie non parallèles aux axes,
  222. C on complète le repère par produit vectoriel des deux premiers
  223. C vecteurs
  224. C
  225. IF (IDIM.EQ.3) THEN
  226. XCOV(3,1) = (XCOV(1,2)*XCOV(2,3))-(XCOV(1,3)*XCOV(2,2))
  227. XCOV(3,2) = (XCOV(1,3)*XCOV(2,1))-(XCOV(1,1)*XCOV(2,3))
  228. XCOV(3,3) = (XCOV(1,1)*XCOV(2,2))-(XCOV(1,2)*XCOV(2,1))
  229. ENDIF
  230. C
  231. ENDIF
  232. C
  233. C Il vaut mieux faire * (1./lambda) que / lambda
  234. C
  235. ZINVL(1) = 1.D0/ZLAM(1)
  236. IF (IANI.EQ.2) THEN
  237. ZINVL(2) = 1.D0/ZLAM(2)
  238. IF (IDIM.EQ.3) ZINVL(3) = 1.D0/ZLAM(3)
  239. ENDIF
  240. C
  241. C Les éléments du MELEME sont transformés en POI1 si nécessaire
  242. C
  243. SEGACT MELEME
  244. IPT1=MELEME
  245. IPT2=MELEME
  246. IF(ITYPEL.NE.1.OR.(LISOUS(/1)).NE.0) THEN
  247. CALL CHANGE(IPT1,1)
  248. IF (IERR.NE.0) RETURN
  249. ENDIF
  250. C
  251. NBLI=IPT1.NUM(/2)
  252. C
  253. C On définit un super-élément correspondant aux noeuds
  254. C
  255. NBSOUS = 0
  256. NBREF = 0
  257. NBNN = NBLI
  258. NBELEM = 1
  259. SEGINI MELEME
  260. ITYPEL = 28
  261. DO 7 I=1,NBLI
  262. NUM(I,1) = IPT1.NUM(1,I)
  263. 7 CONTINUE
  264. SEGDES IPT2
  265. C
  266. C On met la matrice dans un objet rigidité
  267. C
  268. NRIGE = 8
  269. NRIGEL = 1
  270. SEGINI MRIGID
  271. ICHOLE = 0
  272. IMGEO1 = 0
  273. IMGEO2 = 0
  274. IFORIG = IFOUR
  275. ISUPEQ = 0
  276. NELRIG = 1
  277. COERIG(1) = 1.D0
  278. * SEGINI IMATRI
  279. NLIGRD = NBLI
  280. NLIGRP = NBLI
  281. segini xmatri
  282. SEGINI DESCR
  283. IRIGEL(1,1) = MELEME
  284. IRIGEL(2,1) = 0
  285. IRIGEL(3,1) = DESCR
  286. IRIGEL(4,1) = xMATRI
  287. IRIGEL(5,1) = 0
  288. IRIGEL(6,1) = 0
  289. DO 8 I=1,NBLI
  290. NOELEP(I) = I
  291. NOELED(I) = I
  292. LISINC(I) = 'SCAL'
  293. LISDUA(I) = 'SCAL'
  294. 8 CONTINUE
  295. SEGDES DESCR,MRIGID
  296. * SEGINI XMATRI
  297. * IMATTT(1) = XMATRI
  298. * DO 9 I=1,NBLI
  299. * DO 10 J=1,NBLI
  300. * RE(I,J,1) = 0.D0
  301. * 10 CONTINUE
  302. * 9 CONTINUE
  303. C
  304. IF (ILOI.EQ.1) THEN
  305. C
  306. C Loi exponentielle
  307. C
  308. IF (IANI.EQ.1) THEN
  309. C
  310. C Cas isotrope
  311. C
  312. DO 100 I=1,NBLI
  313. IREF=(NUM(I,1)-1)*(IDIM+1)
  314. DO 110 J=1,I
  315. JREF=(NUM(J,1)-1)*(IDIM+1)
  316. VM=0.D0
  317. DO 120 I1=1,IDIM
  318. VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  319. 120 CONTINUE
  320. VM=SQRT(VM)
  321. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  322. 110 CONTINUE
  323. 100 CONTINUE
  324. C
  325. ELSEIF (IANI.EQ.2) THEN
  326. C
  327. C Cas anisotrope :
  328. C directions d'anisotropie parallèles aux axes (IDIR nul)
  329. C
  330. IF (IDIR.EQ.0) THEN
  331. C
  332. DO 200 I=1,NBLI
  333. IREF=(NUM(I,1)-1)*(IDIM+1)
  334. DO 210 J=1,I
  335. JREF=(NUM(J,1)-1)*(IDIM+1)
  336. VM=0.D0
  337. DO 220 I1=1,IDIM
  338. VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  339. 220 CONTINUE
  340. VM=SQRT(VM)
  341. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  342. 210 CONTINUE
  343. 200 CONTINUE
  344. C
  345. C Cas anisotrope :
  346. C directions d'anisotropie non parallèles aux axes (IDIR non nul)
  347. C
  348. ELSEIF (IDIR.NE.0) THEN
  349. C
  350. DO 300 I=1,NBLI
  351. IREF=(NUM(I,1)-1)*(IDIM+1)
  352. DO 310 J=1,I
  353. JREF=(NUM(J,1)-1)*(IDIM+1)
  354. VM=0.D0
  355. DO 320 I1=1,IDIM
  356. PSCAL=0.D0
  357. DO 330 J1=1,IDIM
  358. PSCAL = PSCAL +
  359. * (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  360. 330 CONTINUE
  361. VM=VM+(PSCAL*ZINVL(I1))**2
  362. 320 CONTINUE
  363. VM=SQRT(VM)
  364. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  365. 310 CONTINUE
  366. 300 CONTINUE
  367. C
  368. ENDIF
  369. C
  370. ENDIF
  371. C
  372. ELSEIF (ILOI.EQ.2) THEN
  373. C
  374. C Loi gaussienne
  375. C
  376. IF (IANI.EQ.1) THEN
  377. C
  378. C Cas isotrope
  379. C
  380. DO 400 I=1,NBLI
  381. IREF=(NUM(I,1)-1)*(IDIM+1)
  382. DO 410 J=1,I
  383. JREF=(NUM(J,1)-1)*(IDIM+1)
  384. VM=0.D0
  385. DO 420 I1=1,IDIM
  386. VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  387. 420 CONTINUE
  388. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  389. 410 CONTINUE
  390. 400 CONTINUE
  391. C
  392. ELSEIF (IANI.EQ.2) THEN
  393. C
  394. C Cas anisotrope :
  395. C directions d'anisotropie parallèles aux axes (IDIR nul)
  396. C
  397. IF (IDIR.EQ.0) THEN
  398. C
  399. DO 500 I=1,NBLI
  400. IREF=(NUM(I,1)-1)*(IDIM+1)
  401. DO 510 J=1,I
  402. JREF=(NUM(J,1)-1)*(IDIM+1)
  403. VM=0.D0
  404. DO 520 I1=1,IDIM
  405. VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  406. 520 CONTINUE
  407. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  408. 510 CONTINUE
  409. 500 CONTINUE
  410. C
  411. C Cas anisotrope :
  412. C directions d'anisotropie non parallèles aux axes (IDIR non nul)
  413. C
  414. ELSEIF (IDIR.NE.0) THEN
  415. C
  416. DO 600 I=1,NBLI
  417. IREF=(NUM(I,1)-1)*(IDIM+1)
  418. DO 610 J=1,I
  419. JREF=(NUM(J,1)-1)*(IDIM+1)
  420. VM=0.D0
  421. DO 620 I1=1,IDIM
  422. PSCAL=0.D0
  423. DO 630 J1=1,IDIM
  424. PSCAL = PSCAL +
  425. * (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  426. 630 CONTINUE
  427. VM=VM+(PSCAL*ZINVL(I1))**2
  428. 620 CONTINUE
  429. RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  430. 610 CONTINUE
  431. 600 CONTINUE
  432. C
  433. ENDIF
  434. C
  435. ENDIF
  436. C
  437. ENDIF
  438. C
  439. C Calcul de M telle que (RE =) COV = M Mt
  440. C M est stockée dans la partie inférieure de RE
  441. C
  442. CALL SYCHOM(RE,NBLI)
  443. C
  444. C Création de l'objet rigidité
  445. C
  446. CALL ECROBJ('RIGIDITE',MRIGID)
  447. C
  448. SEGDES XMATRI,MELEME
  449. C
  450. END
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  

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