Télécharger la plaquette Cast3M


fr en

Télécharger iplnu1.eso

Retour à la liste

Numérotation des lignes : 

iplnu1

  1. C IPLNU1    SOURCE    PV090527  25/01/07    14:42:45     12115          
  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 ete stockes 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     NN nombre de composantes du CHPOINT/MCHAML en entree
  12. C     MM nombre de composantes du CHPOINT/MCHAML a calculer
  13. C     NNMM nombe de composantes du nuage (= NN + MM)
  14. C     IADD contient la correspondance entre les numeros des composantes
  15. C          du CHPOINT/MCHAML et celles du nuage :
  16. C          - La composante connue I (entre 1 et NN) du champ est a la position
  17. C            IADD(I) dans le nuage
  18. C          - La composante a calculer J (entre 1 et MM) est a la position
  19. C            IADD(J) dans le nuage
  20. C     MAXV contient les max des composantes du nuage (utilise pour normer
  21. C          le calcul des distances)
  22. C
  23. C Appele par interp qui est le point d'entree de l'operateur IPOL
  24. C
  25. C Auteur A De Gayffier
  26. C Date  05/07/94
  27. C
  28. C Evolution 2015 : nouvelle option GRILL, interpolation dans une grille
  29. C Evolution 2016 : nouvelle option   PID, interpolation par ponderation
  30. C                                         inverse a la distance
  31. C
  32. C Langage esope+fortran 77
  33. C
  34.  
  35.       IMPLICIT INTEGER(I-N)
  36.       IMPLICIT REAL*8(A-H,O-Z)
  37.  
  38.       LOGICAL FLAG
  39.  
  40. -INC SMNUAGE
  41. -INC SMCHAML
  42. -INC CCREEL
  43. -INC SMCOORD
  44.  
  45. -INC PPARAM
  46. -INC CCOPTIO
  47. -INC SMCHPOI
  48. -INC CCASSIS
  49.  
  50. C Introduction d'un COMMON pour la parallelisation
  51. C  Il est a repercuter aussi dans les sources suivantes :
  52. C   - ipln2i.eso
  53.       COMMON/iplnuc/XP1,NBTH1,INUA1,ITR1,IMAX1,IVAL1,N1,NN1,MM1,NNMM1,   
  54.      &              IMPOV1,IMPOV2,
  55.      &              IMCHA1,IMCHA2,N3EL,N3PTEL,
  56.      &              ITR2,XELIM1
  57.  
  58.       EXTERNAL ipln2i
  59.  
  60.  
  61.       SEGMENT MTR1
  62.          INTEGER IADD(NNMM)
  63.       ENDSEGMENT
  64.  
  65.       SEGMENT MTR2
  66.         REAL*8  XI(N,NN)
  67.         REAL*8  YI(N,MM)
  68.       ENDSEGMENT
  69.  
  70.  
  71.       SEGMENT MAXI1
  72.          REAL*8  MAXV(NNMM)
  73.       ENDSEGMENT
  74.  
  75.       CHARACTER*5 CLE(4)
  76.       DATA CLE /'GAUSS','RATIO','PID','GRILL'/
  77.       CHARACTER*4 CLE2(1)
  78.       DATA CLE2 /'ELIM'/
  79.  
  80.       PARAMETER (TINY=1.D-20)
  81.  
  82.  
  83.       INUAG  = INUA
  84.       IMPOV1 = 0
  85.       IMPOV2 = 0
  86.       IMCHA1 = 0
  87.       IMCHA2 = 0
  88.       N3EL   = 0
  89.       N3PTEL = 0
  90.       ITR2   = 0
  91.  
  92.       ITHRD  = 0
  93. C
  94. C   Option pour la fonction de ponderation
  95.       CALL LIRMOT(CLE,4,IVAL,0)
  96.       IF (IVAL .EQ. 0) IVAL = 1
  97. C
  98. C (FDP 2016) Acquisition d'un flottant/entier pour le mot clef PID
  99.       XP=1.D0
  100.       XELIM=XPETIT
  101.       IF (IVAL.EQ.3) THEN
  102.         CALL LIRREE(XP,0,IRETOU)
  103.         IF (XP.LT.0.D0) THEN
  104.           REAERR(1)=REAL(0.D0)
  105.           REAERR(2)=REAL(XP)
  106.           CALL ERREUR(191)
  107.           RETURN
  108.         ENDIF
  109. c       + un flottant supplementaire pour la tolerance a l'ELIMination
  110. c         (i.e. distance en deca de laquelle les noeuds sont confondus)
  111.         CALL LIRMOT(CLE2,1,IVAL2,0)
  112.         IF(IVAL2.EQ.1) THEN
  113.           CALL LIRREE(XELIM,0,IRETOU)
  114.           IF(IRETOU.EQ.0) THEN
  115. c           XP est en fait XELIM ! -> on intervertit
  116.             XELIM=XP
  117.             XP=1.D0
  118.           ELSE          
  119. c           on a lu XP et XELIM (positionnel) -> on teste XELIM
  120.             IF (XELIM.LT.0.D0) THEN
  121.               REAERR(1)=REAL(0.D0)
  122.               REAERR(2)=REAL(XELIM)
  123.               CALL ERREUR(191)
  124.               RETURN
  125.             ENDIF
  126.           ENDIF
  127.         ENDIF
  128.       ENDIF
  129. c       WRITE(IOIMP,*) IVAL,'>>> XP,XELIM=',XP,XELIM
  130. C
  131. C (FDP 2015) Interpolation dans un nuage representant une grille de valeurs
  132. C            La parallelisation est faite dans IPGRIL !
  133.       IF (IVAL.EQ.4) THEN
  134.         CALL IPGRIL(INUA)
  135.         RETURN
  136.       ENDIF
  137.  
  138. C
  139.       MNUAGE = INUA
  140.  
  141.       NNMM = NUANOM(/2)
  142.       SEGINI ,MAXI1
  143.  
  144. C    recherche du maximum de chaque parametre
  145.  
  146.         DO 20 K=1,NNMM
  147.            MAXV(K) = 0
  148.            NUAVFL = NUAPOI(K)
  149.            NBCOUP = NUAFLO(/1)
  150.            MAXV(K)=TINY
  151.            DO 10 J=1,NBCOUP
  152.               IF(MAXV(K).LT.ABS(NUAFLO(J))) THEN
  153.                  MAXV(K) = ABS(NUAFLO(J))
  154.               ENDIF
  155.   10    CONTINUE
  156.   20    CONTINUE
  157. C
  158.       CALL LIROBJ('CHPOINT',IPOI,0,IRETOU)
  159.       IF ( IRETOU .EQ. 1 ) THEN
  160.         CALL ACTOBJ('CHPOINT',IPOI,1)
  161.         GOTO 210
  162.       ENDIF
  163.  
  164.       CALL LIROBJ('MCHAML',ICHML,0,IRETOU)      
  165.       IF ( IRETOU .EQ. 1 ) THEN
  166.         CALL ACTOBJ('MCHAML',ICHML,1)
  167.       ELSE
  168. C       pas d'argument correct trouve: erreur
  169.         CALL QUETYP(MOTERR(1:8),0,IRETOU)
  170.         IF(IRETOU.NE.0) THEN
  171.             CALL ERREUR(39)
  172.         ELSE
  173.             CALL ERREUR(533)
  174.         ENDIF
  175.         RETURN
  176.       ENDIF
  177.  
  178. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  179. C                            Cas d'un MCHAML                           C
  180. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  181.  
  182.       MCHELM = ICHML
  183. C      creation du chapeau du chamelem resultat
  184.       SEGINI ,MCHEL1 = MCHELM
  185. C
  186.       N1 = ICHAML(/1)
  187. C
  188. C boucle sur les chamelems elementaires
  189. C
  190.       DO 200 I=1,N1
  191.         MCHAML = ICHAML(I)
  192. C      calcul des dimensions nn et pp
  193.         NNMM = 0
  194.         SEGINI MTR1
  195.         NN = NOMCHE(/2)
  196.         MM = NUANOM(/2) - NN
  197. C      initialisation du nouveau sous chamelem
  198.         N2 = MM
  199.         SEGINI ,MCHAM1
  200.         MCHEL1.ICHAML(I) = MCHAM1
  201. C
  202. C        on verifie que les noms des composantes du chamelem sont
  203. C        dans le nuage et on rempli iadd
  204. C
  205.         DO 110 J=1,NOMCHE(/2)
  206.           IF ( TYPCHE(J) .NE. 'REAL*8' ) THEN
  207.            MOTERR(1:8) = NOMCHE(J)
  208.            MOTERR(9:13) = 'CHAMP'
  209.            CALL ERREUR(681)
  210. C
  211. C    'la composante',nomche(j),'du chamelem''n est pas scalaire'*
  212.            SEGSUP MTR1,MCHAM1,MCHEL1
  213.            RETURN
  214.           ENDIF
  215. C
  216.           DO 100 K=1,NUANOM(/2)
  217.             IF (NUANOM(K) .EQ. NOMCHE(J)) THEN
  218.                IADD(**)=K
  219.                NNMM = IADD(/1)
  220.              ENDIF
  221.  100      CONTINUE
  222.           IF ( NNMM .NE. J ) THEN
  223. C        une des composantes du champ n'est pas un composante du nuage
  224.            CALL ERREUR(682)
  225.            SEGSUP MTR1,MCHAM1,MCHEL1
  226. C           return
  227.           ENDIF
  228.   110     CONTINUE
  229. C
  230. C         on recupere les composantes du nuage qui ne sont pas celles
  231. C         du chamelem
  232. C
  233.         DO 130 J=1,NUANOM(/2)
  234.           IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
  235.             MOTERR(1:8) = NUANOM(J)
  236.             MOTERR(9:13) = 'NUAGE'
  237.             CALL ERREUR(681)
  238. C      'la composante', nuanom(j), 'du nuage', 'n est pas scalaire'
  239.             SEGSUP MTR1,MCHAM1,MCHEL1
  240.             RETURN
  241.           ENDIF
  242.           FLAG = .TRUE.
  243.           DO 120 K=1,NOMCHE(/2)
  244.             IF (NUANOM(J) .EQ. NOMCHE(K)) THEN
  245.                FLAG = .FALSE.
  246.             ENDIF
  247.   120     CONTINUE
  248.           IF (FLAG) THEN
  249. C          la composante du nuage n'est pas dans le chamelem
  250.             IADD(**) = J
  251.             NNMM = IADD(/1)
  252.            ENDIF
  253.   130   CONTINUE
  254.         IF ( NNMM .NE. NUANOM(/2)) THEN
  255.           CALL ERREUR(682)
  256. C        'incoherence entre les composantes du nuage et du champ'
  257.           SEGSUP MTR1,MCHAM1,MCHEL1
  258.           RETURN
  259.         ENDIF
  260. C
  261. C       Remplissage des types et du nom du nouveau MCHAML
  262.         DO 135 J=1,MM
  263.           MCHAM1.TYPCHE(J)='REAL*8'
  264.           MCHAM1.NOMCHE(J)=NUANOM(IADD(NN+J))
  265.   135   CONTINUE
  266. C
  267. C    on recupere le nombre d'element et le nombre de points par element
  268.         MELVAL = MCHAML.IELVAL(1)
  269. C
  270.         N1PTEL = VELCHE(/1)
  271.         N1EL   = VELCHE(/2)
  272.  
  273. C       Creation des nouveaux MELVAL pour ecrire le resultat
  274.         N2PTEL = 0
  275.         N2EL   = 0
  276.         DO 138 J=1,MM
  277.           SEGINI,MELVA1
  278.           MCHAM1.IELVAL(J) = MELVA1
  279.   138   CONTINUE
  280.  
  281.         N=N1EL*N1PTEL
  282.         NT1= N / (100 * nbthrs)
  283.  
  284. C
  285. C Pour la paralelisation de l'interpolation
  286. C
  287.        ITH = 0
  288.        IF (NBESC .NE. 0) ith=oothrd
  289.        IF ((NT1 .LE. 1) .OR. (NBTHRS .EQ. 1) .OR. (ITH .GT. 0)) THEN
  290.          NBTHR = 1
  291.          ITHRD = 0
  292.        ELSE
  293.          ithrd=1
  294.          NBTHR = MIN(NT1,NBTHRS)
  295.          call threadii
  296.        ENDIF
  297.  
  298.         SEGINI,MTR2
  299.  
  300. C       Remplissage du 'COMMON/iplnuc' apres THREADII : pthread_mutex_lock
  301. C       sinon soucis de cohabitation entre les ASSISTANTS qui ecrivent tous dans le meme COMMON...
  302.         XP1   = XP
  303.         NBTH1 = NBTHR
  304.         INUA1 = INUAG
  305.         ITR1  = MTR1
  306.         IMAX1 = MAXI1
  307.         IVAL1 = IVAL
  308.         N1    = N
  309.         NN1   = NN
  310.         MM1   = MM
  311.         NNMM1 = NNMM
  312.         IMCHA1= MCHAML
  313.         IMCHA2= MCHAM1
  314.         N3EL  = N1EL
  315.         N3PTEL= N1PTEL
  316.         ITR2  = MTR2
  317.         XELIM1= XELIM
  318.  
  319. C       Lancement des Threads
  320.         if ((nbthr .gt. 1) .AND. (ithrd .eq. 1)) then
  321.           do ith=2,nbthr
  322.             call threadid(ith,ipln2i)
  323.           enddo
  324.           call ipln2i(1)
  325.           do ith=2,nbthr
  326.             call threadif(ith)
  327.           enddo
  328.         else
  329.           call ipln2i(1)
  330.         endif
  331.  
  332.        if (ithrd .eq. 1) call threadis
  333.  
  334.  
  335. C        SEGINI ,MTR2
  336. C        ITR2=MTR2
  337.  
  338. C
  339. C    boucle sur chaque point de gauss
  340. C
  341.  
  342. C        DO 170 J=1,N1EL
  343. C        DO 160 K=1,N1PTEL
  344. C
  345. C        boucle sur les composantes pour remplir le tableau xi
  346. C        qui contient les valeurs des variables maitres
  347. C           DO 140 L=1,NN
  348. C             MELVAL = IELVAL(L)
  349. C             XI((J-1)*N1EL+N1PTEL,L) = VELCHE(K,J)
  350. C  140      CONTINUE
  351. C
  352. C        on interpole le nuage pour les valeurs xi
  353. C
  354. C           CALL IPLNU2(INUA,ITR1,ITR2,MAXI1,IVAL)
  355. C           IF ( IERR .NE. 0 ) THEN
  356. C               SEGDES MCHAML,MCHELM
  357. C               DO 145 L=1,MCHAM1.IELVAL(/1)
  358. C                   MELVAL = MCHAM1.IELVAL(L)
  359. C                   SEGSUP MELVAL
  360. C  145          CONTINUE
  361. C               SEGSUP MCHAM1,MCHEL1,MTR1
  362. C               SEGSUP,MTR2
  363. C               if (ithrd.eq.1) call threadis
  364. C               RETURN
  365. C           ENDIF
  366. C
  367. C        remplissage des resultats stockes dans yi
  368. C
  369. C           DO 150 L=1,MM
  370. C             MELVAL = MCHAM1.IELVAL(L)
  371. C             VELCHE(K,J) = YI((J-1)*N1EL+N1PTEL,L)
  372. C  150      CONTINUE
  373. C           SEGSUP ,MTR2
  374. C
  375. C  160    CONTINUE
  376. C  170    CONTINUE
  377. C
  378.  
  379.          SEGSUP ,MTR1,MTR2
  380.   200 CONTINUE
  381. C     Fin de la boucle sur les MCHAML
  382.  
  383.  
  384.        ICH1 = MCHEL1
  385.        SEGSUP MAXI1
  386.  
  387.        CALL ACTOBJ('MCHAML  ',ICH1,1)
  388.        CALL ECROBJ('MCHAML  ',ICH1)
  389.  
  390.        RETURN
  391.  
  392.  
  393.  
  394.   210 CONTINUE
  395. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  396. C                           Cas d'un CHPOINT                           C
  397. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  398.  
  399.       MCHPOI = IPOI
  400.  
  401.       SEGINI ,MCHPO1 = MCHPOI
  402.       NSOUPO = IPCHP(/1)
  403. C
  404. C boucle sur les CHPOINTS elementaire
  405.       DO 280 I=1,NSOUPO
  406.         MSOUPO = IPCHP(I)
  407. C
  408. C  calcul de nn et mm
  409. C
  410.          NNMM=0
  411.          SEGINI,MTR1
  412.          ITR1 = MTR1
  413.          NN = NOCOMP(/2)
  414.          MM = NUANOM(/2) - NN
  415. C   initialisation du nouveau soupo
  416.          NC = MM
  417.          SEGINI MSOUP1
  418.          MCHPO1.IPCHP(I)=MSOUP1
  419.          MSOUP1.IGEOC = IGEOC
  420. C
  421. C  on verifie que les noms des composantes du CHPOINT sont
  422. C        dans le nuage et on rempli iadd
  423. C
  424.         DO 230 J=1,NOCOMP(/2)
  425.           DO 220 K=1,NUANOM(/2)
  426.             IF (NUANOM(K) .EQ. NOCOMP(J)) THEN
  427.                IADD(**)=K
  428.                NNMM = IADD(/1)
  429.              ENDIF
  430.  220      CONTINUE
  431.           IF ( NNMM .NE. J ) THEN
  432. C        une des composantes du champ n'est pas un composante du nuage
  433.            CALL ERREUR(682)
  434. C          'incoherence entre les composantes du nuage et du champ'
  435.            SEGSUP MTR1,MCHPO1,MSOUP1
  436.            RETURN
  437.           ENDIF
  438.   230     CONTINUE
  439. C
  440. C  on recupere les composantes du nuage qui ne sont pas dans le CHPOINT
  441. C
  442.         DO 250 J=1,NUANOM(/2)
  443.           IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
  444. C            la composante du nuanom(j) du nuage n'estpas scalaire
  445.            MOTERR(1:8) = NUANOM(J)
  446.            MOTERR(9:13) = 'NUAGE'
  447.            CALL ERREUR(681)
  448.            SEGSUP MTR1,MCHPO1,MSOUP1
  449.            RETURN
  450.          ENDIF
  451.           FLAG = .TRUE.
  452.           DO 240 K=1,NOCOMP(/2)
  453.             IF (NUANOM(J) .EQ. NOCOMP(K)) THEN
  454.                FLAG = .FALSE.
  455.             ENDIF
  456.   240     CONTINUE
  457.           IF (FLAG) THEN
  458. C          la composante du nuage n'est pas dans le CHPOINT
  459.             IADD(**) = J
  460.             NNMM = NNMM +1
  461.            ENDIF
  462.   250   CONTINUE
  463.         IF ( NNMM .NE. NUANOM(/2)) THEN
  464.           MOTERR='NUAGE'
  465.           CALL ERREUR(980)
  466.           RETURN
  467.         ENDIF
  468. C
  469. C   on rempli le noms des composantes du nouveau champ
  470. C
  471.         DO 255 J=1,MM
  472.           MSOUP1.NOCOMP(J)=NUANOM(IADD(J+NN))
  473.   255   CONTINUE
  474.         MPOVAL = IPOVAL
  475.  
  476. C       Creation du nouveau segment MPOVAL pour ecrire les resultats
  477.         N = VPOCHA(/1)
  478.         NC = MM
  479.         SEGINI,MPOVA1
  480.         MSOUP1.IPOVAL = MPOVA1
  481.  
  482.  
  483. C       Preparation pour le calcul en parallele
  484.         NT1= N / (100 * nbthrs)
  485. C
  486. C Pour la paralelisation de l'interpolation
  487. C
  488.        ITH = 0
  489.        IF (NBESC .NE. 0) ith=oothrd
  490.        IF ((NT1 .LE. 1) .OR. (NBTHRS .EQ. 1) .OR. (ITH .GT. 0)) THEN
  491.          NBTHR = 1
  492.          ITHRD = 0
  493.        ELSE
  494.          ithrd=1
  495.          NBTHR = MIN(NT1,NBTHRS)
  496.          call threadii
  497.        ENDIF
  498.  
  499. C       Remplissage du 'COMMON/iplnuc' apres THREADII : pthread_mutex_lock
  500. C       sinon soucis de cohabitation entre les ASSISTANTS qui ecrivent tous dans le meme COMMON...
  501.         XP1   = XP
  502.         NBTH1 = NBTHR
  503.         INUA1 = INUAG
  504.         ITR1  = MTR1
  505.         IMAX1 = MAXI1
  506.         IVAL1 = IVAL
  507.         N1    = N
  508.         NN1   = NN
  509.         MM1   = MM
  510.         NNMM1 = NNMM        
  511.         IMPOV1= MPOVAL
  512.         IMPOV2= MPOVA1
  513.         XELIM1= XELIM
  514.  
  515. C       Lancement des Threads
  516.         if ((nbthr.gt.1) .AND. (ithrd.eq.1)) then
  517.           do ith=2,nbthr
  518.             call threadid(ith,ipln2i)
  519.           enddo
  520.           call ipln2i(1)
  521.           do ith=2,nbthr
  522.             call threadif(ith)
  523.           enddo
  524.         else
  525.           call ipln2i(1)
  526.         endif
  527.  
  528.        if (ithrd .eq. 1) call threadis
  529.  
  530.  
  531. C        DO 270 J=1,N
  532. C         Boucle sur les POINTS dont les valeurs sont a evaluer
  533.  
  534. C          SEGINI,MTR2
  535. C          ITR2 = MTR2
  536. C          DO 260 K=1,NN
  537. C           Boucle sur les composantes CONNUES du CHPOINT donne en argument
  538. C            XI(K)=VPOCHA(J,K)
  539. C  260     CONTINUE
  540. C
  541. C   Interpolation
  542. C
  543. C          CALL IPLNU2(INUA,ITR1,ITR2,MAXI1,IVAL)
  544. CC         en cas d'erreur on sort proprement
  545. C          IF ( IERR .NE. 0 ) THEN
  546. C               SEGDES MCHPOI,MSOUPO
  547. C               SEGSUP MCHPO1,MSOUP1,MPOVAL,MTR1,MTR2,MAXI1
  548. C               if (ithrd.eq.1) call threadis
  549. C               RETURN
  550. C          ENDIF
  551. C
  552. C          DO 265 K=1,MM
  553. C           Boucle sur les composantes MANQUANTES du CHPOINT resultat
  554. C            MPOVA1.VPOCHA(J,K)=YI(K)
  555. C  265     CONTINUE
  556. C          SEGSUP,MTR2
  557. C
  558. C  270     CONTINUE
  559.  
  560. C   Desactivation des SEGMENTS
  561.           SEGSUP,MTR1
  562. C
  563.   280 CONTINUE
  564. C
  565.       IP1 = MCHPO1
  566.       SEGSUP ,MAXI1
  567.  
  568.       CALL ACTOBJ('CHPOINT ',IP1,1)
  569.       CALL ECROBJ('CHPOINT ',IP1)
  570.  
  571.       END
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  

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