Sources Esope | Cast3M

ipmul
C IPMUL     SOURCE    FD218221  15/10/05    21:15:06     8658
C-----------------------------------------------------------------------
C NOM         : IPMUL
C DESCRIPTION : Interpolation multi-lineaire dans un NUAGE represantant
C               une grille de valeurs
C LANGAGE     : ESOPE
C AUTEUR      : Francois DI PAOLA
C-----------------------------------------------------------------------
C APPELE PAR       : IPMULI via IPGRIL
C APPELE           :
C-----------------------------------------------------------------------
C ENTREES
C     MVAL1  : Tableau des N coordonnees du point ou l'on souhaite faire
C              l'interpolation
C     MLENT1 : Tableau de (N+1) LISTREELs :
C              - N listes des coordonnees des noeuds de la grille
C              - 1 liste des valeurs de la fonction sur chaque noeud
C              Les N listes de coordonnees sont rangees dans le meme
C              ordre que MVAL1
C SORTIES
C     XVAL   : Valeur de la fonction interpolee au point MVAL1
C-----------------------------------------------------------------------
C VERSION    : v1, 05/10/2015, version initiale
C HISTORIQUE : v1, 05/10/2015, creation
C HISTORIQUE :
C HISTORIQUE :
C-----------------------------------------------------------------------
C Priere de PRENDRE LE TEMPS de completer les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C-----------------------------------------------------------------------
C REMARQUES : - L'interpolation est exacte, c'est-a-dire que l'on
C               retrouve les valeurs de la grille si l'on interpole en
C               un noeud de la grille
C             - La grille peut contenir autant de dimensions que
C               souhaitees
C-----------------------------------------------------------------------
C
      SUBROUTINE IPMUL(VAL,NBC,MLENT1,XVAL)
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCNOYAU
-INC SMLENTI
-INC SMLREEL
C     Le tableau VAL peremt de stocker les coordonnees du point
C     d'interpolation
C     VAL(i) --> valeur de la coordonnee i du point d'interpolation
      REAL*8      VAL(NBC)
C
C     Deux tableaux pour reperer les 2^n noeuds encadrant le point ou l'on
C     veut interpoler
C     IND(i,1) --> indice dans la liste de la coordonne i du noeud situe
C                  en dessous
C     IND(i,2) --> indice dans la liste de la coordonne i du noeud situe
C                  au dessus
C     VAL(i,1) --> valeur de la coordonne i du noeud situe en dessous
C     VAL(i,2) --> valeur de la coordonne i du noeud situe au dessus
      INTEGER     IND(NBC,2)
      REAL*8      FVA(NBC,2)
C
C     Un dernier pour stocker la valeur, en base 2, des 2^n entiers
C     IBIT(i) --> valeur du ieme bit (0 ou 1)
      INTEGER     IBIT(NBC)
C
C     ----------------------------------------------------------------
C     PARTIE 1 : Initialisations
C                - determination de la boite des noeuds encadrant le
C                  point d'interpolation
C                - stockage des coordonnes des noeuds encadrant et des
C                  valeurs de la fonction
C     ----------------------------------------------------------------
C     Recherche, pour chaque coordonnee, des 2 noeuds de la grille
C     encadrant le point d'interpolation -> remplissage de IND et FVA
      WTOT=1D0
      DO I=1,NBC
        X=VAL(I)
C       Recuperation, du LISTREEL de la grille correspondant la
C       coordonnee I
        MLREE1=MLENT1.LECT(I)
        SEGACT,MLREE1
        N1=MLREE1.PROG(/1)
C       Si la liste n'est pas de dimension 2, erreur
        IF (N1.LT.2) THEN
          CALL ERREUR(897)
          RETURN
        ENDIF
C       Si le point est en dehors des noeuds, on le projette sur le
C       noeud le plus proche
        XMIN=MLREE1.PROG(1)
        IF (X.LE.XMIN) THEN
          I1=1
          I2=2
          X1=XMIN
          X2=MLREE1.PROG(I2)
          VAL(I)=XMIN
          GOTO 10
        ENDIF
        XMAX=MLREE1.PROG(N1)
        IF (X.GE.XMAX) THEN
          I1=N1-1
          I2=N1
          X1=MLREE1.PROG(I1)
          X2=XMAX
          VAL(I)=XMAX
          GOTO 10
        ENDIF
C       Parcours de la liste a la recherche des valeurs encadrantes
        DO J=1,(N1-1)
          X1=MLREE1.PROG(J)
          X2=MLREE1.PROG(J+1)
          IF ((X1.LE.X).AND.(X.LE.X2)) THEN
            I1=J
            I2=J+1
            GOTO 10
          ENDIF
        ENDDO
10      CONTINUE
        IND(I,1)=I1
        IND(I,2)=I2
        FVA(I,1)=X1
        FVA(I,2)=X2
C       Calcul du "volume" de la boite de noeuds encadrante
        WTOT=WTOT*(MLREE1.PROG(I2)-MLREE1.PROG(I1))
      ENDDO
C
C     -------------------------------------------------------------
C     PARTIE 2 : Interpolation
C                - recuperation des valeurs de F aux 2^n noeuds encadrant
C                - calcul des ponderations associees a chaque noeud
C     -------------------------------------------------------------
C     Initialisation du resultat de l'interpolation
      XVAL=0D0
      N1=2**NBC
      DO I=0,N1-1
C       Conversion de I en base 2 (calcul des NBC bits)
        II=I
        DO J=1,NBC
          IQ=II/2
          IR=II-(IQ*2)
          IBIT(NBC+1-J)=IR
          II=IQ
        ENDDO
C       Recuperation de la valeur de la fonction au noeud encadrant I
        III=1
        DO J=1,NBC
          JB=IBIT(J)+1
          JL=IND(J,JB)
          JP=1
          IF (J.GT.1) THEN
            DO K=1,J-1
              MLREE1=MLENT1.LECT(K)
              JP=JP*(MLREE1.PROG(/1))
            ENDDO
          ENDIF
          III=III+((JL-1)*JP)
        ENDDO
        MLREE1=MLENT1.LECT(NBC+1)
        FI=MLREE1.PROG(III)
C       Calcul de la ponderation du noeud I, c'est le "volume" du
C       polytope diagonalement oppose au noeud I par rapport au point
C       d'interpolation
        WI=1D0
        DO J=1,NBC
          XJ=VAL(J)
C         sur la coord J, le noeud I est il a gauche ou a droite ?
          JB=IBIT(J)
          IF (JB.EQ.0) THEN
            XJD=FVA(J,2)
            WJ=XJD-XJ
          ELSE
            XJG=FVA(J,1)
            WJ=XJ-XJG
          ENDIF
          WI=WI*WJ
        ENDDO
C       La contribution du noeud I dans l'interpolation est FI*WI
        XVAL=XVAL+(FI*WI)
      ENDDO
C     Enfin, on divise par la taille de la boite encadrante
      XVAL=XVAL/WTOT
C
      RETURN
      END