Esope source files

   1 : C IPLNU1    SOURCE    PV        09/03/12    21:25:29     6325           
   2 :       SUBROUTINE IPLNU1(INUA)
   3 : c
   4 : c fonction
   5 : c     cette sous routine interpole une fonction de n variables dont
   6 : c     certains couples (x,f(x)) ont été stockés dans un objet de type
   7 : c     nuage
   8 : c
   9 : c variables
  10 : c     inua pointeur sur l'objet de type nuage
  11 : c     iadd contient la correspondance entre les numeros de mtr2 de
  12 : c          1 à nn+mm et ceux du nuage.
  13 : c     mtr2 contient les valeurs des variables maitres dans xi et les
  14 : c          valeurs esclaves dans yi
  15 : c
  16 : c appelé par interpol qui vient de l'opérateur ipol
  17 : c
  18 : c auteur a de gayffier
  19 : c date  05/07/94
  20 : c
  21 : c langage esope+fortran 77
  22 : c
  23 :       IMPLICIT INTEGER(I-N)
  24 :       IMPLICIT    REAL*8(A-H,O-Z)
  25 :       PARAMETER (TINY=1.D-20)
  26 :       LOGICAL FLAG
  27 : -INC SMNUAGE
  28 : -INC SMCHAML
  29 : -INC CCOPTIO
  30 : -INC SMCHPOI
  31 :       SEGMENT MTR1
  32 :          INTEGER IADD(NNMM)
  33 :       ENDSEGMENT
  34 :       SEGMENT MTR2
  35 :          REAL*8  XI(NN),YI(MM)
  36 :       ENDSEGMENT
  37 :       SEGMENT MAX1
  38 :          REAL*8  MAXV(NNMM)
  39 :       ENDSEGMENT
  40 :       CHARACTER*5 CLE(2)
  41 :       DATA CLE /'GAUSS','RATIO'/
  42 : c
  43 : c   option pour la fonction de ponderation
  44 : c
  44 :       CALL LIRMOT(CLE,2,IVAL,0)
  45 :       IF (IVAL .EQ. 0) IVAL = 1
  46 : c
  47 :       MNUAGE = INUA
  48 :       SEGACT ,MNUAGE
  49 : 
  50 :       NNMM = NUANOM(/2)
  51 :       SEGINI ,MAX1
  52 : 
  53 : c    recherche du maximum de chaque paramètre
  54 : 
  55 :         DO 20 K=1,NNMM
  56 :            MAXV(K) = 0
  57 :            NUAVFL = NUAPOI(K)
  58 :            SEGACT ,NUAVFL
  59 :            NBCOUP = NUAFLO(/1)
  60 :            MAXV(K)=TINY
  61 :            DO 10 J=1,NBCOUP
  62 :               IF(MAXV(K).LT.ABS(NUAFLO(J))) THEN
  63 :                  MAXV(K) = ABS(NUAFLO(J))
  64 :               ENDIF
  65 :   10    CONTINUE
  66 :   20    CONTINUE
  67 : c
  67 :       CALL LIROBJ('CHPOINT',IPOI,0,IRETOU)
  68 :       IF ( IRETOU .EQ. 1 ) GOTO 210
  68 :       CALL LIROBJ('MCHAML',ICHML,0,IRETOU)
  69 :       IF ( IRETOU .EQ. 0 ) THEN
  70 : c   pas d'argument correct trouvé: erreur
  70 :         CALL QUETYP(MOTERR(1:8),0,IRETOU)
  71 :         IF(IRETOU.NE.0) THEN
  71 :             CALL ERREUR(39)
  72 :         ELSE
  72 :             CALL ERREUR(533)
  73 :         ENDIF
  74 :         RETURN
  75 :       ENDIF
  76 : c--------------------cas d'un mchml---------------------
  77 : c
  78 : c
  79 :       MCHELM = ICHML
  80 :       SEGACT MCHELM
  81 : c      creation du chapeau du chamelem resultat
  82 :       SEGINI ,MCHEL1 = MCHELM
  83 : c
  84 :       N1 = ICHAML(/1)   
  85 : c
  86 : c boucle sur les chamelems elementaires
  87 : c
  88 :       DO 200 I=1,N1
  89 :         MCHAML = ICHAML(I)
  90 :         SEGACT MCHAML
  91 : c      calcul des dimensions nn et pp
  92 :         NNMM = 0
  93 :         SEGINI MTR1
  94 :         ITR1 = MTR1
  95 :         NN = NOMCHE(/2)
  96 :         MM = NUANOM(/2) - NN
  97 : c      initialisation du nouveau sous chamelem
  98 :         N2 = MM
  99 :         SEGINI ,MCHAM1
 100 :         MCHEL1.ICHAML(I) = MCHAM1
 101 : c
 102 : c        on verifie que les noms des composantes du chamelem sont
 103 : c        dans le nuage et on rempli iadd
 104 : c
 105 :         DO 110 J=1,NOMCHE(/2)
 106 :           IF ( TYPCHE(J) .NE. 'REAL*8' ) THEN
 107 :            MOTERR(1:8) = NOMCHE(J)
 108 :            MOTERR(9:13) = 'CHAMP'
 108 :            CALL ERREUR(681)
 109 : *
 110 : *    'la composante',nomche(j),'du chamelem''n est pas scalaire'*
 111 :            SEGDES MCHAML,MCHELM,MNUAGE
 112 :            SEGSUP MTR1,MCHAM1,MCHEL1
 113 :            RETURN
 114 :           ENDIF
 115 : *
 116 :           DO 100 K=1,NUANOM(/2)
 117 :             IF (NUANOM(K) .EQ. NOMCHE(J)) THEN
 118 :                IADD(**)=K
 119 :                NNMM = IADD(/1)
 120 :              ENDIF
 121 :  100      CONTINUE
 122 :           IF ( NNMM .NE. J ) THEN
 123 : c        une des composantes du champ n'est pas un composante du nuage
 123 :            CALL ERREUR(682)
 124 :            SEGDES MCHAML,MCHELM,MNUAGE
 125 :            SEGSUP MTR1,MCHAM1,MCHEL1
 126 : *           return
 127 :           ENDIF
 128 :   110     CONTINUE
 129 : c
 130 : c         on recupere les composantes du nuage qui ne sont pas celles
 131 : c         du chamelem
 132 : c
 133 :         DO 130 J=1,NUANOM(/2)
 134 :           IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
 135 :             MOTERR(1:8) = NUANOM(J)
 136 :             MOTERR(9:13) = 'NUAGE'
 136 :             CALL ERREUR(681)
 137 : *      'la composante', nuanom(j), 'du nuage', 'n est pas scalaire'
 138 :             SEGDES MCHAML,MCHELM,MNUAGE
 139 :             SEGSUP MTR1,MCHAM1,MCHEL1
 140 :             RETURN
 141 :           ENDIF
 142 :           FLAG = .TRUE.
 143 :           DO 120 K=1,NOMCHE(/2)
 144 :             IF (NUANOM(J) .EQ. NOMCHE(K)) THEN
 145 :                FLAG = .FALSE.
 146 :             ENDIF
 147 :   120     CONTINUE
 148 :           IF (FLAG) THEN
 149 : c          la composante du nuage n'est pas dans le chamelem
 150 :             IADD(**) = J
 151 :             NNMM = IADD(/1)
 152 :            ENDIF
 153 :   130   CONTINUE
 154 :         IF ( NNMM .NE. NUANOM(/2)) THEN
 154 :           CALL ERREUR(682)
 155 : *        'incoherence entre les composantes du nuage et du champ'
 156 :           SEGDES MCHAML,MCHELM,MNUAGE
 157 :           SEGSUP MTR1,MCHAM1,MCHEL1
 158 :           RETURN
 159 :         ENDIF
 160 : c
 161 : c      on rempli les types et le nom du nouveau sous chamelem
 162 : c
 163 :         DO 135 J=1,MM
 164 :           MCHAM1.TYPCHE(J)='REAL*8'
 165 :           MCHAM1.NOMCHE(J)=NUANOM(IADD(NN+J))
 166 :   135   CONTINUE
 167 : c
 168 : c   activation de tous les melval
 169 : c
 170 :          DO 136 J=1,NN
 171 :           MELVAL = IELVAL(J)
 172 :           SEGACT ,MELVAL
 173 :   136   CONTINUE
 174 : c
 175 : c    on recupere le nombre d'élément et le nombre de points par element
 176 : c
 177 :         N1PTEL = VELCHE(/1)
 178 :         N1EL = VELCHE(/2)
 179 :         N2PTEL = 0
 180 :         N2EL = 0
 181 : c
 182 :         DO 138 J=1,MM
 183 :           SEGINI MELVAL
 184 :           MCHAM1.IELVAL(J) = MELVAL
 185 :   138   CONTINUE
 186 : c
 187 : c    boucle sur chaque point de gauss
 188 : c
 189 :         DO 170 J=1,N1EL
 190 :         DO 160 K=1,N1PTEL
 191 :            SEGINI ,MTR2
 192 :            ITR2 = MTR2
 193 : c
 194 : c        boucle sur les composantes pour remplir le tableau xi
 195 : c        qui contient les valeurs des variables maitres
 196 :            DO 140 L=1,NN
 197 :              MELVAL = IELVAL(L)
 198 :              XI(L) = VELCHE(K,J)
 199 :   140      CONTINUE
 200 : c
 201 : c        on interpole le nuage pour les valeurs xi
 202 : c
 202 :            CALL IPLNU2(INUA,ITR1,ITR2,MAX1,IVAL)
 203 :            IF ( IERR .NE. 0 ) THEN
 204 :                SEGDES MCHAML,MCHELM
 205 :                DO 145 L=1,MCHAM1.IELVAL(/1)
 206 :                    MELVAL = MCHAM1.IELVAL(L)
 207 :                    SEGSUP MELVAL
 208 :   145          CONTINUE
 209 :                SEGSUP MCHAM1,MCHEL1,MTR1,MTR2
 210 :                RETURN
 211 :            ENDIF
 212 : c
 213 : c        remplissage des résultats stockés dasn yi
 214 : c
 215 :            DO 150 L=1,MM
 216 :              MELVAL = MCHAM1.IELVAL(L)
 217 :              VELCHE(K,J) = YI(L)
 218 :   150      CONTINUE
 219 :            SEGSUP ,MTR2
 220 : c
 221 :   160    CONTINUE
 222 :   170    CONTINUE
 223 : c
 224 : c   on desactive tout les segments melval
 225 :          DO 180 J=1,NN
 226 :           MELVAL = IELVAL(J)
 227 :           SEGDES ,MELVAL
 228 :   180   CONTINUE
 229 : c
 230 :         DO 190 J=1,MM
 231 :           MELVAL = MCHAM1.IELVAL(J)
 232 :           SEGDES ,MELVAL
 233 :   190   CONTINUE
 234 :          SEGDES ,MCHAML
 235 :          SEGDES ,MCHAM1
 236 :          SEGSUP ,MTR1
 237 :   200 CONTINUE
 238 : c
 239 : c fin de la boucle sur les sous champs mchaml
 240 :        ICH1 = MCHEL1
 241 :        SEGDES MNUAGE,MCHELM,MCHEL1
 242 :        SEGSUP MAX1
 242 :        CALL ECROBJ('MCHAML  ',ICH1)
 243 :        RETURN
 244 : c
 245 : c---------------cas du chpoint-----------------------
 246 :   210 CONTINUE
 247 : c
 248 :       MCHPOI = IPOI
 249 :       SEGACT MCHPOI
 250 : c
 251 :       SEGINI ,MCHPO1 = MCHPOI
 252 :       NSOUPO = IPCHP(/1)
 253 : c
 254 : c boucle sur les chpoints elementaire
 255 :       DO 280 I=1,NSOUPO
 256 :         MSOUPO = IPCHP(I)
 257 :         SEGACT MSOUPO
 258 : c
 259 : c  calcul de nn et mm
 260 : c
 261 :          NNMM=0
 262 :          SEGINI MTR1
 263 :          ITR1 = MTR1
 264 :          NN = NOCOMP(/2)
 265 :          MM = NUANOM(/2) - NN
 266 : c   initialisation du nouveau soupo
 267 :          NC = MM
 268 :          SEGINI MSOUP1
 269 :          MCHPO1.IPCHP(I)=MSOUP1
 270 :          MSOUP1.IGEOC = IGEOC
 271 : c
 272 : c  on verifie que les noms des composantes du chpoint sont
 273 : c        dans le nuage et on rempli iadd
 274 : c
 275 :         DO 230 J=1,NOCOMP(/2)
 276 :           DO 220 K=1,NUANOM(/2)
 277 :             IF (NUANOM(K) .EQ. NOCOMP(J)) THEN
 278 :                IADD(**)=K
 279 :                NNMM = IADD(/1)
 280 :              ENDIF
 281 :  220      CONTINUE
 282 :           IF ( NNMM .NE. J ) THEN
 283 : c        une des composantes du champ n'est pas un composante du nuage
 283 :            CALL ERREUR(682)
 284 : *          'incoherence entre les composantes du nuage et du champ'
 285 :            SEGDES MCHPOI,MSOUPO,MNUAGE
 286 :            SEGSUP MTR1,MCHPO1,MSOUP1
 287 :            RETURN
 288 :           ENDIF
 289 :   230     CONTINUE
 290 : c
 291 : c  on recupere les composantes du nuage qui ne sont pas dans le chpoint
 292 : c
 293 :         DO 250 J=1,NUANOM(/2)
 294 :           IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
 295 : *            la composante du nuanom(j) du nuage n'estpas scalaire
 296 :            MOTERR(1:8) = NUANOM(J)
 297 :            MOTERR(9:13) = 'NUAGE'
 297 :            CALL ERREUR(681)
 298 :            SEGDES MCHPOI,MSOUPO,MNUAGE
 299 :            SEGSUP MTR1,MCHPO1,MSOUP1
 300 :            RETURN
 301 :          ENDIF
 302 :           FLAG = .TRUE.
 303 :           DO 240 K=1,NOCOMP(/2)
 304 :             IF (NUANOM(J) .EQ. NOCOMP(K)) THEN
 305 :                FLAG = .FALSE.
 306 :             ENDIF
 307 :   240     CONTINUE
 308 :           IF (FLAG) THEN
 309 : c          la composante du nuage n'est pas dans le chamelem
 310 :             IADD(**) = J
 311 :             NNMM = NNMM +1
 312 :            ENDIF
 313 :   250   CONTINUE
 314 :         IF ( NNMM .NE. NUANOM(/2)) THEN
 315 :           PRINT *, 'INCOHERENCE DANS LES COMPOSANTES'
 316 :         ENDIF
 317 : c
 318 : c   on rempli le noms des composantes du nouveau champ
 319 : c
 320 :         DO 255 J=1,MM
 321 :           MSOUP1.NOCOMP(J)=NUANOM(IADD(J+NN))
 322 :   255   CONTINUE
 323 :         MPOVAL = IPOVAL
 324 :         SEGACT MPOVAL
 325 : c      nouveau segment mpoval
 326 :         N = VPOCHA(/1)
 327 :         NC = MM
 328 :         SEGINI MPOVA1
 329 :         MSOUP1.IPOVAL = MPOVA1
 330 : c
 331 : c    boucle sur les points
 332 : c
 333 :         DO 270 J=1,N
 334 :           SEGINI MTR2
 335 :           ITR2 = MTR2
 336 :           DO 260 K=1,NN
 337 :             XI(K)=VPOCHA(J,K)
 338 :   260     CONTINUE
 339 : c
 340 : c   interpolation
 341 : c
 341 :           CALL IPLNU2(INUA,ITR1,ITR2,MAX1,IVAL)
 342 : c         en cas d'erreur on sort proprement
 343 :           IF ( IERR .NE. 0 ) THEN
 344 :                SEGDES MCHPOI,MSOUPO
 345 :                SEGSUP MCHPO1,MSOUP1,MPOVAL,MTR1,MTR2,MAX1
 346 :                RETURN
 347 :           ENDIF
 348 : c
 349 :           DO 265 K=1,MM
 350 :             MPOVA1.VPOCHA(J,K)=YI(K)
 351 :   265     CONTINUE
 352 :           SEGSUP ,MTR2
 353 : c
 354 :   270     CONTINUE
 355 : c   desactivation des segments
 356 :           SEGSUP,MTR1
 357 :           SEGDES ,MPOVAL
 358 :           SEGDES ,MSOUPO
 359 :           SEGDES ,MPOVA1,MSOUP1
 360 : c
 361 :   280 CONTINUE
 362 : c
 363 :       IP1 = MCHPO1
 364 :       SEGDES ,MCHPOI
 365 :       SEGDES ,MCHPO1
 366 :       SEGDES ,MNUAGE
 367 :       SEGSUP ,MAX1
 367 :       CALL ECROBJ('CHPOINT ',IP1)
 368 : c
 369 :       RETURN
 370 :       END
 371 :  
 372 :  
 373 :  
 374 :