Sources Esope | Cast3M

iplnu2
C IPLNU2    SOURCE    BP208322  21/01/28    21:15:11     10868          
 
        SUBROUTINE IPLNU2(INUA,ITR1,XI,YI,MAXI1,IVAL,
     &                    N,INDJ,NN,MM,NNMM,XP,XELIM)
C
C Fonction
C     Cette sous routine interpole une fonction de n variables dont
C     certains couples (x,f(x)) on ete stockes dans un objet de type
C     nuage.
C     2 techniques sont disponibles :
C     - Pour IVAL = 1 ou 2, on utilise une methode des elements diffus.
C       En chaque point on cherche la forme lineaire qui minimise la distance
C       ponderee de l'hyperplan plan au point du nuage. La valeur obtenue est
C       celle de la forme lineaire au point considere.
C     - Pour IVAL = 3, on utilise une interpolation exacte par ponderation
C       inverse a la distance. Le poids attribue a chaque point du nuage est
C       en 1/D**XP ou D est sa distance au point cible et XP le parametre
C       puissance fourni par l'utilisateur.
C Variables
C     INUA pointeur sur l'objet de type nuage
C     NNMM nombe de composantes du nuage
C     NN nombre de composantes connues du point cible
C     MM nombre de composantes a calculer du point cible
C     NBCOUP nombe de points (cad de NNMM-uplets) du nuage
C     ITR1.IADD contient la correspondance entre les numeros des composantes
C               du CHPOINT/MCHAML et celles du nuage :
C          - La composante connue I (entre 1 et NN) du champ est a la position
C            IADD(I) dans le nuage
C          - La composante a calculer J (entre 1 et MM) est a la position
C            IADD(J) dans le nuage
C     MAXI1.MAXV contient les max des composantes du nuage (utilise pour normer
C                le calcul des distances)
C     X tableau contenant les valeurs des NN composantes connues du CHPOINT/MCHAML
C       X(i,j) est la valeur de la composante j au i-eme point support du
C       champ, j est compris entre 1 et NN, la parallelisation est faite sur i
C     Y tableau contenant les valeurs des MM composantes interpolees, c'est
C       le resultat de l'interpolation
C       Y(i,j) est la valeur de la composante j au i-eme point support du
C       champ, j est compris entre 1 et MM, la parallelisation est faite sur i
C     IVAL est l'entier indiquant l'option choisie
C       IVAL=1 : ponderation par fonction gaussienne  EXP(-Z**2) (valeur par defaut)
C       IVAL=2 : ponderation par fonction rationnelle 1/(1+Z)
C       IVAL=3 : interpolation exacte par ponderation inverse a la distance
C
C Auteur A De Gayffier
C Date  05/07/94
C
C Langage esope+fortran 77
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT    REAL*8(A-H,O-Z)
      PARAMETER ( GRAND = 1.D50)
 
-INC PPARAM
-INC CCOPTIO
-INC SMNUAGE
      SEGMENT MTR1
        INTEGER IADD(NNMM)
      ENDSEGMENT
 
      REAL*8  XI,YI
      DIMENSION XI(N,*)
      DIMENSION YI(N,*)
 
C     Ces tableaux contienent la matrice du syteme lineaire
C     A contient les termes et INDX contient les indices de permutation
C     une fois la decomposition effectuee
      REAL*8 A(NN+1,NN+1),D
      INTEGER INDX(NN+1)
 
      REAL*8 B(NN+1,MM)
 
      REAL*8 DX(NN+1)
 
      REAL*8 DIMINI(NN+1)
      INTEGER NNBCP(NN+1)
 
      SEGMENT MAXI1
        REAL*8 MAXV(NNMM)
      ENDSEGMENT
 
C
      MNUAGE = INUA
      NUAVFL = NUAPOI(1)
      NBCOUP = NUAFLO(/1)
 
C (FDP 2016) Nouvelle option 'PID' : interpolation exacte par ponderation
C inverse a la distance
      IF (IVAL.EQ.3) THEN
        CALL IPPID(INUA,NBCOUP,ITR1,XI,YI,N,INDJ,NN,MM,NNMM,XP,XELIM)
        RETURN
      ENDIF
 
 
C      INITIALISATIONS A ZERO
       DO INDICE=1,NN+1
         INDX(INDICE)  =0
         DIMINI(INDICE)=0.D0
         NNBCP(INDICE) =0
         DX(INDICE)    =0
 
         DO JNDICE=1,NN+1
           A(INDICE,JNDICE)=0.D0
         ENDDO
 
         DO JNDICE=1,MM
           B(INDICE,JNDICE)=0.D0
         ENDDO
       ENDDO
 
C
      IF (NBCOUP .EQ. 0) THEN
C       nuage vide operation impossible
        CALL ERREUR(680)
        RETURN
      ENDIF
C
      MTR1 = ITR1
C    initialisation du tableau pour stocker les distances
C    on place des valeurs tres grandes
      DO 1 I=1,NN+1
        DIMINI(I)= GRAND
 1    CONTINUE
C
C   calcul de la distance cararcteristique du nuage
C
      DO 60 I=1,NBCOUP
        SUM = 0.D0
        DO 30 J=1,NN
           NUAVFL = NUAPOI(IADD(J))
C        on divise les distances par les maximums de chaque composante
           DUM = (XI(INDJ,J)-NUAFLO(I))/MAXV(IADD(J))
           SUM = SUM + DUM*DUM
 30     CONTINUE
        DIST = SQRT(SUM)
C   stockage des plus petites valeurs
        DO 50 J=1,NN+1
           IF ( DIST .LT. DIMINI(J) )THEN
             DO 40 K=NN+1,J+1,-1
               DIMINI(K)=DIMINI(K-1)
               NNBCP(K)=NNBCP(K-1)
 40          CONTINUE
             DIMINI(J)=DIST
             NNBCP(J)=I
             GOTO 60
           ENDIF
 50     CONTINUE
 60   CONTINUE
C
      DCARA =0.D0
      DO 65 J=1,NN+1
        DCARA = DCARA + DIMINI(J)
  65  CONTINUE
      DCARA = DCARA / (NN+1.D0)
C
C   affichage de mise au point avec opti impi 9022 ;
C
      IF (IIMPI .EQ. 9022) WRITE(6,1001) DCARA
      IF (IIMPI .EQ. 9022) THEN
           WRITE(ioimp,1002) (NN+1)
           write(ioimp,1005) (DIMINI(I),I=1,NN+1)
C               boucle sur les composantes
           DO 67 I=1,NN
             NUAVFL = NUAPOI(IADD(I))
             WRITE(ioimp,1004) NUANOM(IADD(I))
             write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
             WRITE(ioimp,*)
   67      CONTINUE
           DO 68 I=NN+1,NNMM
             NUAVFL = NUAPOI(IADD(I))
             WRITE(ioimp,1004) NUANOM(IADD(I))
             write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
             WRITE(ioimp,*)
   68      CONTINUE
      ENDIF
 1001 FORMAT('La distance caracteristique vaut ',G10.4)
 1002 FORMAT('Distance et coordonees des',1X,I2,1X,'points
     &        les plus proches sont')
 1003 FORMAT(30(G10.4,1X))
 1004 FORMAT('Composante ',A)
 1005 FORMAT('Distance au point courant:',/,30(G10.4))
C
C   boucle sur les points du nuage (nbcoup)
C
      DO 110 I=1,NBCOUP
C
C     boucle sur les monomes du polynome d'interpolation (1,x,y,...)
C     remplissage complet de DX (segment MDX)
        DX(1)=1.D0
        SUM = 0.D0
        DO 70 J=2,(NN+1)
           NUAVFL = NUAPOI(IADD(J-1))
C        on divise les distances par les maximums de chaque composante
           DX(J)=(XI(INDJ,J-1)-NUAFLO(I))/MAXV(IADD(J-1))
           SUM = SUM + DX(J)*DX(J)
 70     CONTINUE
        DIST = SQRT(SUM)
C
C  fonction de ponderation de la distance
        Z = DIST/DCARA
        IF (IVAL .EQ. 1 ) THEN
          WDED = EXP(-Z*Z)
        ELSE
          WDED = 1.D0/(1.D0+Z)
        ENDIF
C       remplissage de a et de b
        DO 100 J=1,(NN+1)
          r_z = WDED*DX(J)
          DO 80 K=1,(NN+1)
C*            A(J,K) = A(J,K) + WDED*DX(J) * DX(K)
              A(J,K) = A(J,K) + r_z * DX(K)
 80       CONTINUE
          DO 90 K=1,MM
              NUAVFL = NUAPOI(IADD(NN+K))
C*            B(J,K)=B(J,K) + WDED*DX(J)*NUAFLO(I)
              B(J,K)=B(J,K) + r_z * NUAFLO(I)
 90      CONTINUE
C
 100    CONTINUE
C
 110  CONTINUE
C
C      SEGSUP,MDX
C
C  phase de resolution
C
C  on resoud a.x=b pp fois
C
       CALL IPLNU3(A,D,INDX,NN+1)
 
       DO 120 I=1,MM
         CALL IPLNU4(A,D,INDX,NN+1,B,MM,I)
 120   CONTINUE
C
C   on replace le resultat dans mtr2
C   dans les B(*,I) sont stockees les derivees partielles de la fonction
C
        DO 130 I=1,MM
          YI(INDJ,I)=B(1,I)
 130    CONTINUE
C
C
        RETURN
        END