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

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