Sources Esope | Cast3M

rlenco
C RLENCO    SOURCE    OF166741  24/12/13    21:17:23     12097          
      SUBROUTINE RLENCO(MELSOM,MELCEN,MELLI1,MELLI2,INORM,MLEPOI,MLECOE)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  RLENCO
C
C DESCRIPTION       :  Appelle par GRADI2
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  A. BECCANTINI
C
C************************************************************************
C
C     This subroutine computes the coefficients which allow to perform
C     a linear exact reconstruction of each 'sommet' (vertex) value
C     Indeed, given the i-th sommet (V=MELSOM.NUM(1,i)), given the list
C     of its neighbors MLEPOI.LECT(i), we have to compute the matrix of
C     coefficients a(i,j) such that
C
C     VAL(MELSOM.NUM(1,i)) =  \sum_{neig1} a(i,neig1) * VAL(neig1) +
C                             \sum_{neig2} a(i,neig2) * VALNOR(neig2)
C
C     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) is
C     the value of grad(VAL).n is the j-th neighbor.
C
C
C     To compute these coefficients, we impose that VAL is linear with
C     respect to the coordinates of V and its neighbor. Then, since
C     there are less unknowns than equations, we use a least square
C     method. If the neighbor is a 'CENTRE' point or a 'Dirichlet
C     boundary condition' point, we impose that
C
C     A + B (x_{neig} - x_V) + C (y_{neig} - y_V) + D (z_{neig} - z_V)
C     =  VAL(neig)
C
C     If the neighbor is a 'von Neumann boundary condition' point,
C     we impose that
C
C     |r_{neig} - r_V| (B nx_{neig} + C ny_{neig} + D nz_{neig}) =
C     VALNOR(neig) * |r_{neig} - r_V|
C
c     r = (x,y,z)
c     n = (nx, ny, nz) is the normal to the surface
C
C     To determine A, B, C, D we have to solve
C
C     MATA . tran(A, B, C, D) = tran(VAL(1), VAL(2), ... , VALNOR(j) /
C                           (r_{j} - r_V)
C
C     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          ...
C
C    To determine a(i,neig) we have to solve a linear system (for each
C    neig). If the neighbor is a 'CENTRE' point or a 'Dirichlet boundary
C    condition' point we solve
C
C    MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) = e_{neig}
C
C    with e_{neig} = tran (0,0,...,0,1,0,...) (1 in the neig position)
C
C    If the neighbor is a 'von Neumann boundary condition' points,
C    we solve
C
C    MATA . tran(a(i,neig), b(i,neig), c(i,neig), d(i,neig)) =
C       e_{neig}  (r_{neig} - r_V)
C
C    Since each linear system involves MATA, we have to "pseudoinvert"
C    it.
C
C**** Inputs:
C
C     MELSOM = the 'SOMMET' meleme
C     MELCEN = the 'CENTRE' meleme
C     MELLI1 = the support of Dirichlet boundary conditions
C     MELLI2 = the support of von Neumann boundary conditions
C     INOMR  = CHAMPOINT des normales aux faces
C     MLEPOI = list of integers.
C              MLEPOI.LECT(i) points to the list neighbors of
C              MELSOM.NUM(1,i)
C
C**** Output:
C
C     MLECOE =  list of integers.
C
C               MLECOE.LECT(i) points to the MLREEL MLRCOE which
C               contain the a(i,neig)
C
      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 PPARAM
-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.MLENTI
C
C**** We create the MLENTI for the centers
C
      CALL KRIPAD(MELCEN,MLECEN)
      IF(IERR .NE. 0)GOTO 9999
C     En KRIPAD
C     SEGACT MELCEN
C     SEGINI MLECEN
C
C**** We create the MLENTI for the BC
C
      CALL KRIPAD(MELLI1,MLELI1)
      IF(IERR .NE. 0)GOTO 9999
      CALL KRIPAD(MELLI2,MLELI2)
      IF(IERR .NE. 0)GOTO 9999
C     SEGACT MELLI1
C     SEGINI MLELI1
C     SEGACT MELLI2
C     SEGINI MLELI2
C
      SEGACT MELSOM
      NSOMM=MELSOM.NUM(/2)
      JG=NSOMM
      SEGINI MLECOE
      SEGACT MLEPOI
C
C**** Le MPOVAL des normales
C
      CALL LICHT(INORM,MPNORM,TYPE,MELFAC)
C     SEGACT MPNORM
C
C**** We create the MLENTI for the MELFAC
C
      CALL KRIPAD(MELFAC,MLEFAC)
      IF(IERR .NE. 0)GOTO 9999
C     SEGACT MELFAC
C     SEGINI MLEFAC
C
C**** Loop on the sommet to compute NVMAX
C     (maximum number of neighbors)
C
      NVMAX=0
      DO INOEU=1,NSOMM,1
C
C******* The neighbors of NGS
C
         MLENTI=MLEPOI.LECT(INOEU)
         SEGACT MLENTI
         NVOI=MLENTI.LECT(/1)
         NVMAX=MAX(NVMAX,NVOI)
      ENDDO
C
C     MATA = matrix to "pseudoinvert" (NVMAX,IDIM+1)
C     MATU = matrix of the singular right eigenvectors of MATA
C            (NVMAX,IDIM+1)
C     MATV = matrix of the singular left eigenvectors of MATA
C            (IDIM+1,IDIM+1)
C            But in invsvd.eso, MATV dimensions are NVMAX,IDIM+1
C     MLRSIG = singular values of MATA (IDIM+1)
C     MLRB   = vector (NVMAX)
C              MLRB.PROG(j) = 1 if j is a center point or a 'Dirichlet
C                             boundary condition' point
C              MLRB.PROG(j) = (r_j - r_V) id j is a 'von Neumann
C                              boundary condition point.
C
C     N.B. MATA = MATU MATSIG t(MATV)
C          If MATA is non singular,
C          inv(MATA) = MATV inv(MATSIG) t(MATU)
C
C     MLRTRA temporary vector of invsvd.eso
C     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=N2
C
C     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 MLRIN3
C
C**** Loop on the sommet to compute the coefficient
C
 
      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)
C
C******* The neighbors of NGS
C
         MLENTI=MLEPOI.LECT(INOEU)
         NVOI=MLENTI.LECT(/1)
C
C******* The matrix where we store the coefficients
C
         JG=NVOI
         SEGINI MLRCOE
         MLECOE.LECT(INOEU)=MLRCOE
C
C******* Loop involving the neighbors
C        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))THEN
C
C********** NGV is a 'centre' point or a Dirichlet point
C
               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)THEN
C
C**************** NO B.C. specified on this point
C
                  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
         ENDDO
C
C        TEST
C         do ivoi=1,nvoi,1
C            write(*,*)
C     &           'mata.mat(',ivoi,')=',(mata.mat(ivoi,i1),i1=1,nvmin,1)
C            write(*,*) 'b(',ivoi,')=',mlrb.prog(ivoi)
C         enddo
C
C
C********** Now we have to invert this matrix
C
         LOGSVD=.TRUE.
         CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT,
     &        MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD,
     &        MLRTRA.PROG)
         IF(IERSVD.NE.0)THEN
C
C************* SVD decomposition of the matrix does not work
C
            LOGSVD=.FALSE.
         ELSE
C
C************ We check the condition number of MATA
C
            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
         ENDIF
C
C        TEST
C         write(*,*) 'LOGSVD=.FALSE.'
C         LOGSVD=.FALSE.
C
         IF(LOGSVD)THEN
C
C********** INVSVD worked
C
C           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,1
C                    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 0
C           Opération malvenue. Résultat douteux
            CALL ERREUR(22)
C
C********** INVSVD does not worked
C           For each neighbor k we have to compute the solution
C           of
C
C           t(MATA) MATA x = t(MATA) * b
C
C           where b= \sum_l e_l \delta(k,l) = e_k
C
C           To do that, we compute
C
C           X_0 = [inve(t(MATA) MATA)] [t(MATA) * b]
C
C           X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
C                  [t(MATA) MATA] X_0]
C
C
C********** 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
            ENDDO
C
C********** 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 2
C              Tache impossible. Probablement données erronées
               CALL ERREUR(26)
               GOTO 9999
            ENDIF
C
C********** We complete MTAA and MINTAA
C
            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
            ENDDO
C
            DO IVOI=1,NVOI,1
C
C************* We compute [t(MATA) . b]  and we store it in MLRIN1.PROG
C
               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
               ENDDO
C
C************* 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
               ENDDO
C
C************ X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b -
C                  [t(MATA) MATA] X_0]
C
C             [t(MATA) MATA] X_0 into MLRIN3.PROG
C
               DO I2=1,NVMIN,1
                  DO I1=1,NVMIN,1
                     MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+
     &                    (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2))
                  ENDDO
               ENDDO
C
C************* Now we have
C              [t(MATA) . b]  in MLRIN1.PROG
C              X_0(i1)        in MLRIN2.PROG(I1)
C              [t(MATA) MATA] X_0 in MLRIN3.PROG
C
C              X_1(i1) in MLRCOE.MAT(i1,IVOI)
C
               DO I1=1,1,1
C                 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
         ENDIF
CC
CC        TEST
C         write(*,*) 'ngs =', NGS
C         write(*,*) 'invide',LOGSVD
C         write(*,*) 'nvois =',(mlenti.lect(ivoi),ivoi=1,nvoi,1)
C         write(*,*) 'coeff =',(mlrcoe.prog(ivoi),ivoi=1,nvoi,1)
C         dx=0.0D0
C         do ivoi=1,nvoi,1
C            dx=dx+mlrcoe.prog(ivoi)
C         enddo
C         write(*,*) 'sum=',dx
         MLENTI=MLEPOI.LECT(INOEU)
         SEGDES MLENTI
         MLRCOE=MLECOE.LECT(INOEU)
         SEGDES MLRCOE
      ENDDO
C
      SEGDES MELCEN
      SEGSUP MLECEN
C
      IF(MELLI1 .NE. 0) SEGDES MELLI1
      IF(MELLI2 .NE. 0) SEGDES MELLI2
      SEGSUP MLELI1
      SEGSUP MLELI2
C
      SEGDES MELSOM
C
      SEGDES MLECOE
      SEGDES MLEPOI
C
      SEGSUP MATA
      SEGSUP MATU
      SEGSUP MATV
      SEGSUP MLRSIG
      SEGSUP MLRTRA
      SEGSUP MLRB
C
      SEGSUP MTAA
      SEGSUP MINTAA
      SEGSUP MLRIN1
      SEGSUP MLRIN2
      SEGSUP MLRIN3
C
      SEGDES MPNORM
      SEGDES MELFAC
      SEGSUP MLEFAC
C
 9999 CONTINUE
      RETURN
      END