Télécharger la plaquette Cast3M


fr en

Télécharger rlexca.eso

Retour à la liste

Numérotation des lignes : 

rlexca

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

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