Télécharger la plaquette Cast3M


fr en

Télécharger rleca1.eso

Retour à la liste

Numérotation des lignes : 

rleca1

  1. C RLECA1    SOURCE    CHAT      05/01/13    03:01:02     5004
  2.       SUBROUTINE RLECA1(MELFL,MELFP,MLEFSC,MACOE1)
  3. C************************************************************************
  4. C
  5. C PROJET            :  CASTEM 2000
  6. C
  7. C NOM               :  RLECA1
  8. C
  9. C DESCRIPTION       :  Appelle par GRADIA
  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 MELFL  : MELEME FACEL
  20. C
  21. C MELFP  : MELEME FACEP
  22. C
  23. C Outputs:
  24. C
  25. C MLEFSC : list of FACES points and their neighbors (CENTRE and SOMMET
  26. C      points)
  27. C
  28. C MACOE1 : coeff. for linear exact reconstruction of MLEFSC
  29. C
  30. C
  31. C**** Variables de COOPTIO
  32. C
  33. C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  34. C     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
  35. C     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
  36. C     &        ,IECHO,   IIMPI,  IOSPI
  37. C     &        ,IDIM
  38. CC     &        ,MCOORD
  39. C     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  40. C     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  41. C     &        ,NORINC,NORVAL,NORIND,NORVAD
  42. C     &        ,NUCROU, IPSAUV, IFICLE, IPREFI
  43. C
  44.       IMPLICIT INTEGER(I-N)
  45. -INC PPARAM
  46. -INC CCOPTIO
  47. -INC SMCOORD
  48. -INC SMELEME
  49.       POINTEUR MELFL.MELEME,MELFP.MELEME,MELFP1.MELEME
  50.       INTEGER JG
  51. -INC SMLENTI
  52.       POINTEUR MLEFP.MLENTI
  53. -INC SMLREEL
  54.       POINTEUR MLRSIG.MLREEL, MLRTRA.MLREEL
  55.      &     ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL
  56.      &     ,MLRIN3.MLREEL
  57. C
  58.       INTEGER NBL, NBTPOI
  59.       SEGMENT MLELEM
  60.       INTEGER INDEX(NBL+1)
  61.       INTEGER LESPOI(NBTPOI)
  62.       ENDSEGMENT
  63.       POINTEUR MLEFSC.MLELEM
  64. C
  65.       INTEGER N1,N2
  66.       SEGMENT MATRIX
  67.       REAL*8 MAT(N1,N2)
  68.       ENDSEGMENT
  69.       POINTEUR MATA.MATRIX,MATU.MATRIX,MATV.MATRIX
  70.      &     ,MACOE1.MATRIX
  71.       LOGICAL LOGSVD
  72. C
  73.       INTEGER NBSOUS, NVMAX, NELT, NNOE, ISOUS, NVMIN, NVOI, IELEM
  74.      &     ,  IPOS, IELEM1,NGF,NGF1, NGN, IPCOOR, I1, I2, I3, IPVOI
  75.      &     , IERSVD, NGCD, NGCG, IVOI,IERR0,I4
  76.  
  77.       REAL*8 ERRTOL,SMAX,SMIN, XF, YF, ZF
  78.       PARAMETER (ERRTOL=1.0D0-15)
  79. C
  80.       NBL = 0
  81.       NBTPOI = 0
  82.       NVMAX = 0
  83.       SEGACT MELFP
  84.       NBSOUS=MELFP.LISOUS(/1)
  85. C     NBSOUS=0 fais un peux chier!
  86.       JG=MAX(NBSOUS,1)
  87.       SEGINI MLEFP
  88.       IF(NBSOUS .EQ. 0)THEN
  89.          MLEFP.LECT(1)=MELFP
  90.       ELSE
  91.          DO ISOUS=1,NBSOUS,1
  92.             MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS)
  93.          ENDDO
  94.       ENDIF
  95.       SEGDES MELFP
  96.       NBSOUS=JG
  97. C
  98.       DO ISOUS=1,NBSOUS,1
  99.          MELFP1=MLEFP.LECT(ISOUS)
  100.          SEGACT MELFP1
  101.          NNOE=MELFP1.NUM(/1)-1
  102.          NELT=MELFP1.NUM(/2)
  103.          NBTPOI=NBTPOI+((1+2+NNOE)*NELT)
  104. C
  105. C        1 = centre de face
  106. C        2 = centre d'elts "gauche" et "droite" (si differents)
  107. C        NNOE = sommet voisins
  108. C
  109.          NBL=NBL+NELT
  110.          NVMAX=MAX(NVMAX,(2+NNOE))
  111.       ENDDO
  112. C
  113. C     N.B.: NBTPOI surestimé
  114. C
  115. C**** MLEFSC
  116. C     Structure Face-(Sommets-Centres)
  117. C
  118.       SEGINI MLEFSC
  119. C
  120. C**** La matrice de coeff.
  121. C
  122. C     MACOE1(*,I) = coefficients concernants MLEFSC.LESPOI(I)
  123. C
  124.       NVMIN=IDIM+1
  125. C
  126. C**** NVMIN = nombre minimum de voisins au dessus duquel MATA est
  127. C             singulier
  128. C
  129.       N1=NVMIN
  130.       N2=NBTPOI
  131.       SEGINI MACOE1
  132. C
  133. C**** MATA = matrice à inverser (NVMAX,IDIM+1)
  134. C     MATU = matice des vectors propres singuliers droite de MATA
  135. C            (NVMAX,IDIM+1)
  136. C     MATV = matrice des vectors propres singuliers gauche de MATA
  137. C            (IDIM+1,IDIM+1)
  138. C            Mais en INVSVD a dimension NVMAX,IDIM+1
  139. C     MLRSIG =liste des valeurs singuliers de MATA
  140. C
  141. C     N.B. MATA = MATU MATSIG t(MATV)
  142. C          Si MATA non singulier
  143. C          inv(MATA) = MATV inv(MATSIG) t(MATU)
  144. C
  145.       N1=NVMAX
  146.       N2=NVMIN
  147.       SEGINI MATA
  148.       SEGINI MATU
  149.       SEGINI MATV
  150.       JG=NVMIN
  151.       SEGINI MLRSIG
  152.       SEGINI MLRTRA
  153. C     MLRTRA segment de travail de la subroutine invsvd
  154. C
  155. C     MTAA : [t(MATA).MATA]
  156. C     MINTAA : [inve(t(MATA) MATA)]
  157. C     MLRIN1,2,3 : "temporary vectors"
  158. C
  159.       N1=NVMIN
  160.       N2=NVMIN
  161.       SEGINI MTAA
  162.       SEGINI MINTAA
  163.       JG=NVMIN
  164.       SEGINI MLRIN1
  165.       SEGINI MLRIN2
  166.       SEGINI MLRIN3
  167. C
  168.       SEGACT MELFL
  169.       MLEFSC.INDEX(1)=1
  170.       IELEM=0
  171.       DO ISOUS=1,NBSOUS,1
  172.          MELFP1=MLEFP.LECT(ISOUS)
  173.          NNOE=MELFP1.NUM(/1)-1
  174.          NELT=MELFP1.NUM(/2)
  175.          DO IELEM1=1,NELT,1
  176.             IELEM=IELEM+1
  177.             IPOS=MLEFSC.INDEX(IELEM)
  178.             NGF=MELFP1.NUM(NNOE+1,IELEM1)
  179.             NGF1=MELFL.NUM(2,IELEM)
  180.             IF(NGF .NE. NGF1)THEN
  181.                WRITE(IOIMP,*) 'subroutine rleca1.eso'
  182.                CALL ERREUR(5)
  183.                GOTO 9999
  184.             ENDIF
  185.             IPCOOR=(IDIM+1)*(NGF-1)+1
  186.             XF=MCOORD.XCOOR(IPCOOR)
  187.             YF=MCOORD.XCOOR(IPCOOR+1)
  188.             IF(IDIM .EQ. 3) ZF=MCOORD.XCOOR(IPCOOR+2)
  189.             MLEFSC.LESPOI(IPOS)=NGF
  190.             DO I1 = 1, NVMIN, 1
  191.                MACOE1.MAT(I1,IPOS)=0.0D0
  192.             ENDDO
  193. C
  194. C********** Le centre "gauche"
  195. C
  196.             NGCG=MELFL.NUM(1,IELEM)
  197.             NGN=NGCG
  198. C           NVOI=nombre de voisins
  199.             NVOI=1
  200.             MLEFSC.LESPOI(IPOS+1)=NGN
  201.             MATA.MAT(1,1)=1.0D0
  202.             IPCOOR=(IDIM+1)*(NGN-1)+1
  203.             MATA.MAT(1,2)=MCOORD.XCOOR(IPCOOR)-XF
  204.             MATA.MAT(1,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  205.             IF(IDIM.EQ.3) MATA.MAT(1,4)=MCOORD.XCOOR(IPCOOR+2)-ZF
  206. C
  207. C********** Le centre "droite"
  208. C
  209.             NGCD=MELFL.NUM(3,IELEM)
  210.             IF(NGCD .NE. NGCG)THEN
  211.                NVOI=NVOI+1
  212.                NGN=NGCD
  213.                MLEFSC.LESPOI(IPOS+2)=NGN
  214.                MATA.MAT(2,1)=1.0D0
  215.                IPCOOR=(IDIM+1)*(NGN-1)+1
  216.                MATA.MAT(2,2)=MCOORD.XCOOR(IPCOOR)-XF
  217.                MATA.MAT(2,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  218.                IF(IDIM.EQ.3) MATA.MAT(2,4)=MCOORD.XCOOR(IPCOOR+2)-ZF
  219.             ENDIF
  220. C
  221. C********** Les sommets
  222. C
  223.             DO I1 = 1, NNOE, 1
  224.                NGN=MELFP1.NUM(I1,IELEM1)
  225.                MLEFSC.LESPOI(IPOS+NVOI+I1)=NGN
  226.                MATA.MAT(NVOI+I1,1)=1.0D0
  227.                IPCOOR=(IDIM+1)*(NGN-1)+1
  228.                MATA.MAT(NVOI+I1,2)=MCOORD.XCOOR(IPCOOR)-XF
  229.                MATA.MAT(NVOI+I1,3)=MCOORD.XCOOR(IPCOOR+1)-YF
  230.                IF(IDIM.EQ.3) MATA.MAT(NVOI+I1,4)=MCOORD.XCOOR(IPCOOR+2)
  231.      &              -ZF
  232.             ENDDO
  233.             NVOI=NVOI+NNOE
  234.             MLEFSC.INDEX(IELEM+1)=IPOS+NVOI+1
  235. C
  236.             LOGSVD=.TRUE.
  237.             CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
  238.      &           MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
  239.      &           MLRTRA.PROG)
  240.             IF(IERSVD.NE.0)THEN
  241.                LOGSVD=.FALSE.
  242.             ELSE
  243.                SMAX=0.0D0
  244.                DO I2=1,NVMIN,1
  245.                   SMAX=MAX(SMAX,MLRSIG.PROG(I2))
  246.                ENDDO
  247.                SMIN=SMAX
  248.                DO I2=1,NVMIN,1
  249.                   SMIN=MIN(SMIN,MLRSIG.PROG(I2))
  250.                ENDDO
  251.             ENDIF
  252.             IF((SMIN/SMAX) .LT. ERRTOL)THEN
  253.                LOGSVD=.FALSE.
  254.             ENDIF
  255. C
  256. C********** We have two different cases:
  257. C           LOGSVD = .TRUE. -> we have the coeff. we are looking for
  258. C           LOGSVD = .FALSE. -> we have not
  259. C
  260. C            LOGSVD=.FALSE.
  261.             IF(LOGSVD)THEN
  262.                DO I4=1,NVMIN,1
  263.                   DO I3=1,NVOI,1
  264.                      IPVOI=IPOS+I3
  265.                      DO I2=1,NVMIN,1
  266.                         MACOE1.MAT(I2,IPVOI)=MACOE1.MAT(I2,IPVOI)+
  267.      &                       (MATV.MAT(I2,I4)*MATU.MAT(I3,I4)/
  268.      &                       MLRSIG.PROG(I4))
  269.                      ENDDO
  270.                   ENDDO
  271.                ENDDO
  272.             ELSE
  273.                WRITE (IOIMP,*) 'rleca1.eso'
  274. C              22 0
  275. C              Opération malvenue. Résultat douteux
  276.                CALL ERREUR(22)
  277. C
  278. C************* INVSVD does not worked
  279. C              For each neighbor k we have to compute the solution
  280. C              of
  281. C
  282. C              t(MATA) MATA x = t(MATA) * b
  283. C
  284. C              where b= \sum_l e_l \delta(k,l) = e_k
  285. C
  286. C              To do that, we compute
  287. C
  288. C              X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  289. C
  290. C              X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  291. C                   [t(MATA) MATA] X_0]
  292. C
  293. C
  294. C************* We compute [t(MATA) MATA]
  295. C              We store it in the upper triangle of MTAA(NVMIN,NVMIN)
  296. C
  297.                DO I1=1,NVMIN,1
  298.                   DO I2=I1,NVMIN,1
  299.                      MTAA.MAT(I1,I2)=0.0D0
  300.                   ENDDO
  301.                ENDDO
  302. C
  303.                DO I1=1,NVMIN,1
  304.                   DO I2=I1,NVMIN,1
  305.                      DO I3=1,NVOI,1
  306.                         MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+
  307.      &                       (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))
  308.                      ENDDO
  309.                   ENDDO
  310.                ENDDO
  311. C
  312. C************* We compute [inve(t(MATA) MATA)]
  313. C              CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)
  314. C
  315.                CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,
  316.      &              MLRIN2.PROG,IERR0)
  317.                IF(IERR0.NE.0)THEN
  318.                   WRITE(IOIMP,*) 'subroutine rleca1.eso.'
  319. C                 26 2
  320. C                 Tache impossible. Probablement données erronées
  321.                   CALL ERREUR(26)
  322.                   GOTO 9999
  323.                ENDIF
  324. C
  325. C************* We complete MTAA and MINTAA
  326. C
  327.                DO I1=1,NVMIN,1
  328.                   DO I2=I1+1,NVMIN,1
  329.                      MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)
  330.                      MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)
  331.                   ENDDO
  332.                ENDDO
  333. C
  334.                DO IVOI=1,NVOI,1
  335. C
  336. C**************** We compute [t(MATA) . b]  and we store it in MLRIN1.PROG
  337. C
  338.                   DO I1=1,NVMIN,1
  339.                      MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)
  340.                      MLRIN2.PROG(I1)=0.0D0
  341.                      MLRIN3.PROG(I1)=0.0D0
  342.                   ENDDO
  343. C
  344. C**************** X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
  345. C                 X_0(i1) into MLRIN2.PROG(I1)
  346. C
  347.                   DO I2=1,NVMIN,1
  348.                      DO I1=1,NVMIN,1
  349.                         MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+
  350.      &                       (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))
  351.                      ENDDO
  352.                   ENDDO
  353. C
  354. C**************** X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
  355. C                       [t(MATA) MATA] X_0]
  356. C
  357. C                 [t(MATA) MATA] X_0 into MLRIN3.PROG
  358. C
  359.                   DO I2=1,NVMIN,1
  360.                      DO I1=1,NVMIN,1
  361.                         MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
  362.      &                       (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
  363.                      ENDDO
  364.                   ENDDO
  365. C
  366. C**************** Now we have
  367. C                 [t(MATA) . b]  in MLRIN1.PROG
  368. C                 X_0(i1)        in MLRIN2.PROG(I1)
  369. C                 [t(MATA) MATA] X_0 in MLRIN3.PROG
  370. C
  371. C                 X_1(i1) in MLRCOE.MAT(i1,IVOI)
  372. C
  373.                   IPVOI=IPOS+IVOI
  374.                   DO I1=1,NVMIN,1
  375.                      DO I2=1,NVMIN,1
  376.                         MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+
  377.      &                       (MINTAA.MAT(I1,I2)*
  378.      &                       (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))
  379.                      ENDDO
  380.                      MACOE1.MAT(I1,IPVOI)=MACOE1.MAT(I1,IPVOI)+
  381.      &                    MLRIN2.PROG(I1)
  382.                   ENDDO
  383.                ENDDO
  384.             ENDIF
  385. C
  386. C********** If we are here, we have the coeff. in the 'FACE' frame.
  387. C
  388.          ENDDO
  389.          SEGDES MELFP1
  390.       ENDDO
  391. C
  392.       NBTPOI=MLEFSC.INDEX(NBL+1)-1
  393.       SEGADJ MLEFSC
  394.       N1=NVMIN
  395.       N2=NBTPOI
  396.       SEGADJ MACOE1
  397. C
  398.       SEGSUP MLEFP
  399.       SEGDES MLEFSC
  400.       SEGDES MACOE1
  401.       SEGSUP MATA
  402.       SEGSUP MATU
  403.       SEGSUP MATV
  404.       SEGSUP MLRSIG
  405.       SEGSUP MLRTRA
  406. C
  407.       SEGDES MELFL
  408.       SEGSUP MTAA
  409.       SEGSUP MINTAA
  410.       SEGSUP MLRIN1
  411.       SEGSUP MLRIN2
  412.       SEGSUP MLRIN3
  413. C
  414.  9999 RETURN
  415.       END
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  

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