Télécharger rleca1.eso

Retour à la liste

Numérotation des lignes :

  1. C RLECA1 SOURCE CHAT 05/01/13 03:01:02 5004
  2. SUBROUTINE RLECA1(MELFL,MELFP,MLEFSC,MACOE1)
  3. C************************************************************************
  4. C
  5. C PROJET : CASTEM 2000
  6. C
  7. C NOM : RLECA1
  8. C
  9. C DESCRIPTION : Appelle par GRADIA
  10. C
  11. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  12. C
  13. C AUTEUR : A. BECCANTINI
  14. C
  15. C************************************************************************
  16. C
  17. C Inputs:
  18. C
  19. C MELFL : MELEME FACEL
  20. C
  21. C MELFP : MELEME FACEP
  22. C
  23. C Outputs:
  24. C
  25. C MLEFSC : list of FACES points and their neighbors (CENTRE and SOMMET
  26. C points)
  27. C
  28. C MACOE1 : coeff. for linear exact reconstruction of MLEFSC
  29. C
  30. C
  31. C**** Variables de COOPTIO
  32. C
  33. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  34. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  35. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  36. C & ,IECHO, IIMPI, IOSPI
  37. C & ,IDIM
  38. CC & ,MCOORD
  39. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  40. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  41. C & ,NORINC,NORVAL,NORIND,NORVAD
  42. C & ,NUCROU, IPSAUV, IFICLE, IPREFI
  43. C
  44. IMPLICIT INTEGER(I-N)
  45. -INC CCOPTIO
  46. -INC SMCOORD
  47. -INC SMELEME
  48. POINTEUR MELFL.MELEME,MELFP.MELEME,MELFP1.MELEME
  49. INTEGER JG
  50. -INC SMLENTI
  51. POINTEUR MLEFP.MLENTI
  52. -INC SMLREEL
  53. POINTEUR MLRSIG.MLREEL, MLRTRA.MLREEL
  54. & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
  55. & ,MLRIN3.MLREEL
  56. C
  57. INTEGER NBL, NBTPOI
  58. SEGMENT MLELEM
  59. INTEGER INDEX(NBL+1)
  60. INTEGER LESPOI(NBTPOI)
  61. ENDSEGMENT
  62. POINTEUR MLEFSC.MLELEM
  63. C
  64. INTEGER N1,N2
  65. SEGMENT MATRIX
  66. REAL*8 MAT(N1,N2)
  67. ENDSEGMENT
  68. POINTEUR MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
  69. & ,MACOE1.MATRIX
  70. LOGICAL LOGSVD
  71. C
  72. INTEGER NBSOUS, NVMAX, NELT, NNOE, ISOUS, NVMIN, NVOI, IELEM
  73. & , IPOS, IELEM1,NGF,NGF1, NGN, IPCOOR, I1, I2, I3, IPVOI
  74. & , IERSVD, NGCD, NGCG, IVOI,IERR0,I4
  75.  
  76. REAL*8 ERRTOL,SMAX,SMIN, XF, YF, ZF
  77. PARAMETER (ERRTOL=1.0D0-15)
  78. C
  79. NBL = 0
  80. NBTPOI = 0
  81. NVMAX = 0
  82. SEGACT MELFP
  83. NBSOUS=MELFP.LISOUS(/1)
  84. C NBSOUS=0 fais un peux chier!
  85. JG=MAX(NBSOUS,1)
  86. SEGINI MLEFP
  87. IF(NBSOUS .EQ. 0)THEN
  88. MLEFP.LECT(1)=MELFP
  89. ELSE
  90. DO ISOUS=1,NBSOUS,1
  91. MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS)
  92. ENDDO
  93. ENDIF
  94. SEGDES MELFP
  95. NBSOUS=JG
  96. C
  97. DO ISOUS=1,NBSOUS,1
  98. MELFP1=MLEFP.LECT(ISOUS)
  99. SEGACT MELFP1
  100. NNOE=MELFP1.NUM(/1)-1
  101. NELT=MELFP1.NUM(/2)
  102. NBTPOI=NBTPOI+((1+2+NNOE)*NELT)
  103. C
  104. C 1 = centre de face
  105. C 2 = centre d'elts "gauche" et "droite" (si differents)
  106. C NNOE = sommet voisins
  107. C
  108. NBL=NBL+NELT
  109. NVMAX=MAX(NVMAX,(2+NNOE))
  110. ENDDO
  111. C
  112. C N.B.: NBTPOI surestimé
  113. C
  114. C**** MLEFSC
  115. C Structure Face-(Sommets-Centres)
  116. C
  117. SEGINI MLEFSC
  118. C
  119. C**** La matrice de coeff.
  120. C
  121. C MACOE1(*,I) = coefficients concernants MLEFSC.LESPOI(I)
  122. C
  123. NVMIN=IDIM+1
  124. C
  125. C**** NVMIN = nombre minimum de voisins au dessus duquel MATA est
  126. C singulier
  127. C
  128. N1=NVMIN
  129. N2=NBTPOI
  130. SEGINI MACOE1
  131. C
  132. C**** MATA = matrice à inverser (NVMAX,IDIM+1)
  133. C MATU = matice des vectors propres singuliers droite de MATA
  134. C (NVMAX,IDIM+1)
  135. C MATV = matrice des vectors propres singuliers gauche de MATA
  136. C (IDIM+1,IDIM+1)
  137. C Mais en INVSVD a dimension NVMAX,IDIM+1
  138. C MLRSIG =liste des valeurs singuliers de MATA
  139. C
  140. C N.B. MATA = MATU MATSIG t(MATV)
  141. C Si MATA non singulier
  142. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  143. C
  144. N1=NVMAX
  145. N2=NVMIN
  146. SEGINI MATA
  147. SEGINI MATU
  148. SEGINI MATV
  149. JG=NVMIN
  150. SEGINI MLRSIG
  151. SEGINI MLRTRA
  152. C MLRTRA segment de travail de la subroutine invsvd
  153. C
  154. C MTAA : [t(MATA).MATA]
  155. C MINTAA : [inve(t(MATA) MATA)]
  156. C MLRIN1,2,3 : "temporary vectors"
  157. C
  158. N1=NVMIN
  159. N2=NVMIN
  160. SEGINI MTAA
  161. SEGINI MINTAA
  162. JG=NVMIN
  163. SEGINI MLRIN1
  164. SEGINI MLRIN2
  165. SEGINI MLRIN3
  166. C
  167. SEGACT MELFL
  168. MLEFSC.INDEX(1)=1
  169. IELEM=0
  170. DO ISOUS=1,NBSOUS,1
  171. MELFP1=MLEFP.LECT(ISOUS)
  172. NNOE=MELFP1.NUM(/1)-1
  173. NELT=MELFP1.NUM(/2)
  174. DO IELEM1=1,NELT,1
  175. IELEM=IELEM+1
  176. IPOS=MLEFSC.INDEX(IELEM)
  177. NGF=MELFP1.NUM(NNOE+1,IELEM1)
  178. NGF1=MELFL.NUM(2,IELEM)
  179. IF(NGF .NE. NGF1)THEN
  180. WRITE(IOIMP,*) 'subroutine rleca1.eso'
  181. CALL ERREUR(5)
  182. GOTO 9999
  183. ENDIF
  184. IPCOOR=(IDIM+1)*(NGF-1)+1
  185. XF=MCOORD.XCOOR(IPCOOR)
  186. YF=MCOORD.XCOOR(IPCOOR+1)
  187. IF(IDIM .EQ. 3) ZF=MCOORD.XCOOR(IPCOOR+2)
  188. MLEFSC.LESPOI(IPOS)=NGF
  189. DO I1 = 1, NVMIN, 1
  190. MACOE1.MAT(I1,IPOS)=0.0D0
  191. ENDDO
  192. C
  193. C********** Le centre "gauche"
  194. C
  195. NGCG=MELFL.NUM(1,IELEM)
  196. NGN=NGCG
  197. C NVOI=nombre de voisins
  198. NVOI=1
  199. MLEFSC.LESPOI(IPOS+1)=NGN
  200. MATA.MAT(1,1)=1.0D0
  201. IPCOOR=(IDIM+1)*(NGN-1)+1
  202. MATA.MAT(1,2)=MCOORD.XCOOR(IPCOOR)-XF
  203. MATA.MAT(1,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  204. IF(IDIM.EQ.3) MATA.MAT(1,4)=MCOORD.XCOOR(IPCOOR+2)-ZF
  205. C
  206. C********** Le centre "droite"
  207. C
  208. NGCD=MELFL.NUM(3,IELEM)
  209. IF(NGCD .NE. NGCG)THEN
  210. NVOI=NVOI+1
  211. NGN=NGCD
  212. MLEFSC.LESPOI(IPOS+2)=NGN
  213. MATA.MAT(2,1)=1.0D0
  214. IPCOOR=(IDIM+1)*(NGN-1)+1
  215. MATA.MAT(2,2)=MCOORD.XCOOR(IPCOOR)-XF
  216. MATA.MAT(2,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  217. IF(IDIM.EQ.3) MATA.MAT(2,4)=MCOORD.XCOOR(IPCOOR+2)-ZF
  218. ENDIF
  219. C
  220. C********** Les sommets
  221. C
  222. DO I1 = 1, NNOE, 1
  223. NGN=MELFP1.NUM(I1,IELEM1)
  224. MLEFSC.LESPOI(IPOS+NVOI+I1)=NGN
  225. MATA.MAT(NVOI+I1,1)=1.0D0
  226. IPCOOR=(IDIM+1)*(NGN-1)+1
  227. MATA.MAT(NVOI+I1,2)=MCOORD.XCOOR(IPCOOR)-XF
  228. MATA.MAT(NVOI+I1,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  229. IF(IDIM.EQ.3) MATA.MAT(NVOI+I1,4)=MCOORD.XCOOR(IPCOOR+2)
  230. & -ZF
  231. ENDDO
  232. NVOI=NVOI+NNOE
  233. MLEFSC.INDEX(IELEM+1)=IPOS+NVOI+1
  234. C
  235. LOGSVD=.TRUE.
  236. CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  237. & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
  238. & MLRTRA.PROG)
  239. IF(IERSVD.NE.0)THEN
  240. LOGSVD=.FALSE.
  241. ELSE
  242. SMAX=0.0D0
  243. DO I2=1,NVMIN,1
  244. SMAX=MAX(SMAX,MLRSIG.PROG(I2))
  245. ENDDO
  246. SMIN=SMAX
  247. DO I2=1,NVMIN,1
  248. SMIN=MIN(SMIN,MLRSIG.PROG(I2))
  249. ENDDO
  250. ENDIF
  251. IF((SMIN/SMAX) .LT. ERRTOL)THEN
  252. LOGSVD=.FALSE.
  253. ENDIF
  254. C
  255. C********** We have two different cases:
  256. C LOGSVD = .TRUE. -> we have the coeff. we are looking for
  257. C LOGSVD = .FALSE. -> we have not
  258. C
  259. C LOGSVD=.FALSE.
  260. IF(LOGSVD)THEN
  261. DO I4=1,NVMIN,1
  262. DO I3=1,NVOI,1
  263. IPVOI=IPOS+I3
  264. DO I2=1,NVMIN,1
  265. MACOE1.MAT(I2,IPVOI)=MACOE1.MAT(I2,IPVOI)+
  266. & (MATV.MAT(I2,I4)*MATU.MAT(I3,I4)/
  267. & MLRSIG.PROG(I4))
  268. ENDDO
  269. ENDDO
  270. ENDDO
  271. ELSE
  272. WRITE (IOIMP,*) 'rleca1.eso'
  273. C 22 0
  274. C Opération malvenue. Résultat douteux
  275. CALL ERREUR(22)
  276. C
  277. C************* INVSVD does not worked
  278. C For each neighbor k we have to compute the solution
  279. C of
  280. C
  281. C t(MATA) MATA x = t(MATA) * b
  282. C
  283. C where b= \sum_l e_l \delta(k,l) = e_k
  284. C
  285. C To do that, we compute
  286. C
  287. C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  288. C
  289. C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  290. C [t(MATA) MATA] X_0]
  291. C
  292. C
  293. C************* We compute [t(MATA) MATA]
  294. C We store it in the upper triangle of MTAA(NVMIN,NVMIN)
  295. C
  296. DO I1=1,NVMIN,1
  297. DO I2=I1,NVMIN,1
  298. MTAA.MAT(I1,I2)=0.0D0
  299. ENDDO
  300. ENDDO
  301. C
  302. DO I1=1,NVMIN,1
  303. DO I2=I1,NVMIN,1
  304. DO I3=1,NVOI,1
  305. MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
  306. & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
  307. ENDDO
  308. ENDDO
  309. ENDDO
  310. C
  311. C************* We compute [inve(t(MATA) MATA)]
  312. C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
  313. C
  314. CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
  315. & MLRIN2.PROG,IERR0)
  316. IF(IERR0.NE.0)THEN
  317. WRITE(IOIMP,*) 'subroutine rleca1.eso.'
  318. C 26 2
  319. C Tache impossible. Probablement données erronées
  320. CALL ERREUR(26)
  321. GOTO 9999
  322. ENDIF
  323. C
  324. C************* We complete MTAA and MINTAA
  325. C
  326. DO I1=1,NVMIN,1
  327. DO I2=I1+1,NVMIN,1
  328. MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
  329. MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
  330. ENDDO
  331. ENDDO
  332. C
  333. DO IVOI=1,NVOI,1
  334. C
  335. C**************** We compute [t(MATA) . b] and we store it in MLRIN1.PROG
  336. C
  337. DO I1=1,NVMIN,1
  338. MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)
  339. MLRIN2.PROG(I1)=0.0D0
  340. MLRIN3.PROG(I1)=0.0D0
  341. ENDDO
  342. C
  343. C**************** X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  344. C X_0(i1) into MLRIN2.PROG(I1)
  345. C
  346. DO I2=1,NVMIN,1
  347. DO I1=1,NVMIN,1
  348. MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
  349. & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
  350. ENDDO
  351. ENDDO
  352. C
  353. C**************** X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  354. C [t(MATA) MATA] X_0]
  355. C
  356. C [t(MATA) MATA] X_0 into MLRIN3.PROG
  357. C
  358. DO I2=1,NVMIN,1
  359. DO I1=1,NVMIN,1
  360. MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
  361. & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
  362. ENDDO
  363. ENDDO
  364. C
  365. C**************** Now we have
  366. C [t(MATA) . b] in MLRIN1.PROG
  367. C X_0(i1) in MLRIN2.PROG(I1)
  368. C [t(MATA) MATA] X_0 in MLRIN3.PROG
  369. C
  370. C X_1(i1) in MLRCOE.MAT(i1,IVOI)
  371. C
  372. IPVOI=IPOS+IVOI
  373. DO I1=1,NVMIN,1
  374. DO I2=1,NVMIN,1
  375. MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+
  376. & (MINTAA.MAT(I1,I2)*
  377. & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
  378. ENDDO
  379. MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+
  380. & MLRIN2.PROG(I1)
  381. ENDDO
  382. ENDDO
  383. ENDIF
  384. C
  385. C********** If we are here, we have the coeff. in the 'FACE' frame.
  386. C
  387. ENDDO
  388. SEGDES MELFP1
  389. ENDDO
  390. C
  391. NBTPOI=MLEFSC.INDEX(NBL+1)-1
  392. SEGADJ MLEFSC
  393. N1=NVMIN
  394. N2=NBTPOI
  395. SEGADJ MACOE1
  396. C
  397. SEGSUP MLEFP
  398. SEGDES MLEFSC
  399. SEGDES MACOE1
  400. SEGSUP MATA
  401. SEGSUP MATU
  402. SEGSUP MATV
  403. SEGSUP MLRSIG
  404. SEGSUP MLRTRA
  405. C
  406. SEGDES MELFL
  407. SEGSUP MTAA
  408. SEGSUP MINTAA
  409. SEGSUP MLRIN1
  410. SEGSUP MLRIN2
  411. SEGSUP MLRIN3
  412. C
  413. 9999 RETURN
  414. END
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  

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