Numérotation des lignes :

C RLEXCE    SOURCE    CHAT      05/01/13    03:02:22     5004      SUBROUTINE RLEXCE(MELEME,MELCEN,MELFAC,INORM,ICHCL,MLECOE)C************************************************************************CC PROJET            :  CASTEM 2000CC NOM               :  RLEXCECC DESCRIPTION       :  Appelle par GRADGECC 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 'CENTRE' value.C     Indeed, given the i-th centre (NC=MELEME.NUM(NBNN,i)), given the listC     of its neighbors MELEME.NUM(j,i),j=1,NBNN-1,C     we have to compute the matrix of  coefficients a(i,j) such thatCC     VAL(MELEME.NUM(NBNN,i)) = \sum_{neig1} a(i,neig1) * VAL(neig1)CC     where neig1 is a 'CENTRE' point or a 'wall boundary condition'C     point.CC     To compute these coefficients, we impose that VAL is linear withC     respect to the coordinates of NC and its neighbor. Then, sinceC     there are less unknowns than equations, we use a least squareC     method. We impose thatCC     A + B (x_{neig} - x_NC) + C (y_{neig} - y_NC) + D (z_{neig} - z_NC)C     =  VAL(neig)CC     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_NC)CC     with MATA(1,*) = (1,x_1-x_NC,y_1-y_NC,z_1-z_NC)C          MATA(2,*) = (1,x_2-x_NC,y_2-y_NC,z_2-z_NC)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 'wall 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    N.B. Concerning the neighbor, two different cases are taken intoC         account.C    1) Reflecting boundary conditions:C       The neighbor does not belong to the geometrical support ofC       ICHCL -> It is a wall pointC       Then the reconstruction is performed on a virtual point,C       the symmetric point of the centre with respect to the wallC    2) Non_reflecting boundary conditions:C       The neighbor belongs to the geometrical support of ICHCLC       -> The reconstruction is performed on such FACE pointCC**** Inputs:CC     MELEME = the MELEME which contains the stencil of the gradientC              moleculeC     MELCEN = the 'CENTRE' melemeC     MELFAC = the 'FACE' melemeC     INORM  = the CHAMPOINT containing the face normalsC     ICHCL  = the CHAMPOINT containing the non-reflecting boundaryC              conditionsCC**** Output:CC     MLECOE =  list of integers.CC               MLECOE.LECT(i) points to the MLREEL MATCOE whichC               contain the a(i,neig)CC***********************************************************************CC Modified the 25-th of february to take into account the 'MODE' 'AXIS'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, IPREFI, IREFOR, ISAFORCC      IMPLICIT INTEGER(I-N)      INTEGER INORM, ICHCL, IGEOM, NSOU, JG, NTCEN, NVMAX, NVMIN     &     ,NBNN, NBELEM, NGV     &     , NGC, IPCOOR, NLV, NLV1, NLF, IERSVD, IERR0     &     ,ISOUS , ICEN, IELEM, IVOI, I1,  I2, I3, I4, NAXI      REAL*8 XC ,YC, ZC, XV, YV, ZV, DX, DY, DZ, ERRTOL, SMAX, SMIN     &     ,RNX, RNY, RNZ, ALPHA      PARAMETER(ERRTOL=1.0D-16)      CHARACTER*8 TYPE      LOGICAL LOGSVD, LOGAXI-INC CCOPTIO-INC SMCOORD-INC SMELEME-INC SMLREEL-INC SMLENTI-INC SMCHPOIC      INTEGER N1,N2      SEGMENT MATRIX      REAL*8 MAT(N1,N2)      ENDSEGMENTC      POINTEUR MELFAC.MELEME, MLEFAC.MLENTI, MLEBC.MLENTI, MPNORM.MPOVAL     &     ,MLESOU.MLENTI, MELCEN.MELEME ,MLECEN.MLENTI     &     ,MLECOE.MLENTI, MATCOE.MATRIX     &     ,MATA.MATRIX,MATU.MATRIX,MATV.MATRIX,MLRB.MLREEL     &     ,MLRSIG.MLREEL,MLRTRA.MLREEL,MTAA.MATRIX, MINTAA.MATRIX     &     ,MLRIN1.MLREEL,MLRIN2.MLREEL     &     ,MLRIN3.MLREELCC**** Axis-symmetrical caseC      IF((IFOMOD .EQ. 0) .AND. (IDIM .EQ. 2))THEN         LOGAXI=.TRUE.         NAXI=2      ELSE         NAXI=0         LOGAXI=.FALSE.      ENDIFCC**** We create the MLENTI for the centersC      CALL KRIPAD(MELCEN,MLECEN)      IF(IERR .NE. 0)GOTO 9999C     En KRIPADC     SEGINI MLECENCC**** We create the MLENTI for the facesC      CALL KRIPAD(MELFAC,MLEFAC)      IF(IERR .NE. 0)GOTO 9999C     En KRIPADC     SEGINI MLEFACCC**** We create the MLENTI for the BCC      IF(ICHCL.NE.0)THEN         CALL LICHT(ICHCL,MPOVAL,TYPE,IGEOM)C        In LICHTC        SEGACT MPOVAL*MOD         SEGDES MPOVAL      ELSE         IGEOM=0      ENDIF      CALL KRIPAD(IGEOM,MLEBC)      IF(IERR .NE. 0)GOTO 9999C     SEGINI MLEBCCC**** Le MPOVAL des normalesC      CALL LICHT(INORM,MPNORM,TYPE,IGEOM)C     SEGACT MPNORMCC**** We recover the elemenary mesh of MELEMEC     We compute the number of maximum number of neighborsC     We compute the number of centersC      SEGACT MELEME      NSOU=MAX(MELEME.LISOUS(/1),1)      JG=NSOU      SEGINI MLESOU      IF (NSOU.EQ.1)THEN         MLESOU.LECT(1)=MELEME         NBNN=MELEME.NUM(/1)         NTCEN=MELEME.NUM(/2)         NVMAX=NBNN      ELSE         NVMAX=0         NTCEN=0         DO ISOUS=1,NSOU,1            IPT1=MELEME.LISOUS(ISOUS)            MLESOU.LECT(ISOUS)=IPT1            SEGACT IPT1            NBNN=IPT1.NUM(/1)            NVMAX=MAX(NVMAX,NBNN)            NTCEN=NTCEN+IPT1.NUM(/2)         ENDDO      ENDIF      NVMAX=NVMAX+NAXICC**** The outputC      JG=NTCEN      SEGINI MLECOECC     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) = 1CC     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      ICEN=0      DO ISOUS=1,NSOU,1         IPT1=MLESOU.LECT(ISOUS)         NBNN=IPT1.NUM(/1)         NBELEM=IPT1.NUM(/2)         DO IELEM=1,NBELEM,1C            NGC=IPT1.NUM(NBNN,IELEM)            IPCOOR=((NGC-1)*(IDIM+1))+1            XC=MCOORD.XCOOR(IPCOOR)            YC=MCOORD.XCOOR(IPCOOR+1)            IF(IDIM.EQ.3) ZC=MCOORD.XCOOR(IPCOOR+2)CC********** We create the matrix of coefficients.C            N2=NBNN+NAXI            N1=NVMIN            ICEN=ICEN+1            SEGINI MATCOE            MLECOE.LECT(ICEN)=MATCOECC********** Loop involving the neighbors (and the centre itself)C           We create the matrix to "pseudoinvert"C            DO IVOI=1,NBNN,1               NGV=IPT1.NUM(IVOI,IELEM)               NLV=MLECEN.LECT(NGV)               NLV1=MLEBC.LECT(NGV)               IF((NLV .NE. 0) .OR. (NLV1 .NE. 0))THENCC**************** NGV is a 'centre' point or a ICHCL 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                  MATA.MAT(IVOI,3)=YV                  IF(IDIM.EQ.3)THEN                     ZV=MCOORD.XCOOR(IPCOOR+2)                     MATA.MAT(IVOI,4)=ZV                  ENDIF               ELSECC**************** Reflecting BC reconstructionC                  NLF=MLEFAC.LECT(NGV)                  RNX=MPNORM.VPOCHA(NLF,1)                  RNY=MPNORM.VPOCHA(NLF,2)                  IPCOOR=((NGV-1)*(IDIM+1))+1                  XV=MCOORD.XCOOR(IPCOOR)                  YV=MCOORD.XCOOR(IPCOOR+1)                  ALPHA=(XV-XC)*RNX+(YV-YC)*RNY                  DX=ALPHA*2.0D0*RNX                  DY=ALPHA*2.0D0*RNY                  MATA.MAT(IVOI,1)=1                  MATA.MAT(IVOI,2)=XC+DX                  MATA.MAT(IVOI,3)=YC+DY                  IF(IDIM.EQ.3)THEN                     RNZ=MPNORM.VPOCHA(NLF,3)                     ZV=MCOORD.XCOOR(IPCOOR+2)                     ALPHA=ALPHA+(ZV-ZC)*RNZ                     DZ=ALPHA*2.0D0*RNZ                     MATA.MAT(IVOI,4)=ZC+DZ                  ENDIF               ENDIF               MLRB.PROG(IVOI)=1.0D0            ENDDOC            DO IVOI=NBNN+1,NBNN+NAXI,1               MATA.MAT(IVOI,1)=1               MATA.MAT(IVOI,2)=XC               MATA.MAT(IVOI,3)=YC               MLRB.PROG(IVOI)=1.0D0            ENDDOCCCCC           TESTCCC            do ivoi=1,nbnn,1C               write(*,*) 'ngv =', ipt1.num(ivoi,ielem)C               write(*,*)C     &              'mata.mat(',ivoi,')=',(mata.mat(ivoi,i1),i1=1,nvminC     \$              ,1)C               write(*,*) 'b(',ivoi,')=',mlrb.prog(ivoi)C            enddoCCC********** Now we have to invert this matrixC            LOGSVD=.TRUE.            CALL INVSVD( NVMAX, NBNN+NAXI, 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,NBNN+NAXI,1                     DO I2=1,IDIM+1,1C                       I2=1 is the only coefficient we are not interestedC                       in. But we computed it to verify thatC                       sum_ivoi MATCOE.MAT(ivoi,1) = 1                        MATCOE.MAT(I2,IVOI)=MATCOE.MAT(I2,IVOI)+     &                       (MATV.MAT(I2,I4)*MATU.MAT(IVOI,I4)     &                       /MLRSIG.PROG(I4))                     ENDDO                  ENDDO               ENDDO            ELSE               WRITE (IOIMP,*) 'rlexce.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               ENDDOC               DO I1=1,NVMIN,1                  DO I2=I1,NVMIN,1                     DO I3=1,NBNN+NAXI,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 rlexce.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,NBNN+NAXI,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,IDIM+1,1C                    The only unuseful one is I1=1                     DO I2=1,NVMIN,1                        MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+     &                       (MINTAA.MAT(I1,I2)*     &                       (MLRIN1.PROG(I2)-MLRIN3.PROG(I2)))                     ENDDO                     MATCOE.MAT(I1,IVOI)=MATCOE.MAT(I1,IVOI)+     &                    MLRIN2.PROG(I1)                  ENDDO               ENDDO            ENDIFCC           LOGAXI -> We eliminate the false neighborsC            IF(LOGAXI)THEN               DO I1=1,NVMIN,1                  MATCOE.MAT(I1,NBNN)=MATCOE.MAT(I1,NBNN)+     &                 MATCOE.MAT(I1,NBNN+1)+MATCOE.MAT(I1,NBNN+2)               ENDDO               N2=NBNN               N1=NVMIN               SEGADJ MATCOE            ENDIFCCCC           TESTC            write(*,*) 'ngc =', NGCC            write(*,*) 'invide',LOGSVDC            write(*,*) 'nvois =',(ipt1.num(ivoi,ielem),ivoi=1,nbnn,1)C            write(*,*) 'coeff(1) =',(matcoe.mat(1,ivoi),ivoi=1,nbnn,1)C            write(*,*) 'coeff(2) =',(matcoe.mat(2,ivoi),ivoi=1,nbnn,1)C            write(*,*) 'coeff(3) =',(matcoe.mat(3,ivoi),ivoi=1,nbnn,1)C            if(idim.eq.3) write(*,*) 'coeff(4)=',C     &           (matcoe.mat(4,ivoi),ivoi=1,nbnn,1)C            xv=0.0D0C            yv=0.0D0C            zv=0.0D0C            xc=0.0D0C            do ivoi=1,nbnn,1C               xv=xv+matcoe.mat(1,ivoi)C               yv=yv+matcoe.mat(2,ivoi)C               zv=zv+matcoe.mat(3,ivoi)C               if(idim.eq.3) xc=xc+matcoe.mat(4,ivoi)C            enddoC            write(*,*) 'sum_1=',xvC            write(*,*) 'sum_2=',yvC            write(*,*) 'sum_3=',zvC            if(idim.eq.3) write(*,*) 'sum_4=',xcCC            SEGDES MATCOE         ENDDO         SEGDES IPT1      ENDDOC      SEGSUP MLECEN      SEGSUP MLEFAC      SEGSUP MLEBCC      SEGDES MPNORM      IF(NSOU .GT. 1) SEGDES MELEME      SEGSUP MLESOUC      SEGDES MLECOEC      SEGSUP MATA      SEGSUP MATU      SEGSUP MATV      SEGSUP MLRSIG      SEGSUP MLRTRA      SEGSUP MLRBC      SEGSUP MTAA      SEGSUP MINTAA      SEGSUP MLRIN1      SEGSUP MLRIN2      SEGSUP MLRIN3CC      write(ioimp,*) 'FINITO'C      stopC 9999 CONTINUE      RETURN      END

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