Télécharger rlencf.eso

Retour à la liste

Numérotation des lignes :

  1. C RLENCF SOURCE CHAT 05/01/13 03:01:14 5004
  2. SUBROUTINE RLENCF(MELFL,MELFP,MLEPOI,MLECOE)
  3. C************************************************************************
  4. C
  5. C PROJET : CASTEM 2000
  6. C
  7. C NOM : RLENCF
  8. C
  9. C DESCRIPTION : Appelle par GRADI2
  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 MLESCF : list of SOMMET points and their CENTRE neighbors
  20. C
  21. C MATCOE : coeff. for linear exact reconstruction of MLESCF
  22. C
  23. C MLEFSC : list of FACES points and their neighbors (CENTRE and SOMMET
  24. C points)
  25. C
  26. C MACOE1 : coeff. for linear exact reconstruction of MLEFSC
  27. C
  28. C Output
  29. C
  30. C MLEFC : list of FACES points and their neighbors (CENTRE points only)
  31. C
  32. C MACOE2 : coeff. for linear exact reconstruction of MLEFC
  33. C
  34. C This subroutine computes the coefficients to compute the gradient
  35. C at intefaces with respect to the values on its neighbors.
  36. C The neighbors are 'CENTRE' points (in FACEL) ore 'SOMMET' points
  37. C (in FACEP). which allow to perform
  38. C A linear exact reconstruction is performed via a least square
  39. C method.
  40. C
  41. C**** Inputs:
  42. C
  43. C MELFL = the 'FACEL' meleme
  44. C MELFP = the 'FACEP' meleme
  45. C
  46. C**** Output:
  47. C
  48. C MLEPOI = list of integers.
  49. C MLEPOI.LECT(i) points to the list neighbors of
  50. C MELFL.NUM(2,i)
  51. C MLECOE = list of integers.
  52. C MLECOE.LECT(i) points to the matrix of coeffients to
  53. C compute the gradient
  54. C
  55. C**** Variables de COOPTIO
  56. C
  57. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  58. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  59. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  60. C & ,IECHO, IIMPI, IOSPI
  61. C & ,IDIM
  62. CC & ,MCOORD
  63. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  64. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  65. C & ,NORINC,NORVAL,NORIND,NORVAD
  66. C & ,NUCROU, IPSAUV, IFICLE, IPREFI
  67. C
  68. IMPLICIT INTEGER(I-N)
  69. INTEGER NBSOUS,ISOUS,NBELEM,NBNO,NVMAX,NVMIN,NFAC,IFAC,IELEM
  70. & ,NGF,NGF1,NGC1,NGC2,INOEU,NGS,IPCOOR,NVOI,IVOI,NGV
  71. & ,IERSVD,IERR0
  72. & ,I1,I2,I3,I4
  73.  
  74. REAL*8 XF,YF,ZF, XV, YV, ZV,SMAX,SMIN,ERRTOL
  75. PARAMETER(ERRTOL=1.0D-16)
  76. LOGICAL LOGSVD
  77. -INC CCOPTIO
  78. -INC SMCOORD
  79. -INC SMELEME
  80. -INC SMLREEL
  81. -INC SMLENTI
  82. -INC SMCHPOI
  83. INTEGER JG
  84. INTEGER N1,N2
  85. SEGMENT MATRIX
  86. REAL*8 MAT(N1,N2)
  87. ENDSEGMENT
  88. POINTEUR MELFL.MELEME, MELFP.MELEME, MLEPOI.MLENTI,MLECOE.MLENTI
  89. & , MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
  90. & ,MATCOE.MATRIX,MLRSIG.MLREEL,MLRTRA.MLREEL
  91. & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
  92. & ,MLRIN3.MLREEL,MLEFP.MLENTI,MELFP1.MELEME
  93. C
  94. SEGACT MELFP
  95. C
  96. C
  97. C**** En 2D MELFP est un maillage elementaire
  98. C En 3D pas à priori
  99. C -> MLEFP contains the list of LISOUS
  100. C
  101. NBSOUS=MELFP.LISOUS(/1)
  102. C NBSOUS=0 fais un peux chier!
  103. JG=MAX(NBSOUS,1)
  104. SEGINI MLEFP
  105. IF(NBSOUS .EQ. 0)THEN
  106. MLEFP.LECT(1)=MELFP
  107. ELSE
  108. DO ISOUS=1,NBSOUS,1
  109. MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS)
  110. ENDDO
  111. ENDIF
  112. SEGDES MELFP
  113. SEGACT MELFL
  114. NFAC=MELFL.NUM(/2)
  115. JG=NFAC
  116. SEGINI MLEPOI
  117. SEGINI MLECOE
  118. C
  119. C**** Loop on the faces to compute NVMAX
  120. C (maximum number of neighbors)
  121. C
  122. NBSOUS=MLEFP.LECT(/1)
  123. NVMAX=0
  124. DO ISOUS = 1, NBSOUS, 1
  125. MELFP1=MLEFP.LECT(ISOUS)
  126. SEGACT MELFP1
  127. NBELEM=MELFP1.NUM(/2)
  128. NBNO=MELFP1.NUM(/1) - 1
  129. C
  130. C The ISOUS elementary meleme of 'FACEP' has NBELEM elements
  131. C which contain NBNO sommets and one point face. Each face has
  132. C NBNO 'SOMMET' neighbors and 2 'CENTRE' neighbors.
  133. C
  134. NVMAX=MAX(NVMAX,NBNO+2)
  135. ENDDO
  136. C
  137. C MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1)
  138. C MATU = matrix of the singular right eigenvectors of MATA
  139. C (NVMAX,IDIM+1)
  140. C MATV = matrix of the singular left eigenvectors of MATA
  141. C (IDIM+1,IDIM+1)
  142. C But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1
  143. C MLRSIG = singular values of MATA (IDIM+1)
  144. C
  145. C N.B. MATA = MATU MATSIG t(MATV)
  146. C If MATA is non singular,
  147. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  148. C
  149. C MLRTRA temporary vector of invsvd.eso
  150. C NVMIN = IDIM + 1 (the most little dimension of the matrices)
  151. C
  152. N1=NVMAX
  153. N2=IDIM+1
  154. SEGINI MATA
  155. SEGINI MATU
  156. SEGINI MATV
  157. JG=IDIM+1
  158. SEGINI MLRSIG
  159. SEGINI MLRTRA
  160. NVMIN=N2
  161. C
  162. C MTAA : [t(MATA).MATA]
  163. C MINTAA : [inve(t(MATA) MATA)]
  164. C MLRIN1,2,3 : "temporary vectors"
  165. C
  166. N1=NVMIN
  167. N2=NVMIN
  168. SEGINI MTAA
  169. SEGINI MINTAA
  170. JG=NVMIN
  171. SEGINI MLRIN1
  172. SEGINI MLRIN2
  173. SEGINI MLRIN3
  174. C
  175. C**** Loop on faces to compute the coefficient
  176. C
  177. NBSOUS=MLEFP.LECT(/1)
  178. IFAC=0
  179. DO ISOUS = 1, NBSOUS, 1
  180. MELFP1=MLEFP.LECT(ISOUS)
  181. NBELEM=MELFP1.NUM(/2)
  182. NBNO=MELFP1.NUM(/1) - 1
  183. DO IELEM=1,NBELEM,1
  184. IFAC=IFAC+1
  185. NGF=MELFP1.NUM(NBNO+1,IELEM)
  186. NGF1=MELFL.NUM(2,IFAC)
  187. NGC1=MELFL.NUM(1,IFAC)
  188. NGC2=MELFL.NUM(3,IFAC)
  189. IF(NGF .NE. NGF1)THEN
  190. WRITE(IOIMP,*) 'Subroutine rlencf.eso'
  191. CALL ERREUR(5)
  192. GOTO 9999
  193. ENDIF
  194. IF(NGC1 .EQ. NGC2)THEN
  195. JG=NBNO+1
  196. ELSE
  197. JG=NBNO+2
  198. ENDIF
  199. C
  200. C********** We create the list of neighbors.
  201. C
  202. SEGINI MLENTI
  203. MLEPOI.LECT(IFAC)=MLENTI
  204. C
  205. C********** We create the matrix of coefficients.
  206. C
  207. N2=JG
  208. N1=NVMIN
  209. SEGINI MATCOE
  210. MLECOE.LECT(IFAC)=MATCOE
  211. C
  212. C********** We fill the list of neighbors.
  213. C
  214. DO INOEU=1,NBNO,1
  215. NGS=MELFP1.NUM(INOEU,IELEM)
  216. MLENTI.LECT(INOEU)=NGS
  217. ENDDO
  218. MLENTI.LECT(NBNO+1)=NGC1
  219. IF(NGC1 .NE. NGC2) MLENTI.LECT(NBNO+2)=NGC2
  220. ENDDO
  221. SEGDES MELFP1
  222. ENDDO
  223. C
  224. C**** Test for MLEPOI
  225. C
  226. C do ifac=1,nfac,1
  227. C mlenti=mlepoi.lect(ifac)
  228. C nbno=mlenti.lect(/1)
  229. C write(*,*) 'ngf=',melfl.num(2,ifac)
  230. C write(*,*) 'nvoi=',nbno
  231. C write(*,*) (mlenti.lect(inoeu),inoeu=1,nbno,1)
  232. C enddo
  233. C
  234. C**** We have to fill the matrix coefficient
  235. C
  236. NFAC=MELFL.NUM(/2)
  237. DO IFAC=1,NFAC,1
  238. NGF=MELFL.NUM(2,IFAC)
  239. IPCOOR=((NGF-1)*(IDIM+1))+1
  240. XF=MCOORD.XCOOR(IPCOOR)
  241. YF=MCOORD.XCOOR(IPCOOR+1)
  242. IF(IDIM .EQ. 3) ZF=MCOORD.XCOOR(IPCOOR+2)
  243. C
  244. C******* The neighbors of NGF
  245. C
  246. MLENTI=MLEPOI.LECT(IFAC)
  247. NVOI=MLENTI.LECT(/1)
  248. C
  249. C******* The matrix where we store the coefficients
  250. C
  251. MATCOE=MLECOE.LECT(IFAC)
  252. C
  253. C******* Loop involving the neighbors
  254. C We create the matrix to "pseudoinvert"
  255. C
  256. DO IVOI=1,NVOI,1
  257. NGV=MLENTI.LECT(IVOI)
  258. IPCOOR=((NGV-1)*(IDIM+1))+1
  259. XV=MCOORD.XCOOR(IPCOOR)
  260. YV=MCOORD.XCOOR(IPCOOR+1)
  261. MATA.MAT(IVOI,1)=1
  262. MATA.MAT(IVOI,2)=XV-XF
  263. MATA.MAT(IVOI,3)=YV-YF
  264. IF(IDIM.EQ.3)THEN
  265. ZV=MCOORD.XCOOR(IPCOOR+2)
  266. MATA.MAT(IVOI,4)=ZV-ZF
  267. ENDIF
  268. ENDDO
  269. C
  270. C********** Now we have to invert this matrix
  271. C
  272. LOGSVD=.TRUE.
  273. CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  274. & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
  275. & MLRTRA.PROG)
  276. IF(IERSVD.NE.0)THEN
  277. C
  278. C************* SVD decomposition of the matrix does not work
  279. C
  280. LOGSVD=.FALSE.
  281. ELSE
  282. C
  283. C************ We check the condition number of MATA
  284. C
  285. SMAX=0.0D0
  286. DO I1=1,NVMIN,1
  287. SMAX=MAX(SMAX,MLRSIG.PROG(I1))
  288. ENDDO
  289. SMIN=SMAX
  290. DO I1=1,NVMIN,1
  291. SMIN=MIN(SMIN,MLRSIG.PROG(I1))
  292. ENDDO
  293. IF((SMIN/SMAX) .LT. ERRTOL)THEN
  294. LOGSVD=.FALSE.
  295. ENDIF
  296. ENDIF
  297. C
  298. CC TEST
  299. C write(*,*) 'LOGSVD=.FALSE.'
  300. C LOGSVD=.FALSE.
  301. C
  302. IF(LOGSVD)THEN
  303. C
  304. C********** INVSVD worked
  305. C
  306. C MATA = MATU MATSIG t(MATV)
  307. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  308. C
  309. DO I4=1,NVMIN,1
  310. DO IVOI=1,NVOI,1
  311. DO I2=1,IDIM+1,1
  312. C I2=1 is the only coefficient we are not interested
  313. C in. But we computed it to verify that
  314. C sum_ivoi MATCOE.MAT(ivoi,1) = 1
  315. C
  316. MATCOE.MAT(I2,IVOI)=MATCOE.MAT(I2,IVOI)+
  317. & (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4)
  318. & /MLRSIG.PROG(I4))
  319. ENDDO
  320. ENDDO
  321. ENDDO
  322. ELSE
  323. WRITE (IOIMP,*) 'rlencf.eso'
  324. C 22 0
  325. C Opération malvenue. Résultat douteux
  326. CALL ERREUR(22)
  327. C
  328. C********** INVSVD does not worked
  329. C For each neighbor k we have to compute the solution
  330. C of
  331. C
  332. C t(MATA) MATA x = t(MATA) * b
  333. C
  334. C where b= \sum_l e_l \delta(k,l) = e_k
  335. C
  336. C To do that, we compute
  337. C
  338. C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  339. C
  340. C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  341. C [t(MATA) MATA] X_0]
  342. C
  343. C
  344. C********** We compute [t(MATA) MATA]
  345. C We store it in the upper triangle of MTAA(NVMIN,NVMIN)
  346. C
  347. DO I1=1,NVMIN,1
  348. DO I2=I1,NVMIN,1
  349. MTAA.MAT(I1,I2)=0.0D0
  350. ENDDO
  351. ENDDO
  352. C
  353. DO I1=1,NVMIN,1
  354. DO I2=I1,NVMIN,1
  355. DO I3=1,NVOI,1
  356. MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
  357. & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
  358. ENDDO
  359. ENDDO
  360. ENDDO
  361. C
  362. C********** We compute [inve(t(MATA) MATA)]
  363. C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
  364. C
  365. CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
  366. & MLRIN2.PROG,IERR0)
  367. IF(IERR0.NE.0)THEN
  368. WRITE(IOIMP,*) 'subroutine rlencf.eso.'
  369. C 26 2
  370. C Tache impossible. Probablement données erronées
  371. CALL ERREUR(26)
  372. GOTO 9999
  373. ENDIF
  374. C
  375. C********** We complete MTAA and MINTAA
  376. C
  377. DO I1=1,NVMIN,1
  378. DO I2=I1+1,NVMIN,1
  379. MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
  380. MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
  381. ENDDO
  382. ENDDO
  383. C
  384. DO IVOI=1,NVOI,1
  385. C
  386. C************* We compute [t(MATA) . b] and we store it in MLRIN1.PROG
  387. C
  388. DO I1=1,NVMIN,1
  389. MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)
  390. MLRIN2.PROG(I1)=0.0D0
  391. MLRIN3.PROG(I1)=0.0D0
  392. ENDDO
  393. C
  394. C************* X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  395. C X_0(i1) into MLRIN2.PROG(I1)
  396. C
  397. DO I2=1,NVMIN,1
  398. DO I1=1,NVMIN,1
  399. MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
  400. & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
  401. ENDDO
  402. ENDDO
  403. C
  404. C************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  405. C [t(MATA) MATA] X_0]
  406. C
  407. C [t(MATA) MATA] X_0 into MLRIN3.PROG
  408. C
  409. DO I2=1,NVMIN,1
  410. DO I1=1,NVMIN,1
  411. MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
  412. & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
  413. ENDDO
  414. ENDDO
  415. C
  416. C************* Now we have
  417. C [t(MATA) . b] in MLRIN1.PROG
  418. C X_0(i1) in MLRIN2.PROG(I1)
  419. C [t(MATA) MATA] X_0 in MLRIN3.PROG
  420. C
  421. C X_1(i1) in MATCOE.MAT(IVOI,i1)
  422. C
  423. DO I1=1,IDIM+1,1
  424. C The only unuseful one is I1=1
  425. DO I2=1,NVMIN,1
  426. MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+
  427. & (MINTAA.MAT(I1,I2)*
  428. & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
  429. ENDDO
  430. MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+
  431. & MLRIN2.PROG(I1)
  432. ENDDO
  433. ENDDO
  434. ENDIF
  435. CC
  436. CC TEST
  437. C write(*,*) 'ngf =', NGF
  438. C write(*,*) 'invide',LOGSVD
  439. C write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1)
  440. C write(*,*) 'coeff(1) =',(matcoe.mat(1,ivoi),ivoi=1,nvoi,1)
  441. C write(*,*) 'coeff(2) =',(matcoe.mat(2,ivoi),ivoi=1,nvoi,1)
  442. C write(*,*) 'coeff(3) =',(matcoe.mat(3,ivoi),ivoi=1,nvoi,1)
  443. C if(idim.eq.3) write(*,*) 'coeff(4)=',
  444. C & (matcoe.mat(4,ivoi),ivoi=1,nvoi,1)
  445. C xv=0.0D0
  446. C do ivoi=1,nvoi,1
  447. C xv=xv+matcoe.mat(1,ivoi)
  448. C enddo
  449. C write(*,*) 'sum=',xv
  450. CC
  451. MLENTI=MLEPOI.LECT(IFAC)
  452. SEGDES MLENTI
  453. MATCOE=MLECOE.LECT(IFAC)
  454. SEGDES MATCOE
  455. ENDDO
  456. C
  457. SEGDES MATCOE
  458. SEGDES MLEPOI
  459. C
  460. SEGDES MELFL
  461. SEGSUP MLEFP
  462. C
  463. SEGSUP MATA
  464. SEGSUP MATU
  465. SEGSUP MATV
  466. SEGSUP MLRSIG
  467. SEGSUP MLRTRA
  468. C
  469. SEGSUP MTAA
  470. SEGSUP MINTAA
  471. SEGSUP MLRIN1
  472. SEGSUP MLRIN2
  473. SEGSUP MLRIN3
  474. C
  475. 9999 CONTINUE
  476. RETURN
  477. END
  478.  
  479.  
  480.  
  481.  

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