Télécharger rlexca.eso

Retour à la liste

Numérotation des lignes :

  1. C RLEXCA SOURCE CHAT 05/01/13 03:02:11 5004
  2. SUBROUTINE RLEXCA(MELLIM,MLELSC,MLESBC,MLRDIS,
  3. & MLESCF,MATCOE)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : RLEXCA
  9. C
  10. C DESCRIPTION : Appellé par GRADIA
  11. C
  12. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  13. C
  14. C AUTEUR : A. BECCANTINI
  15. C
  16. C************************************************************************
  17. C
  18. C Inputs:
  19. C
  20. C MELLIM : MELEME spg boundary conditions
  21. C
  22. C MLELSC : neighbors of SOMMETS
  23. C
  24. C MLESBC : neighbors of SOMMETS belonging to the boundary
  25. C
  26. C MLRDIS : list of distances of a boundary SOMMET and its neighbors
  27. C
  28. C Output :
  29. C
  30. C MLESCF : final list of SOMMET points and their CENTRE neighbors
  31. C
  32. C MATCOE : coefficient to compute the reconstruction
  33. C
  34. C
  35. C**** Variables de COOPTIO
  36. C
  37. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  38. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  39. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  40. C & ,IECHO, IIMPI, IOSPI
  41. C & ,IDIM
  42. CC & ,MCOORD
  43. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  44. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  45. C & ,NORINC,NORVAL,NORIND,NORVAD
  46. C & ,NUCROU, IPSAUV, IFICLE, IPREFI
  47. CC
  48. IMPLICIT INTEGER(I-N)
  49. -INC CCOPTIO
  50. -INC SMCOORD
  51. -INC SMELEME
  52. POINTEUR MELLIM.MELEME
  53. C
  54. INTEGER NBL, NBTPOI
  55. SEGMENT MLELEM
  56. INTEGER INDEX(NBL+1)
  57. INTEGER LESPOI(NBTPOI)
  58. ENDSEGMENT
  59. POINTEUR MLELSC.MLELEM, MLESBC.MLELEM, MLESCF.MLELEM
  60. C
  61. INTEGER JG
  62. -INC SMLREEL
  63. POINTEUR MLRDIS.MLREEL, MLRSIG.MLREEL, MLRTRA.MLREEL
  64. -INC SMLENTI
  65. POINTEUR MLELIM.MLENTI, MLESOM.MLENTI
  66. C
  67. INTEGER N1,N2
  68. SEGMENT MATRIX
  69. REAL*8 MAT(N1,N2)
  70. ENDSEGMENT
  71. POINTEUR MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
  72. & ,MATCOE.MATRIX
  73. & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
  74. & ,MLRIN3.MLREEL
  75. C
  76. INTEGER NTP,NBTPO1,NBL1,ISOMM,IPOS,NGS,NLS,IPOS1,NVMAX,IPOS2
  77. & ,IPOS3,NVOI,NVMIN,NSOMM,I1,I2,I3,I4,NGC,IERSVD,NGS1,NVOIS1
  78. & ,NLLIM,IPCOOR,IPVOI,ICEL,NVOI0,IERR0,IVOI
  79. Real*8 XC,YC,ZC,SMAX,SMIN,DIST2,DIST21,ERRTOL,XS,YS,ZS,ERRTO1
  80. PARAMETER(ERRTOL=1.0D-15,ERRTO1=1.0D-7)
  81. LOGICAL LOGSVD
  82. C
  83. NTP=MCOORD.XCOOR(/1)/(IDIM+1)
  84. C
  85. C**** Le MELEME support de CL
  86. C
  87. IF(MELLIM.NE.0) THEN
  88. CALL KRIPAD(MELLIM,MLELIM)
  89. C
  90. C******* En KRIPAD
  91. C SEGACT MELLIM
  92. C SEGINI MLELIM
  93. C
  94. SEGDES MELLIM
  95. ELSE
  96. JG=NTP
  97. SEGINI MLELIM
  98. ENDIF
  99. C
  100. SEGACT MLESBC
  101. NBL1=MLESBC.INDEX(/1)-1
  102. NBTPO1=MLESBC.LESPOI(/1)
  103. JG=NTP
  104. SEGINI MLESOM
  105. DO I1 = 1, NBL1, 1
  106. IPOS=MLESBC.INDEX(I1)
  107. NGS=MLESBC.LESPOI(IPOS)
  108. MLESOM.LECT(NGS)=I1
  109. ENDDO
  110. C
  111. SEGACT MLELSC
  112. NBL=MLELSC.INDEX(/1)-1
  113. NSOMM=NBL
  114. NBTPOI=MLELSC.LESPOI(/1)
  115. NBTPOI=NBTPOI+NBTPO1-NBL1
  116. C
  117. C**** MLESCF
  118. C Structure Sommet-Centres Finale
  119. C NBTPOI(MLESCF) <= NBTPOI (SEGADJ à la fin)
  120. C
  121. SEGINI MLESCF
  122. C
  123. C**** La matrice de coeff.
  124. C
  125. C MATCOE(*,I) = coefficients concernants MLESCF.LESPOI(I)
  126. C
  127. C NVMIN : nombre minimum de voisins que un sommet
  128. C doit avoir pour une reconstruction lineaire
  129. C exacte
  130. C N2(MATCOE) <= N2 (SEGADJ à la fin)
  131. C
  132. NVMIN=IDIM+1
  133. N1=NVMIN
  134. N2=NBTPOI
  135. SEGINI MATCOE
  136. SEGACT MLRDIS
  137. C
  138. C**** Definition des matrices pour calculer MATCOE
  139. C Les dimensions
  140. C
  141. IPOS1=MLELSC.INDEX(1)
  142. NVMAX=0
  143. DO ISOMM = 1, NSOMM, 1
  144. IPOS=IPOS1
  145. IPOS1=MLELSC.INDEX(ISOMM+1)
  146. NVOI=IPOS1-IPOS-1
  147. NGS=MLELSC.LESPOI(IPOS)
  148. NLS=MLESOM.LECT(NGS)
  149. IF(NLS .GT. 0)THEN
  150. IPOS2=MLESBC.INDEX(NLS)
  151. IPOS3=MLESBC.INDEX(NLS+1)
  152. NVOI=NVOI+IPOS3-IPOS2-1
  153. ENDIF
  154. NVMAX=MAX(NVMAX,NVOI)
  155. ENDDO
  156. C
  157. C**** MATA = matrice à inverser (NVMAX,IDIM+1)
  158. C MATU = matice des vectors propres singuliers droite de MATA
  159. C (NVMAX,IDIM+1)
  160. C MATV = matrice des vectors propres singuliers gauche de MATA
  161. C (IDIM+1,IDIM+1)
  162. C Mais en INVSVD a dimension NVMAX,IDIM+1
  163. C MLRSIG =liste des valeurs singuliers de MATA
  164. C
  165. C N.B. MATA = MATU MATSIG t(MATV)
  166. C Si MATA non singulier
  167. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  168. C
  169. N1=NVMAX
  170. N2=NVMIN
  171. SEGINI MATA
  172. SEGINI MATU
  173. SEGINI MATV
  174. JG=NVMIN
  175. SEGINI MLRSIG
  176. SEGINI MLRTRA
  177. C MLRTRA segment de travail de la subroutine invsvd
  178. C
  179. C
  180. C MTAA : [t(MATA).MATA]
  181. C MINTAA : [inve(t(MATA) MATA)]
  182. C MLRIN1,2,3 : "temporary vectors"
  183. C
  184. N1=NVMIN
  185. N2=NVMIN
  186. SEGINI MTAA
  187. SEGINI MINTAA
  188. JG=NVMIN
  189. SEGINI MLRIN1
  190. SEGINI MLRIN2
  191. SEGINI MLRIN3
  192. C
  193. MLESCF.INDEX(1)=1
  194. NBTPOI=0
  195. DO ISOMM = 1, NSOMM, 1
  196. IPOS=MLESCF.INDEX(ISOMM)
  197. IPOS1=MLELSC.INDEX(ISOMM)
  198. NGS=MLELSC.LESPOI(IPOS1)
  199. IPCOOR=(NGS-1)*(IDIM+1)+1
  200. XS=MCOORD.XCOOR(IPCOOR)
  201. YS=MCOORD.XCOOR(IPCOOR+1)
  202. IF(IDIM .EQ. 3) ZS=MCOORD.XCOOR(IPCOOR+2)
  203. MLESCF.LESPOI(IPOS)=NGS
  204. NLLIM=MLELIM.LECT(NGS)
  205. IF(NLLIM.GT.0)THEN
  206. C
  207. C********** NGS est un noeud du MELEME des conditions
  208. C limites
  209. C
  210. MLESCF.INDEX(ISOMM+1)=IPOS+1
  211. MATCOE.MAT(1,IPOS)=1.0D0
  212. DO I2 = 2, NVMIN, 1
  213. MATCOE.MAT(I2,IPOS)=0.0D0
  214. ENDDO
  215. NBTPOI=NBTPOI+1
  216. C
  217. C********** On passe au sommet suivant (i.e. ISOMM->ISOMM+1)
  218. C
  219. ELSE
  220. DO I2 = 1, NVMIN, 1
  221. MATCOE.MAT(I2,IPOS)=0.0D0
  222. ENDDO
  223. NLS=MLESOM.LECT(NGS)
  224. NVOI=MLELSC.INDEX(ISOMM+1)-IPOS1-1
  225. NVOI0=NVOI
  226. C NVOI0 = nombre de centre voisins directe!
  227. DO I2=1,NVOI,1
  228. NGC=MLELSC.LESPOI(IPOS1+I2)
  229. MLESCF.LESPOI(IPOS+I2)=NGC
  230. IPCOOR=(NGC-1)*(IDIM+1)+1
  231. XC=MCOORD.XCOOR(IPCOOR)
  232. YC=MCOORD.XCOOR(IPCOOR+1)
  233. MATA.MAT(I2,1)=1
  234. MATA.MAT(I2,2)=XC-XS
  235. MATA.MAT(I2,3)=YC-YS
  236. IF(IDIM.EQ.3)THEN
  237. ZC=MCOORD.XCOOR(IPCOOR+2)
  238. MATA.MAT(I2,4)=ZC-ZS
  239. ENDIF
  240. ENDDO
  241. IF(NLS.EQ.0)THEN
  242. C
  243. C************* NGS est un noeud À l'interieur du domaine
  244. C
  245. LOGSVD=.TRUE.
  246. CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  247. & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
  248. & MLRTRA.PROG)
  249. IF(IERSVD.NE.0)THEN
  250. C INVSVD did not work
  251. LOGSVD=.FALSE.
  252. ELSE
  253. SMAX=0.0D0
  254. DO I2=1,NVMIN,1
  255. SMAX=MAX(SMAX,MLRSIG.PROG(I2))
  256. ENDDO
  257. SMIN=SMAX
  258. DO I2=1,NVMIN,1
  259. SMIN=MIN(SMIN,MLRSIG.PROG(I2))
  260. ENDDO
  261. ENDIF
  262. IF((SMIN/SMAX) .LT. ERRTOL)THEN
  263. C INVSVD did not work
  264. LOGSVD=.FALSE.
  265. ENDIF
  266. ELSE
  267. C
  268. C************* NGS est un noeud sur le bord mais pas un noeud des
  269. C conditions limites
  270. C
  271. IPVOI=MLESBC.INDEX(NLS)
  272. NGS1=MLESBC.LESPOI(IPVOI)
  273. NVOIS1=MLESBC.INDEX(NLS+1)-IPVOI-1
  274. IF((NVOIS1.LE.0).OR.(NGS .NE. NGS1))THEN
  275. WRITE(IOIMP,*) 'Subroutine rlexca.eso.'
  276. CALL ERREUR(5)
  277. GOTO 9999
  278. ENDIF
  279. DIST2=MLRDIS.PROG(IPVOI+1)*1.10D0
  280. ICEL=1
  281. C
  282. C************* Repeat until
  283. C
  284. 10 CONTINUE
  285. NGC=MLESBC.LESPOI(IPVOI+ICEL)
  286. DIST21=MLRDIS.PROG(IPVOI+ICEL)
  287. IF((DIST21.LT.DIST2) .OR. (NVOI .LT. NVMIN))THEN
  288. NVOI=NVOI+1
  289. MLESCF.LESPOI(IPOS+NVOI)=NGC
  290. IPCOOR=(NGC-1)*(IDIM+1)+1
  291. XC=MCOORD.XCOOR(IPCOOR)
  292. YC=MCOORD.XCOOR(IPCOOR+1)
  293. MATA.MAT(NVOI,1)=1
  294. MATA.MAT(NVOI,2)=XC-XS
  295. MATA.MAT(NVOI,3)=YC-YS
  296. IF(IDIM.EQ.3)THEN
  297. ZC=MCOORD.XCOOR(IPCOOR+2)
  298. MATA.MAT(NVOI,4)=ZC-ZS
  299. ENDIF
  300. ICEL=ICEL+1
  301. DIST2=DIST21*1.10D0
  302. IF(ICEL .LE. NVOIS1) GOTO 10
  303. ENDIF
  304. LOGSVD=.TRUE.
  305. CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  306. & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT
  307. & ,IERSVD,MLRTRA.PROG)
  308. IF(IERSVD.NE.0)THEN
  309. SMIN=0.0D0
  310. SMAX=1.0D0
  311. ELSE
  312. SMAX=0.0D0
  313. DO I2=1,NVMIN,1
  314. SMAX=MAX(SMAX,MLRSIG.PROG(I2))
  315. ENDDO
  316. SMIN=SMAX
  317. DO I2=1,NVMIN,1
  318. SMIN=MIN(SMIN,MLRSIG.PROG(I2))
  319. ENDDO
  320. ENDIF
  321. IF((SMIN/SMAX) .LT. ERRTO1)THEN
  322. C
  323. C**************** We add other neighbors
  324. C
  325. ICEL=NVOI-NVOI0
  326. IF(ICEL.LT.NVOIS1)THEN
  327. ICEL=ICEL+1
  328. NGC=MLESBC.LESPOI(IPVOI+ICEL)
  329. DIST2=MLRDIS.PROG(IPVOI+ICEL)*1.10D0
  330. GOTO 10
  331. ELSE
  332. IF((SMIN/SMAX) .LT. (ERRTOL))THEN
  333. LOGSVD=.FALSE.
  334. ENDIF
  335. ENDIF
  336. ENDIF
  337. C
  338. C************* Fin repeat until
  339. C
  340. ENDIF
  341. C
  342. C********** We have two different cases:
  343. C LOGSVD = .TRUE. -> we have the coeff. we are looking for
  344. C LOGSVD = .FALSE. -> we have not
  345. C
  346. IF(LOGSVD)THEN
  347. DO I4=1,NVMIN,1
  348. DO I3=1,NVOI,1
  349. IPVOI=IPOS+I3
  350. DO I2=1,NVMIN,1
  351. MATCOE.MAT(I2,IPVOI)=MATCOE.MAT(I2,IPVOI)+
  352. & (MATV.MAT(I2,I4)*MATU.MAT(I3,I4)/
  353. & MLRSIG.PROG(I4))
  354. ENDDO
  355. ENDDO
  356. ENDDO
  357. ELSE
  358. WRITE (IOIMP,*) 'rlexca.eso'
  359. C 22 0
  360. C Opération malvenue. Résultat douteux
  361. CALL ERREUR(22)
  362. C
  363. C************* INVSVD does not worked
  364. C For each neighbor k we have to compute the solution
  365. C of
  366. C
  367. C t(MATA) MATA x = t(MATA) * b
  368. C
  369. C where b= \sum_l e_l \delta(k,l) = e_k
  370. C
  371. C To do that, we compute
  372. C
  373. C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  374. C
  375. C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  376. C [t(MATA) MATA] X_0]
  377. C
  378. C
  379. C************* We compute [t(MATA) MATA]
  380. C We store it in the upper triangle of MTAA(NVMIN,NVMIN)
  381. C
  382. DO I1=1,NVMIN,1
  383. DO I2=I1,NVMIN,1
  384. MTAA.MAT(I1,I2)=0.0D0
  385. ENDDO
  386. ENDDO
  387. C
  388. DO I1=1,NVMIN,1
  389. DO I2=I1,NVMIN,1
  390. DO I3=1,NVOI,1
  391. MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
  392. & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
  393. ENDDO
  394. ENDDO
  395. ENDDO
  396. C
  397. C************* We compute [inve(t(MATA) MATA)]
  398. C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
  399. C
  400. CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
  401. & MLRIN2.PROG,IERR0)
  402. IF(IERR0.NE.0)THEN
  403. WRITE(IOIMP,*) 'subroutine rlexca.eso.'
  404. C 26 2
  405. C Tache impossible. Probablement données erronées
  406. CALL ERREUR(26)
  407. GOTO 9999
  408. ENDIF
  409. C
  410. C************* We complete MTAA and MINTAA
  411. C
  412. DO I1=1,NVMIN,1
  413. DO I2=I1+1,NVMIN,1
  414. MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
  415. MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
  416. ENDDO
  417. ENDDO
  418. C
  419. DO IVOI=1,NVOI,1
  420. C
  421. C**************** We compute [t(MATA) . b] and we store it in MLRIN1.PROG
  422. C
  423. DO I1=1,NVMIN,1
  424. MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)
  425. MLRIN2.PROG(I1)=0.0D0
  426. MLRIN3.PROG(I1)=0.0D0
  427. ENDDO
  428. C
  429. C**************** X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  430. C X_0(i1) into MLRIN2.PROG(I1)
  431. C
  432. DO I2=1,NVMIN,1
  433. DO I1=1,NVMIN,1
  434. MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
  435. & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
  436. ENDDO
  437. ENDDO
  438. C
  439. C**************** X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  440. C [t(MATA) MATA] X_0]
  441. C
  442. C [t(MATA) MATA] X_0 into MLRIN3.PROG
  443. C
  444. DO I2=1,NVMIN,1
  445. DO I1=1,NVMIN,1
  446. MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
  447. & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
  448. ENDDO
  449. ENDDO
  450. C
  451. C**************** Now we have
  452. C [t(MATA) . b] in MLRIN1.PROG
  453. C X_0(i1) in MLRIN2.PROG(I1)
  454. C [t(MATA) MATA] X_0 in MLRIN3.PROG
  455. C
  456. C X_1(i1) in MLRCOE.MAT(i1,IVOI)
  457. C
  458. IPVOI=IPOS+IVOI
  459. DO I1=1,NVMIN,1
  460. DO I2=1,NVMIN,1
  461. MATCOE.MAT(I1,IPVOI)=MATCOE.MAT(I1,IPVOI)+
  462. & (MINTAA.MAT(I1,I2)*
  463. & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
  464. ENDDO
  465. MATCOE.MAT(I1,IPVOI)=MATCOE.MAT(I1,IPVOI)+
  466. & MLRIN2.PROG(I1)
  467. ENDDO
  468. ENDDO
  469. ENDIF
  470. CC TEST
  471. C if(isomm .eq. (nsomm - 1))then
  472. C write(*,*) 'ngs =', NGS
  473. C write(*,*) 'invide',LOGSVD
  474. C write(*,*) 'coeff(1) =',(matcoe.mat(1,ipos+ivoi),ivoi=1
  475. C $ ,nvoi,1)
  476. C write(*,*) 'coeff(2) =',(matcoe.mat(2,ipos+ivoi),ivoi=1
  477. C $ ,nvoi,1)
  478. C write(*,*) 'coeff(3) =',(matcoe.mat(3,ipos+ivoi),ivoi=1
  479. C $ ,nvoi,1)
  480. C if(idim.eq.3) write(*,*) 'coeff(3)='
  481. C & ,(matcoe.mat(4,ipos+ivoi),ivoi=1,nvoi,1)
  482. C dist2=0.0D0
  483. C do ivoi=1,nvoi,1
  484. C dist2=dist2+matcoe.mat(1,ipos+ivoi)
  485. C enddo
  486. C write(*,*) 'sum=',dist2
  487. CC
  488. C call erreur(21)
  489. C goto 9999
  490. C endif
  491. C
  492. C********** If we are here, we have the coeff. in the 'SOMMET' frame.
  493. C
  494. C
  495. NBTPOI=NBTPOI+NVOI+1
  496. MLESCF.INDEX(ISOMM+1)=IPOS+NVOI+1
  497. ENDIF
  498. ENDDO
  499. C
  500. SEGADJ MLESCF
  501. SEGDES MLESCF
  502. C
  503. N2=NBTPOI
  504. SEGADJ MATCOE
  505. SEGDES MATCOE
  506. C
  507. IF(MELLIM .NE. 0) SEGDES MELLIM
  508. SEGSUP MLELIM
  509. SEGSUP MLESBC
  510. SEGSUP MLESOM
  511. SEGSUP MLELSC
  512. SEGSUP MLRDIS
  513. SEGSUP MATA
  514. SEGSUP MATU
  515. SEGSUP MATV
  516. SEGSUP MLRSIG
  517. SEGSUP MLRTRA
  518. C
  519. SEGSUP MTAA
  520. SEGSUP MINTAA
  521. SEGSUP MLRIN1
  522. SEGSUP MLRIN2
  523. SEGSUP MLRIN3
  524. C
  525. 9999 RETURN
  526. END
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  

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