Télécharger rleca1.eso

Retour à la liste

Numérotation des lignes :

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

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