Télécharger rlenco.eso

Retour à la liste

Numérotation des lignes :

  1. C RLENCO SOURCE CHAT 05/01/13 03:01:24 5004
  2. SUBROUTINE RLENCO(MELSOM,MELCEN,MELLI1,MELLI2,INORM,MLEPOI,MLECOE)
  3. C************************************************************************
  4. C
  5. C PROJET : CASTEM 2000
  6. C
  7. C NOM : RLENCO
  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 This subroutine computes the coefficients which allow to perform
  18. C a linear exact reconstruction of each 'sommet' (vertex) value
  19. C Indeed, given the i-th sommet (V=MELSOM.NUM(1,i)), given the list
  20. C of its neighbors MLEPOI.LECT(i), we have to compute the matrix of
  21. C coefficients a(i,j) such that
  22. C
  23. C VAL(MELSOM.NUM(1,i)) = \sum_{neig1} a(i,neig1) * VAL(neig1) +
  24. C \sum_{neig2} a(i,neig2) * VALNOR(neig2)
  25. C
  26. C where neig1 is a 'CENTRE' point or a 'Dirichlet boundary condition'
  27. C point, neig2 is a 'Von Neumann boundary condition' point, VAL(j)
  28. C is the value of the function in the j-th neighbor, VALNOR(j) is
  29. C the value of grad(VAL).n is the j-th neighbor.
  30. C
  31. C
  32. C To compute these coefficients, we impose that VAL is linear with
  33. C respect to the coordinates of V and its neighbor. Then, since
  34. C there are less unknowns than equations, we use a least square
  35. C method. If the neighbor is a 'CENTRE' point or a 'Dirichlet
  36. C boundary condition' point, we impose that
  37. C
  38. C A + B (x_{neig} - x_V) + C (y_{neig} - y_V) + D (z_{neig} - z_V)
  39. C = VAL(neig)
  40. C
  41. C If the neighbor is a 'von Neumann boundary condition' point,
  42. C we impose that
  43. C
  44. C |r_{neig} - r_V| (B nx_{neig} + C ny_{neig} + D nz_{neig}) =
  45. C VALNOR(neig) * |r_{neig} - r_V|
  46. C
  47. c r = (x,y,z)
  48. c n = (nx, ny, nz) is the normal to the surface
  49. C
  50. C To determine A, B, C, D we have to solve
  51. C
  52. C MATA . tran(A, B, C, D) = tran(VAL(1), VAL(2), ... , VALNOR(j) /
  53. C (r_{j} - r_V)
  54. C
  55. C with MATA(1,*) = (1,x_1-x_V,y_1-y_V,z_1-z_V)
  56. C MATA(2,*) = (1,x_2-x_V,y_2-y_V,z_2-z_V)
  57. C ...
  58. C MATA(j,*) = (0,nx_{j}, ny_{j}, nz_{j})
  59. C ...
  60. C
  61. C To determine a(i,neig) we have to solve a linear system (for each
  62. C neig). If the neighbor is a 'CENTRE' point or a 'Dirichlet boundary
  63. C condition' point we solve
  64. C
  65. C MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) = e_{neig}
  66. C
  67. C with e_{neig} = tran (0,0,...,0,1,0,...) (1 in the neig position)
  68. C
  69. C If the neighbor is a 'von Neumann boundary condition' points,
  70. C we solve
  71. C
  72. C MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) =
  73. C e_{neig} (r_{neig} - r_V)
  74. C
  75. C Since each linear system involves MATA, we have to "pseudoinvert"
  76. C it.
  77. C
  78. C**** Inputs:
  79. C
  80. C MELSOM = the 'SOMMET' meleme
  81. C MELCEN = the 'CENTRE' meleme
  82. C MELLI1 = the support of Dirichlet boundary conditions
  83. C MELLI2 = the support of von Neumann boundary conditions
  84. C INOMR = CHAMPOINT des normales aux faces
  85. C MLEPOI = list of integers.
  86. C MLEPOI.LECT(i) points to the list neighbors of
  87. C MELSOM.NUM(1,i)
  88. C
  89. C**** Output:
  90. C
  91. C MLECOE = list of integers.
  92. C
  93. C MLECOE.LECT(i) points to the MLREEL MLRCOE which
  94. C contain the a(i,neig)
  95. C
  96. C**** Variables de COOPTIO
  97. C
  98. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  99. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  100. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  101. C & ,IECHO, IIMPI, IOSPI
  102. C & ,IDIM
  103. CC & ,MCOORD
  104. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  105. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  106. C & ,NORINC,NORVAL,NORIND,NORVAD
  107. C & ,NUCROU, IPSAUV, IFICLE, IPREFI
  108. C
  109. IMPLICIT INTEGER(I-N)
  110. INTEGER INORM, NSOMM, INOEU, NGS, IPCOOR, IVOI, NVOI, NGV, NLV
  111. & ,NLV1,NVMAX,NVMIN,IERSVD,IERR0,NLN
  112. & ,I1,I2,I3,I4
  113. REAL*8 XS,YS,ZS,XV,YV,ZV,DX,ERRTOL,SMAX,SMIN
  114. PARAMETER(ERRTOL=1.0D-16)
  115. CHARACTER*8 TYPE
  116. LOGICAL LOGSVD
  117. -INC CCOPTIO
  118. -INC SMCOORD
  119. -INC SMELEME
  120. -INC SMLREEL
  121. -INC SMLENTI
  122. -INC SMCHPOI
  123. INTEGER JG
  124. INTEGER N1,N2
  125. SEGMENT MATRIX
  126. REAL*8 MAT(N1,N2)
  127. ENDSEGMENT
  128. POINTEUR MELSOM.MELEME, MLEPOI.MLENTI,MELLI1.MELEME,MELLI2.MELEME
  129. & ,MLELI1.MLENTI,MLELI2.MLENTI,MELCEN.MELEME,MLECEN.MLENTI
  130. & ,MPNORM.MPOVAL, MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
  131. & ,MLRCOE.MLREEL,MLRSIG.MLREEL,MLRTRA.MLREEL
  132. & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
  133. & ,MLRIN3.MLREEL, MLECOE.MLENTI, MLRB.MLREEL
  134. & ,MELFAC.MELEME,MLEFAC.MLENTI
  135. C
  136. C**** We create the MLENTI for the centers
  137. C
  138. CALL KRIPAD(MELCEN,MLECEN)
  139. IF(IERR .NE. 0)GOTO 9999
  140. C En KRIPAD
  141. C SEGACT MELCEN
  142. C SEGINI MLECEN
  143. C
  144. C**** We create the MLENTI for the BC
  145. C
  146. CALL KRIPAD(MELLI1,MLELI1)
  147. IF(IERR .NE. 0)GOTO 9999
  148. CALL KRIPAD(MELLI2,MLELI2)
  149. IF(IERR .NE. 0)GOTO 9999
  150. C SEGACT MELLI1
  151. C SEGINI MLELI1
  152. C SEGACT MELLI2
  153. C SEGINI MLELI2
  154. C
  155. SEGACT MELSOM
  156. NSOMM=MELSOM.NUM(/2)
  157. JG=NSOMM
  158. SEGINI MLECOE
  159. SEGACT MLEPOI
  160. C
  161. C**** Le MPOVAL des normales
  162. C
  163. CALL LICHT(INORM,MPNORM,TYPE,MELFAC)
  164. C SEGACT MPNORM
  165. C
  166. C**** We create the MLENTI for the MELFAC
  167. C
  168. CALL KRIPAD(MELFAC,MLEFAC)
  169. IF(IERR .NE. 0)GOTO 9999
  170. C SEGACT MELFAC
  171. C SEGINI MLEFAC
  172. C
  173. C**** Loop on the sommet to compute NVMAX
  174. C (maximum number of neighbors)
  175. C
  176. NVMAX=0
  177. DO INOEU=1,NSOMM,1
  178. C
  179. C******* The neighbors of NGS
  180. C
  181. MLENTI=MLEPOI.LECT(INOEU)
  182. SEGACT MLENTI
  183. NVOI=MLENTI.LECT(/1)
  184. NVMAX=MAX(NVMAX,NVOI)
  185. ENDDO
  186. C
  187. C MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1)
  188. C MATU = matrix of the singular right eigenvectors of MATA
  189. C (NVMAX,IDIM+1)
  190. C MATV = matrix of the singular left eigenvectors of MATA
  191. C (IDIM+1,IDIM+1)
  192. C But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1
  193. C MLRSIG = singular values of MATA (IDIM+1)
  194. C MLRB = vector (NVMAX)
  195. C MLRB.PROG(j) = 1 if j is a center point or a 'Dirichlet
  196. C boundary condition' point
  197. C MLRB.PROG(j) = (r_j - r_V) id j is a 'von Neumann
  198. C boundary condition point.
  199. C
  200. C N.B. MATA = MATU MATSIG t(MATV)
  201. C If MATA is non singular,
  202. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  203. C
  204. C MLRTRA temporary vector of invsvd.eso
  205. C NVMIN = IDIM + 1 (the most little dimension of the matrices)
  206. C
  207. N1=NVMAX
  208. N2=IDIM+1
  209. SEGINI MATA
  210. SEGINI MATU
  211. SEGINI MATV
  212. JG=NVMAX
  213. SEGINI MLRB
  214. JG=IDIM+1
  215. SEGINI MLRSIG
  216. SEGINI MLRTRA
  217. NVMIN=N2
  218. C
  219. C MTAA : [t(MATA).MATA]
  220. C MINTAA : [inve(t(MATA) MATA)]
  221. C MLRIN1,2,3 : "temporary vectors"
  222. C
  223. N1=NVMIN
  224. N2=NVMIN
  225. SEGINI MTAA
  226. SEGINI MINTAA
  227. JG=NVMIN
  228. SEGINI MLRIN1
  229. SEGINI MLRIN2
  230. SEGINI MLRIN3
  231. C
  232. C**** Loop on the sommet to compute the coefficient
  233. C
  234.  
  235. DO INOEU=1,NSOMM,1
  236. NGS=MELSOM.NUM(1,INOEU)
  237. IPCOOR=((NGS-1)*(IDIM+1))+1
  238. XS=MCOORD.XCOOR(IPCOOR)
  239. YS=MCOORD.XCOOR(IPCOOR+1)
  240. IF(IDIM .EQ. 3) ZS=MCOORD.XCOOR(IPCOOR+2)
  241. C
  242. C******* The neighbors of NGS
  243. C
  244. MLENTI=MLEPOI.LECT(INOEU)
  245. NVOI=MLENTI.LECT(/1)
  246. C
  247. C******* The matrix where we store the coefficients
  248. C
  249. JG=NVOI
  250. SEGINI MLRCOE
  251. MLECOE.LECT(INOEU)=MLRCOE
  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. NLV=MLECEN.LECT(NGV)
  259. NLV1=MLELI1.LECT(NGV)
  260. IF((NLV .NE. 0) .OR. (NLV1 .NE. 0))THEN
  261. C
  262. C********** NGV is a 'centre' point or a Dirichlet point
  263. C
  264. IPCOOR=((NGV-1)*(IDIM+1))+1
  265. XV=MCOORD.XCOOR(IPCOOR)
  266. YV=MCOORD.XCOOR(IPCOOR+1)
  267. MATA.MAT(IVOI,1)=1
  268. MATA.MAT(IVOI,2)=XV-XS
  269. MATA.MAT(IVOI,3)=YV-YS
  270. IF(IDIM.EQ.3)THEN
  271. ZV=MCOORD.XCOOR(IPCOOR+2)
  272. MATA.MAT(IVOI,4)=ZV-ZS
  273. ENDIF
  274. MLRB.PROG(IVOI)=1.0D0
  275. ELSE
  276. NLV=MLELI2.LECT(NGV)
  277. IF(NLV .EQ. 0)THEN
  278. C
  279. C**************** NO B.C. specified on this point
  280. C
  281. WRITE(IOIMP,*) 'Boundary conditions ???'
  282. CALL ERREUR(21)
  283. GOTO 9999
  284. ELSE
  285. IPCOOR=((NGV-1)*(IDIM+1))+1
  286. XV=MCOORD.XCOOR(IPCOOR)
  287. YV=MCOORD.XCOOR(IPCOOR+1)
  288. DX=((XV - XS)**2)+((YV - YS)**2)
  289. IF(IDIM .EQ. 3)THEN
  290. ZV=MCOORD.XCOOR(IPCOOR+2)
  291. DX=DX+((ZV - ZS)**2)
  292. ENDIF
  293. NLN=MLEFAC.LECT(NGV)
  294. DX=DX**0.5D0
  295. MATA.MAT(IVOI,1)=0
  296. MATA.MAT(IVOI,2)=DX*MPNORM.VPOCHA(NLN,1)
  297. MATA.MAT(IVOI,3)=DX*MPNORM.VPOCHA(NLN,2)
  298. IF(IDIM .EQ. 3)THEN
  299. MATA.MAT(IVOI,4)=DX*MPNORM.VPOCHA(NLN,3)
  300. ENDIF
  301. MLRB.PROG(IVOI)=DX
  302. ENDIF
  303. ENDIF
  304. ENDDO
  305. C
  306. C TEST
  307. C do ivoi=1,nvoi,1
  308. C write(*,*)
  309. C & 'mata.mat(',ivoi,')=',(mata.mat(ivoi,i1),i1=1,nvmin,1)
  310. C write(*,*) 'b(',ivoi,')=',mlrb.prog(ivoi)
  311. C enddo
  312. C
  313. C
  314. C********** Now we have to invert this matrix
  315. C
  316. LOGSVD=.TRUE.
  317. CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  318. & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
  319. & MLRTRA.PROG)
  320. IF(IERSVD.NE.0)THEN
  321. C
  322. C************* SVD decomposition of the matrix does not work
  323. C
  324. LOGSVD=.FALSE.
  325. ELSE
  326. C
  327. C************ We check the condition number of MATA
  328. C
  329. SMAX=0.0D0
  330. DO I1=1,NVMIN,1
  331. SMAX=MAX(SMAX,MLRSIG.PROG(I1))
  332. ENDDO
  333. SMIN=SMAX
  334. DO I1=1,NVMIN,1
  335. SMIN=MIN(SMIN,MLRSIG.PROG(I1))
  336. ENDDO
  337. IF((SMIN/SMAX) .LT. ERRTOL)THEN
  338. LOGSVD=.FALSE.
  339. ENDIF
  340. ENDIF
  341. C
  342. C TEST
  343. C write(*,*) 'LOGSVD=.FALSE.'
  344. C LOGSVD=.FALSE.
  345. C
  346. IF(LOGSVD)THEN
  347. C
  348. C********** INVSVD worked
  349. C
  350. C MATA = MATU MATSIG t(MATV)
  351. C inv(MATA) = MATV inv(MATSIG) t(MATU)
  352. C
  353. DO I4=1,NVMIN,1
  354. DO IVOI=1,NVOI,1
  355. DO I2=1,1,1
  356. C I2=1 is the only coefficient we are interested in
  357. MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
  358. & (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4)
  359. & *MLRB.PROG(IVOI)/MLRSIG.PROG(I4))
  360. ENDDO
  361. ENDDO
  362. ENDDO
  363. ELSE
  364. WRITE (IOIMP,*) 'rlenco.eso'
  365. C 22 0
  366. C Opération malvenue. Résultat douteux
  367. CALL ERREUR(22)
  368. C
  369. C********** INVSVD does not worked
  370. C For each neighbor k we have to compute the solution
  371. C of
  372. C
  373. C t(MATA) MATA x = t(MATA) * b
  374. C
  375. C where b= \sum_l e_l \delta(k,l) = e_k
  376. C
  377. C To do that, we compute
  378. C
  379. C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  380. C
  381. C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  382. C [t(MATA) MATA] X_0]
  383. C
  384. C
  385. C********** We compute [t(MATA) MATA]
  386. C We store it in the upper triangle of MTAA(NVMIN,NVMIN)
  387. C
  388. DO I1=1,NVMIN,1
  389. DO I2=I1,NVMIN,1
  390. MTAA.MAT(I1,I2)=0.0D0
  391. ENDDO
  392. ENDDO
  393.  
  394. DO I1=1,NVMIN,1
  395. DO I2=I1,NVMIN,1
  396. DO I3=1,NVOI,1
  397. MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
  398. & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
  399. ENDDO
  400. ENDDO
  401. ENDDO
  402. C
  403. C********** We compute [inve(t(MATA) MATA)]
  404. C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
  405. C
  406. CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
  407. & MLRIN2.PROG,IERR0)
  408. IF(IERR0.NE.0)THEN
  409. WRITE(IOIMP,*) 'subroutine rlenco.eso.'
  410. C 26 2
  411. C Tache impossible. Probablement données erronées
  412. CALL ERREUR(26)
  413. GOTO 9999
  414. ENDIF
  415. C
  416. C********** We complete MTAA and MINTAA
  417. C
  418. DO I1=1,NVMIN,1
  419. DO I2=I1+1,NVMIN,1
  420. MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
  421. MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
  422. ENDDO
  423. ENDDO
  424. C
  425. DO IVOI=1,NVOI,1
  426. C
  427. C************* We compute [t(MATA) . b] and we store it in MLRIN1.PROG
  428. C
  429. DO I1=1,NVMIN,1
  430. MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)*MLRB.PROG(IVOI)
  431. MLRIN2.PROG(I1)=0.0D0
  432. MLRIN3.PROG(I1)=0.0D0
  433. ENDDO
  434. C
  435. C************* X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  436. C X_0(i1) into MLRIN2.PROG(I1)
  437. C
  438. DO I2=1,NVMIN,1
  439. DO I1=1,NVMIN,1
  440. MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
  441. & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
  442. ENDDO
  443. ENDDO
  444. C
  445. C************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  446. C [t(MATA) MATA] X_0]
  447. C
  448. C [t(MATA) MATA] X_0 into MLRIN3.PROG
  449. C
  450. DO I2=1,NVMIN,1
  451. DO I1=1,NVMIN,1
  452. MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
  453. & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
  454. ENDDO
  455. ENDDO
  456. C
  457. C************* Now we have
  458. C [t(MATA) . b] in MLRIN1.PROG
  459. C X_0(i1) in MLRIN2.PROG(I1)
  460. C [t(MATA) MATA] X_0 in MLRIN3.PROG
  461. C
  462. C X_1(i1) in MLRCOE.MAT(i1,IVOI)
  463. C
  464. DO I1=1,1,1
  465. C The only interesting one is I1=1
  466. DO I2=1,NVMIN,1
  467. MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
  468. & (MINTAA.MAT(I1,I2)*
  469. & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
  470. ENDDO
  471. MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+
  472. & MLRIN2.PROG(I1)
  473. ENDDO
  474. ENDDO
  475. ENDIF
  476. CC
  477. CC TEST
  478. C write(*,*) 'ngs =', NGS
  479. C write(*,*) 'invide',LOGSVD
  480. C write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1)
  481. C write(*,*) 'coeff =',(mlrcoe.prog(ivoi),ivoi=1,nvoi,1)
  482. C dx=0.0D0
  483. C do ivoi=1,nvoi,1
  484. C dx=dx+mlrcoe.prog(ivoi)
  485. C enddo
  486. C write(*,*) 'sum=',dx
  487. MLENTI=MLEPOI.LECT(INOEU)
  488. SEGDES MLENTI
  489. MLRCOE=MLECOE.LECT(INOEU)
  490. SEGDES MLRCOE
  491. ENDDO
  492. C
  493. SEGDES MELCEN
  494. SEGSUP MLECEN
  495. C
  496. IF(MELLI1 .NE. 0) SEGDES MELLI1
  497. IF(MELLI2 .NE. 0) SEGDES MELLI2
  498. SEGSUP MLELI1
  499. SEGSUP MLELI2
  500. C
  501. SEGDES MELSOM
  502. C
  503. SEGDES MLECOE
  504. SEGDES MLEPOI
  505. C
  506. SEGSUP MATA
  507. SEGSUP MATU
  508. SEGSUP MATV
  509. SEGSUP MLRSIG
  510. SEGSUP MLRTRA
  511. SEGSUP MLRB
  512. C
  513. SEGSUP MTAA
  514. SEGSUP MINTAA
  515. SEGSUP MLRIN1
  516. SEGSUP MLRIN2
  517. SEGSUP MLRIN3
  518. C
  519. SEGDES MPNORM
  520. SEGDES MELFAC
  521. SEGSUP MLEFAC
  522. C
  523. 9999 CONTINUE
  524. RETURN
  525. END
  526.  
  527.  
  528.  
  529.  

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