Numérotation des lignes :

C RLENCO    SOURCE    CHAT      05/01/13    03:01:24     5004      SUBROUTINE RLENCO(MELSOM,MELCEN,MELLI1,MELLI2,INORM,MLEPOI,MLECOE)C************************************************************************CC PROJET            :  CASTEM 2000CC NOM               :  RLENCOCC DESCRIPTION       :  Appelle par GRADI2CC LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)CC AUTEUR            :  A. BECCANTINICC************************************************************************CC     This subroutine computes the coefficients which allow to performC     a linear exact reconstruction of each 'sommet' (vertex) valueC     Indeed, given the i-th sommet (V=MELSOM.NUM(1,i)), given the listC     of its neighbors MLEPOI.LECT(i), we have to compute the matrix ofC     coefficients a(i,j) such thatCC     VAL(MELSOM.NUM(1,i)) =  \sum_{neig1} a(i,neig1) * VAL(neig1) +C                             \sum_{neig2} a(i,neig2) * VALNOR(neig2)CC     where neig1 is a 'CENTRE' point or a 'Dirichlet boundary condition'C     point, neig2 is a 'Von Neumann boundary condition' point, VAL(j)C     is the value of the function in the j-th neighbor, VALNOR(j) isC     the value of grad(VAL).n is the j-th neighbor.CCC     To compute these coefficients, we impose that VAL is linear withC     respect to the coordinates of V and its neighbor. Then, sinceC     there are less unknowns than equations, we use a least squareC     method. If the neighbor is a 'CENTRE' point or a 'DirichletC     boundary condition' point, we impose thatCC     A + B (x_{neig} - x_V) + C (y_{neig} - y_V) + D (z_{neig} - z_V)C     =  VAL(neig)CC     If the neighbor is a 'von Neumann boundary condition' point,C     we impose thatCC     |r_{neig} - r_V| (B nx_{neig} + C ny_{neig} + D nz_{neig}) =C     VALNOR(neig) * |r_{neig} - r_V|Cc     r = (x,y,z)c     n = (nx, ny, nz) is the normal to the surfaceCC     To determine A, B, C, D we have to solveCC     MATA . tran(A, B, C, D) = tran(VAL(1), VAL(2), ... , VALNOR(j) /C                           (r_{j} - r_V)CC     with MATA(1,*) = (1,x_1-x_V,y_1-y_V,z_1-z_V)C          MATA(2,*) = (1,x_2-x_V,y_2-y_V,z_2-z_V)C          ...C          MATA(j,*) = (0,nx_{j}, ny_{j}, nz_{j})C          ...CC    To determine a(i,neig) we have to solve a linear system (for eachC    neig). If the neighbor is a 'CENTRE' point or a 'Dirichlet boundaryC    condition' point we solveCC    MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) = e_{neig}CC    with e_{neig} = tran (0,0,...,0,1,0,...) (1 in the neig position)CC    If the neighbor is a 'von Neumann boundary condition' points,C    we solveCC    MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) =C       e_{neig}  (r_{neig} - r_V)CC    Since each linear system involves MATA, we have to "pseudoinvert"C    it.CC**** Inputs:CC     MELSOM = the 'SOMMET' melemeC     MELCEN = the 'CENTRE' melemeC     MELLI1 = the support of Dirichlet boundary conditionsC     MELLI2 = the support of von Neumann boundary conditionsC     INOMR  = CHAMPOINT des normales aux facesC     MLEPOI = list of integers.C              MLEPOI.LECT(i) points to the list neighbors ofC              MELSOM.NUM(1,i)CC**** Output:CC     MLECOE =  list of integers.CC               MLECOE.LECT(i) points to the MLREEL MLRCOE whichC               contain the a(i,neig)CC**** Variables de COOPTIOCC      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERRC     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQC     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORESC     &        ,IECHO,   IIMPI,  IOSPIC     &        ,IDIMCC     &        ,MCOORDC     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVEC     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLUC     &        ,NORINC,NORVAL,NORIND,NORVADC     &        ,NUCROU, IPSAUV, IFICLE, IPREFIC      IMPLICIT INTEGER(I-N)      INTEGER INORM, NSOMM, INOEU, NGS, IPCOOR, IVOI, NVOI, NGV, NLV     &     ,NLV1,NVMAX,NVMIN,IERSVD,IERR0,NLN     &     ,I1,I2,I3,I4      REAL*8 XS,YS,ZS,XV,YV,ZV,DX,ERRTOL,SMAX,SMIN      PARAMETER(ERRTOL=1.0D-16)      CHARACTER*8 TYPE      LOGICAL LOGSVD-INC CCOPTIO-INC SMCOORD-INC SMELEME-INC SMLREEL-INC SMLENTI-INC SMCHPOI      INTEGER JG      INTEGER N1,N2      SEGMENT MATRIX      REAL*8 MAT(N1,N2)      ENDSEGMENT      POINTEUR MELSOM.MELEME, MLEPOI.MLENTI,MELLI1.MELEME,MELLI2.MELEME     &     ,MLELI1.MLENTI,MLELI2.MLENTI,MELCEN.MELEME,MLECEN.MLENTI     &     ,MPNORM.MPOVAL, MATA.MATRIX,MATU.MATRIX,MATV.MATRIX     &     ,MLRCOE.MLREEL,MLRSIG.MLREEL,MLRTRA.MLREEL     &     ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL     &     ,MLRIN3.MLREEL, MLECOE.MLENTI, MLRB.MLREEL     &     ,MELFAC.MELEME,MLEFAC.MLENTICC**** We create the MLENTI for the centersC      CALL KRIPAD(MELCEN,MLECEN)      IF(IERR .NE. 0)GOTO 9999C     En KRIPADC     SEGACT MELCENC     SEGINI MLECENCC**** We create the MLENTI for the BCC      CALL KRIPAD(MELLI1,MLELI1)      IF(IERR .NE. 0)GOTO 9999      CALL KRIPAD(MELLI2,MLELI2)      IF(IERR .NE. 0)GOTO 9999C     SEGACT MELLI1C     SEGINI MLELI1C     SEGACT MELLI2C     SEGINI MLELI2C      SEGACT MELSOM      NSOMM=MELSOM.NUM(/2)      JG=NSOMM      SEGINI MLECOE      SEGACT MLEPOICC**** Le MPOVAL des normalesC      CALL LICHT(INORM,MPNORM,TYPE,MELFAC)C     SEGACT MPNORMCC**** We create the MLENTI for the MELFACC      CALL KRIPAD(MELFAC,MLEFAC)      IF(IERR .NE. 0)GOTO 9999C     SEGACT MELFACC     SEGINI MLEFACCC**** Loop on the sommet to compute NVMAXC     (maximum number of neighbors)C      NVMAX=0      DO INOEU=1,NSOMM,1CC******* The neighbors of NGSC         MLENTI=MLEPOI.LECT(INOEU)         SEGACT MLENTI         NVOI=MLENTI.LECT(/1)         NVMAX=MAX(NVMAX,NVOI)      ENDDOCC     MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1)C     MATU = matrix of the singular right eigenvectors of MATAC            (NVMAX,IDIM+1)C     MATV = matrix of the singular left eigenvectors of MATAC            (IDIM+1,IDIM+1)C            But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1C     MLRSIG = singular values of MATA (IDIM+1)C     MLRB   = vector (NVMAX)C              MLRB.PROG(j) = 1 if j is a center point or a 'DirichletC                             boundary condition' pointC              MLRB.PROG(j) = (r_j - r_V) id j is a 'von NeumannC                              boundary condition point.CC     N.B. MATA = MATU MATSIG t(MATV)C          If MATA is non singular,C          inv(MATA) = MATV inv(MATSIG) t(MATU)CC     MLRTRA temporary vector of invsvd.esoC     NVMIN = IDIM + 1 (the most little dimension of the matrices)C      N1=NVMAX      N2=IDIM+1      SEGINI MATA      SEGINI MATU      SEGINI MATV      JG=NVMAX      SEGINI MLRB      JG=IDIM+1      SEGINI MLRSIG      SEGINI MLRTRA      NVMIN=N2CC     MTAA : [t(MATA).MATA]C     MINTAA : [inve(t(MATA) MATA)]C     MLRIN1,2,3 : "temporary vectors"C      N1=NVMIN      N2=NVMIN      SEGINI MTAA      SEGINI MINTAA      JG=NVMIN      SEGINI MLRIN1      SEGINI MLRIN2      SEGINI MLRIN3CC**** Loop on the sommet to compute the coefficientC       DO INOEU=1,NSOMM,1         NGS=MELSOM.NUM(1,INOEU)         IPCOOR=((NGS-1)*(IDIM+1))+1         XS=MCOORD.XCOOR(IPCOOR)         YS=MCOORD.XCOOR(IPCOOR+1)         IF(IDIM .EQ. 3) ZS=MCOORD.XCOOR(IPCOOR+2)CC******* The neighbors of NGSC         MLENTI=MLEPOI.LECT(INOEU)         NVOI=MLENTI.LECT(/1)CC******* The matrix where we store the coefficientsC         JG=NVOI         SEGINI MLRCOE         MLECOE.LECT(INOEU)=MLRCOECC******* Loop involving the neighborsC        We create the matrix to "pseudoinvert"C         DO IVOI=1,NVOI,1            NGV=MLENTI.LECT(IVOI)            NLV=MLECEN.LECT(NGV)            NLV1=MLELI1.LECT(NGV)            IF((NLV .NE. 0) .OR. (NLV1 .NE. 0))THENCC********** NGV is a 'centre' point or a Dirichlet pointC               IPCOOR=((NGV-1)*(IDIM+1))+1               XV=MCOORD.XCOOR(IPCOOR)               YV=MCOORD.XCOOR(IPCOOR+1)               MATA.MAT(IVOI,1)=1               MATA.MAT(IVOI,2)=XV-XS               MATA.MAT(IVOI,3)=YV-YS               IF(IDIM.EQ.3)THEN                  ZV=MCOORD.XCOOR(IPCOOR+2)                  MATA.MAT(IVOI,4)=ZV-ZS               ENDIF               MLRB.PROG(IVOI)=1.0D0            ELSE               NLV=MLELI2.LECT(NGV)               IF(NLV .EQ. 0)THENCC**************** NO B.C. specified on this pointC                  WRITE(IOIMP,*) 'Boundary conditions ???'                  CALL ERREUR(21)                  GOTO 9999               ELSE                  IPCOOR=((NGV-1)*(IDIM+1))+1                  XV=MCOORD.XCOOR(IPCOOR)                  YV=MCOORD.XCOOR(IPCOOR+1)                  DX=((XV - XS)**2)+((YV - YS)**2)                  IF(IDIM .EQ. 3)THEN                     ZV=MCOORD.XCOOR(IPCOOR+2)                     DX=DX+((ZV - ZS)**2)                  ENDIF                  NLN=MLEFAC.LECT(NGV)                  DX=DX**0.5D0                  MATA.MAT(IVOI,1)=0                  MATA.MAT(IVOI,2)=DX*MPNORM.VPOCHA(NLN,1)                  MATA.MAT(IVOI,3)=DX*MPNORM.VPOCHA(NLN,2)                  IF(IDIM .EQ. 3)THEN                     MATA.MAT(IVOI,4)=DX*MPNORM.VPOCHA(NLN,3)                  ENDIF                  MLRB.PROG(IVOI)=DX               ENDIF            ENDIF         ENDDOCC        TESTC         do ivoi=1,nvoi,1C            write(*,*)C     &           'mata.mat(',ivoi,')=',(mata.mat(ivoi,i1),i1=1,nvmin,1)C            write(*,*) 'b(',ivoi,')=',mlrb.prog(ivoi)C         enddoCCC********** Now we have to invert this matrixC         LOGSVD=.TRUE.         CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,     &        MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,     &        MLRTRA.PROG)         IF(IERSVD.NE.0)THENCC************* SVD decomposition of the matrix does not workC            LOGSVD=.FALSE.         ELSECC************ We check the condition number of MATAC            SMAX=0.0D0            DO I1=1,NVMIN,1               SMAX=MAX(SMAX,MLRSIG.PROG(I1))            ENDDO            SMIN=SMAX            DO I1=1,NVMIN,1               SMIN=MIN(SMIN,MLRSIG.PROG(I1))            ENDDO            IF((SMIN/SMAX) .LT. ERRTOL)THEN               LOGSVD=.FALSE.            ENDIF         ENDIFCC        TESTC         write(*,*) 'LOGSVD=.FALSE.'C         LOGSVD=.FALSE.C         IF(LOGSVD)THENCC********** INVSVD workedCC           MATA = MATU MATSIG t(MATV)C           inv(MATA) = MATV inv(MATSIG) t(MATU)C            DO I4=1,NVMIN,1               DO IVOI=1,NVOI,1                  DO I2=1,1,1C                    I2=1 is the only coefficient we are interested in                     MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+     &                    (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4)     &                    *MLRB.PROG(IVOI)/MLRSIG.PROG(I4))                  ENDDO               ENDDO            ENDDO         ELSE            WRITE (IOIMP,*) 'rlenco.eso'C           22 0C           Opération malvenue. Résultat douteux            CALL ERREUR(22)CC********** INVSVD does not workedC           For each neighbor k we have to compute the solutionC           ofCC           t(MATA) MATA x = t(MATA) * bCC           where b= \sum_l e_l \delta(k,l) = e_kCC           To do that, we computeCC           X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]CC           X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -C                  [t(MATA) MATA] X_0]CCC********** We compute [t(MATA) MATA]C           We store it in the upper triangle of MTAA(NVMIN,NVMIN)C            DO I1=1,NVMIN,1               DO I2=I1,NVMIN,1                  MTAA.MAT(I1,I2)=0.0D0               ENDDO            ENDDO             DO I1=1,NVMIN,1               DO I2=I1,NVMIN,1                  DO I3=1,NVOI,1                     MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+     &                    (MATA.MAT(I3,I1)*MATA.MAT(I3,I2))                  ENDDO               ENDDO            ENDDOCC********** We compute [inve(t(MATA) MATA)]C           CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN)C            CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT,     &           MLRIN2.PROG,IERR0)            IF(IERR0.NE.0)THEN               WRITE(IOIMP,*) 'subroutine rlenco.eso.'C              26 2C              Tache impossible. Probablement données erronées               CALL ERREUR(26)               GOTO 9999            ENDIFCC********** We complete MTAA and MINTAAC            DO I1=1,NVMIN,1               DO I2=I1+1,NVMIN,1                  MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2)                  MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2)               ENDDO            ENDDOC            DO IVOI=1,NVOI,1CC************* We compute [t(MATA) . b]  and we store it in MLRIN1.PROGC               DO I1=1,NVMIN,1                  MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1)*MLRB.PROG(IVOI)                  MLRIN2.PROG(I1)=0.0D0                  MLRIN3.PROG(I1)=0.0D0               ENDDOCC************* X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]C              X_0(i1) into MLRIN2.PROG(I1)C               DO I2=1,NVMIN,1                  DO I1=1,NVMIN,1                     MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+     &                    (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2))                  ENDDO               ENDDOCC************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -C                  [t(MATA) MATA] X_0]CC             [t(MATA) MATA] X_0 into MLRIN3.PROGC               DO I2=1,NVMIN,1                  DO I1=1,NVMIN,1                     MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+     &                    (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))                  ENDDO               ENDDOCC************* Now we haveC              [t(MATA) . b]  in MLRIN1.PROGC              X_0(i1)        in MLRIN2.PROG(I1)C              [t(MATA) MATA] X_0 in MLRIN3.PROGCC              X_1(i1) in MLRCOE.MAT(i1,IVOI)C               DO I1=1,1,1C                 The only interesting one is I1=1                  DO I2=1,NVMIN,1                     MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+     &                    (MINTAA.MAT(I1,I2)*     &                    (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))                  ENDDO                  MLRCOE.PROG(IVOI)=MLRCOE.PROG(IVOI)+     &                 MLRIN2.PROG(I1)               ENDDO            ENDDO         ENDIFCCCC        TESTC         write(*,*) 'ngs =', NGSC         write(*,*) 'invide',LOGSVDC         write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1)C         write(*,*) 'coeff =',(mlrcoe.prog(ivoi),ivoi=1,nvoi,1)C         dx=0.0D0C         do ivoi=1,nvoi,1C            dx=dx+mlrcoe.prog(ivoi)C         enddoC         write(*,*) 'sum=',dx         MLENTI=MLEPOI.LECT(INOEU)         SEGDES MLENTI         MLRCOE=MLECOE.LECT(INOEU)         SEGDES MLRCOE      ENDDOC      SEGDES MELCEN      SEGSUP MLECENC      IF(MELLI1 .NE. 0) SEGDES MELLI1      IF(MELLI2 .NE. 0) SEGDES MELLI2      SEGSUP MLELI1      SEGSUP MLELI2C      SEGDES MELSOMC      SEGDES MLECOE      SEGDES MLEPOIC      SEGSUP MATA      SEGSUP MATU      SEGSUP MATV      SEGSUP MLRSIG      SEGSUP MLRTRA      SEGSUP MLRBC      SEGSUP MTAA      SEGSUP MINTAA      SEGSUP MLRIN1      SEGSUP MLRIN2      SEGSUP MLRIN3C      SEGDES MPNORM      SEGDES MELFAC      SEGSUP MLEFACC 9999 CONTINUE      RETURN      END

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