Télécharger la plaquette Cast3M


fr en

Télécharger rlencf.eso

Retour à la liste

Numérotation des lignes : 

rlencf

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

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