Télécharger rlexca.eso

Retour à la liste

Numérotation des lignes :

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

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