Télécharger la plaquette Cast3M


fr en

Télécharger fuchpo.eso

Retour à la liste

Numérotation des lignes : 

fuchpo

  1. C FUCHPO    SOURCE    SP204843  24/10/09    21:15:05     12027          
  2.       SUBROUTINE FUCHPO(IP1,IP2,IRET)
  3. C======================================================================
  4. C  fonction:
  5. C  sous routine pour fusionner deux champs par points diffus
  6. C
  7. C  arguments:
  8. C  ip1 (E) pointeur sur le premier des deux champ par point
  9. C  ip2 (E) pointeur sur le second des deux champ par point
  10. C  iret (S) pointeur sur le champ par point resultat
  11. C
  12. C  variables:
  13. C
  14. C * mtrav : - bb(i,j) est la valeur de la ieme inconnue de champ pour
  15. C             le jieme noeud du tableau igeo .
  16. C           - inco(nnin) contient le nom des nnin inconnues differentes
  17. C           - ibin(i,j)=1 ou 0 indique que la ieme inconnue du champ
  18. C             existe pour le jieme noeud du tableau igeo .
  19. C           - igeo(i) est le numero a mettre dans un objet meleme pour
  20. C          referencer le ieme noeud .
  21. C
  22. C * mtra  : - nopoin(i) adresse de colonne dans bb et ibin des valeurs
  23. C          correspondant au noeud i .
  24. C
  25. C * mtr1  : - ipcom liste des noms des inconnues permet de creer inco .
  26. C
  27. C * mtr2  : - iicom adresse dans ipcom des inconnues correspondant au
  28. C          2ieme ch point .
  29. C
  30. C * mtr3  : - index tableau de correspondance entre les supports geome-
  31. C          triques du 1er chpoint et du 2ieme chpoint .
  32. C
  33. C  auteur: A de Gayffier 13/06/94
  34. C======================================================================
  35.       IMPLICIT INTEGER(I-N)
  36.       IMPLICIT REAL*8(A-H,O-Z)
  37. -INC SMCHPOI
  38.  
  39. -INC PPARAM
  40. -INC CCOPTIO
  41. -INC SMELEME
  42. -INC SMCOORD
  43. -INC TMTRAV
  44.       SEGMENT/MTRA/(NOPOIN(nbpts))
  45.       SEGMENT MTR1
  46.        CHARACTER*(LOCOMP) IPCOM(0)
  47.       ENDSEGMENT
  48.       SEGMENT/MTR2/(IICOM(0))
  49.       SEGMENT/MTR3/(INDEX(0))
  50.       SEGMENT/MTR4/(IPHAR(0))
  51. C     ordre de grandeur des composantes
  52.       SEGMENT/MTR5/(DMOY(NNIN))
  53. C
  54.       DIMENSION IPO(2)
  55.       CHARACTER*8 MOT
  56.       DATA UN,ZERO/1.D0,0.D0/
  57.       character*4 mcle(1)
  58.       data mcle/'NOER'/
  59. C
  60.       IRET = 0
  61.  
  62.       noer=0
  63.       call lirmot(mcle,1,noer,0)
  64.       if (ierr.ne.0) return
  65.  
  66.       MCHPO1 = IP1
  67.       MCHPO2 = IP2
  68.       SEGACT,MCHPO1,MCHPO2
  69.  
  70.       NSOUP1 = MCHPO1.IPCHP(/1)
  71.       NSOUP2 = MCHPO2.IPCHP(/1)
  72.  
  73.       NAT1 = MCHPO1.JATTRI(1)
  74.       NAT2 = MCHPO2.JATTRI(1)
  75.  
  76. *     Si CHPOINT vide, on renvoie l'autre si il est non vide:
  77.       IF (NSOUP1.EQ.0) THEN
  78.         IRET = MCHPO2
  79.         RETURN
  80.       ENDIF
  81.       IF (NSOUP2.EQ.0) then
  82.         IRET = MCHPO1
  83.         RETURN
  84.       ENDIF
  85. C
  86. C    verification de la compatibilité des natures
  87. C
  88.       IF ( (NAT1*NAT2) .EQ. 0) THEN
  89. C    une des deux natures est indeterminée
  90.          CALL ERREUR(650)
  91.          RETURN
  92.       ELSE
  93.         IF ((NAT1 .EQ. 2) .AND. (NAT2 .EQ. 2)) THEN
  94. C      les deux champ sont discrets: on somme les composantes communes
  95.            CALL ADCHPO(IP1,IP2,IRET,1D0,1D0)
  96.            RETURN
  97.         ENDIF
  98.         IF ((NAT1 .NE. 1) .OR. (NAT2 .NE. 1)) THEN
  99. C      les natures ne sont pas compatibles
  100.            CALL ERREUR(649)
  101.            RETURN
  102.         ENDIF
  103.       ENDIF
  104. C
  105. C Petite verification sur les modes de calcul
  106.       ifo1 = MCHPO1.IFOPOI
  107.       ifo2 = MCHPO2.IFOPOI
  108.       ifos = ifo1
  109.       IF (ifo1 .NE. ifo2) THEN
  110.         interr(1)=ifo1
  111.         interr(2)=ifo2
  112.         interr(3)=IFOUR
  113. c-dbg      write(ioimp,*) '1132 FUCHPO',ip1,ip2
  114.         call erreur(1132)
  115.         ifos = IFOUR
  116.       END IF
  117.  
  118. C   les deux champs sont de nature diffuse
  119. C   on moyenne les composantes communes
  120. C
  121.       IF ( IP1 .NE. IP2) GOTO 60
  122. C
  123. C *** cas ou les 2 pointeurs ip1 et ip2 sont egaux
  124. C
  125. c*      SEGACT MCHPO1
  126.       NSOUPO=NSOUP1
  127.       NAT   =NAT1
  128.       SEGINI MCHPOI
  129.       DO 10 I=1,NAT
  130.         JATTRI(I)=MCHPO1.JATTRI(I)
  131.  10   CONTINUE
  132.       MOCHDE=MCHPO1.MOCHDE
  133.       MTYPOI=MCHPO1.MTYPOI
  134.       IFOPOI=ifos
  135.       DO 50 IA=1,NSOUPO
  136.         MSOUP1=MCHPO1.IPCHP(IA)
  137.         SEGACT MSOUP1
  138.         NC=MSOUP1.NOCOMP(/2)
  139.         SEGINI MSOUPO
  140.         IPCHP(IA)=MSOUPO
  141.         DO 20 IB=1,NC
  142.           NOCOMP(IB)=MSOUP1.NOCOMP(IB)
  143.           NOHARM(IB)=MSOUP1.NOHARM(IB)
  144.   20    CONTINUE
  145.         IGEOC=MSOUP1.IGEOC
  146.         MPOVA1=MSOUP1.IPOVAL
  147.         SEGACT MPOVA1
  148.         N  =MPOVA1.VPOCHA(/1)
  149.         NC1=MPOVA1.VPOCHA(/2)
  150. C
  151. C     erreur pb dimension tableau voir routine adchpo
  152.         IF (NC1.NE.NC) THEN
  153.           IRET=0
  154.           SEGSUP MSOUPO,MCHPOI
  155.           CALL ERREUR(114)
  156.           RETURN
  157.         ENDIF
  158. C
  159.         SEGINI MPOVAL
  160.         IPOVAL=MPOVAL
  161.         DO 40 IC=1,NC
  162.           DO 41 IB=1,N
  163.             VPOCHA(IB,IC)=MPOVA1.VPOCHA(IB,IC)
  164.    41     CONTINUE
  165.    40   CONTINUE
  166.    50 CONTINUE
  167. C
  168. C      on sort de la sous routine
  169.       IRET=MCHPOI
  170.       GOTO 666
  171. C
  172. C *** cas  ou les pointeurs ip1 et ip2 sont differents
  173. C
  174.   60  CONTINUE
  175.       IPO(1)=IP1
  176.       IPO(2)=IP2
  177.       MOT=MCHPO1.MTYPOI
  178.       IF(MOT.NE.MCHPO2.MTYPOI) THEN
  179.         MOT='ET OU +'
  180.       ENDIF
  181. C
  182. C     on verifie que les nbres de sous paquets sont egaux
  183. C
  184.       IF(NSOUP1.EQ.NSOUP2) GO TO 75
  185. C      traitement par la methode générale
  186.       GO TO 300
  187. C
  188. C     on regarde si les supports geometriques sont identiques a une
  189. C     permutation pres
  190. C
  191.    75 SEGINI MTR3
  192.       DO 100 IA=1,NSOUP1
  193.         MSOUP1=MCHPO1.IPCHP(IA)
  194.         SEGACT MSOUP1
  195.         DO 80 IB=1,NSOUP2
  196.           MSOUP2=MCHPO2.IPCHP(IB)
  197.           SEGACT MSOUP2
  198.           IF(MSOUP1.IGEOC.EQ.MSOUP2.IGEOC) GO TO 90
  199.    80   CONTINUE
  200. C
  201. C     il n y a pas egalite des supports geometriques a une permutation
  202. C     pres
  203. C
  204.         SEGSUP MTR3
  205. C      traitement par la methode générale
  206.         GO TO 300
  207. C
  208.    90   CONTINUE
  209. C    la permutation est rangée dans  index
  210.         INDEX(**)=IB
  211.   100 CONTINUE
  212. C
  213. C *** cas ou il y a egalite des supports geometriques a une permutation
  214. C     pres
  215. C
  216.       NSOUPO=NSOUP1
  217.       NAT=MAX(NAT1,NAT2,1)
  218.       SEGINI MCHPOI
  219.       JATTRI(1) = 1
  220.       IRET=MCHPOI
  221.       MTYPOI=MOT
  222.       MOCHDE=MCHPO1.MOCHDE
  223.       IFOPOI=ifos
  224. C
  225.       DO 250 IA=1,NSOUP1
  226.         MSOUP1=MCHPO1.IPCHP(IA)
  227.         MSOUP2=MCHPO2.IPCHP(INDEX(IA))
  228.         SEGACT MSOUP1,MSOUP2
  229. C
  230. C     comparaison des noms des composantes
  231. C
  232.         SEGINI MTR1,MTR4
  233.         NC1=MSOUP1.NOCOMP(/2)
  234.         NC2=MSOUP2.NOCOMP(/2)
  235.         DO 130 IB=1,NC1
  236.          IPCOM(**)=MSOUP1.NOCOMP(IB)
  237.          IPHAR(**)=MSOUP1.NOHARM(IB)
  238.   130   CONTINUE
  239.         SEGINI MTR2
  240.         DO 160 IB=1,NC2
  241.           DO 140 IC=1,NC1
  242.            IF(MSOUP2.NOCOMP(IB).NE.MSOUP1.NOCOMP(IC)) GOTO 140
  243.            IF(MSOUP2.NOHARM(IB).EQ.MSOUP1.NOHARM(IC)) GOTO 150
  244.   140    CONTINUE
  245. C         la composante du IB n'est pas commune
  246.          IPCOM(**)=MSOUP2.NOCOMP(IB)
  247.          IPHAR(**)=MSOUP2.NOHARM(IB)
  248.          IICOM(**)=IPCOM(/2)
  249.          GO TO 160
  250.   150    CONTINUE
  251. C          la composante est commune
  252.          IICOM(**)=IC
  253.   160   CONTINUE
  254. C
  255.         MPOVA1=MSOUP1.IPOVAL
  256.         MPOVA2=MSOUP2.IPOVAL
  257.         SEGACT MPOVA1,MPOVA2
  258.         N1=MPOVA1.VPOCHA(/1)
  259.         N2=MPOVA2.VPOCHA(/1)
  260.         NCX1=MPOVA1.VPOCHA(/2)
  261.         NCX2=MPOVA2.VPOCHA(/2)
  262.         IF(NCX1.NE.NC1) GOTO 170
  263.         IF(NCX2.NE.NC2) GOTO 170
  264.         IF(N1.NE.N2)    GOTO 170
  265.         GOTO 180
  266.   170   CONTINUE
  267.         SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  268. C
  269. C     pb avec les dimensions des chpoints
  270. C
  271.         CALL ERREUR(114)
  272.         RETURN
  273.         IRET=0
  274. C      on sort de la sous routine
  275.         GOTO 666
  276. C
  277.   180   CONTINUE
  278.         N=N1
  279.         NC=IPCOM(/2)
  280.         SEGINI MPOVAL
  281.         NNIN = NC
  282.         SEGINI MTR5
  283. C
  284. C     mise a 0 de vpocha
  285. C
  286. *        DO 190 IB=1,N
  287. *        DO 190 IC=1,NC
  288. *         VPOCHA(IB,IC)=ZERO
  289. *  190   CONTINUE
  290. C
  291. C     addition des chpoints
  292. C
  293. C     on place les valeurs de MCHPO1
  294.         DO 210 IC=1,NC1
  295.           DO 200 IB=1,N
  296.             VPOCHA(IB,IC) = MPOVA1.VPOCHA(IB,IC)+VPOCHA(IB,IC)
  297.             DMOY(IC) = DMOY(IC) + ABS(VPOCHA(IB,IC)/N)
  298.   200     CONTINUE
  299.           IF (IIMPI.EQ.123)
  300.      &      write (IOIMP,*) ' ic dmoy(ic) ',ic,dmoy(ic)
  301.   210   CONTINUE
  302. C
  303.         DO 230 IC=1,NC2
  304.           IIC=IICOM(IC)
  305.           DO 220 IB=1,N
  306.             IF (IIC .LE. NC1 ) THEN
  307. C      il s'agit d'ne composante commune
  308.               XX1 = MPOVA2.VPOCHA(IB,IC)
  309.               XX2 = VPOCHA(IB,IIC)
  310.               DXX = ABS ( XX2 - XX1)
  311. *             SXX = MIN(ABS ( XX1 + XX2 ) / 2.D0,1.D-50)
  312.               SXX = DMOY(IIC)
  313.               IF (DXX .LE. (1.D-4*SXX) .or.noer.eq.1) THEN
  314.                 VPOCHA(IB,IIC)= ( XX1 + XX2 ) / 2.D0
  315.               ELSE
  316. C        les valeurs des champ diffus au meme point sont différentes
  317.                IF (IIMPI.EQ.123)
  318.      &           write (IOIMP,*) xx1,xx2,SXX,DXX
  319.                  CALL ERREUR(651)
  320.                  RETURN
  321.                SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  322. C         on sort
  323.                GOTO 666
  324.               ENDIF
  325.             ELSE
  326. C    composantes non communes
  327.               VPOCHA(IB,IIC) = MPOVA2.VPOCHA(IB,IC)+VPOCHA(IB,IIC)
  328.             ENDIF
  329.   220     CONTINUE
  330.   230   CONTINUE
  331. C
  332.         SEGINI MSOUPO
  333.         DO 240 IB=1,NC
  334.          NOCOMP(IB)=IPCOM(IB)
  335.          NOHARM(IB)=IPHAR(IB)
  336.   240   CONTINUE
  337.         SEGSUP MTR1,MTR2,MTR4
  338.         IPOVAL=MPOVAL
  339.         IPT2=MSOUP1.IGEOC
  340. ****    SEGINI,IPT1=IPT2
  341.         IPT1 = IPT2
  342.         IGEOC=IPT1
  343.         IPCHP(IA)=MSOUPO
  344.         SEGSUP MTR5
  345.   250 CONTINUE
  346. C
  347.       SEGSUP MTR3
  348. C     on sort
  349.       GOTO 666
  350. C
  351. C *** cas ou les supports geometriques ne sont pas identiques
  352. C     a une permutation pres
  353. C
  354.   300 CONTINUE
  355. C
  356. C ****  a-t-on affaires a des champoints vides?
  357. C
  358.       MCHPOI=IP1
  359. c*      SEGACT MCHPOI
  360.       NS1=IPCHP(/1)
  361.       MCHPO2=IP2
  362. c*      SEGACT MCHPO2
  363.       NS2=MCHPO2.IPCHP(/1)
  364.       IF(NS1*NS2.NE.0) GO TO 3001
  365.       IF(NS1+NS2.NE.0) THEN
  366. C un seul des chpoints est vide
  367.           IF(NS1.EQ.0)  IP1=IP2
  368.           CALL ECRCHA('GEOM')
  369.           CALL ECROBJ('CHPOINT ',IP1)
  370.           CALL COPIER
  371.           CALL LIROBJ('CHPOINT',IRET,1,IRETOU)
  372.       ELSE
  373. C  les deux chpoints sont vides
  374.          NSOUPO=0
  375.          NAT=1
  376.          SEGINI MCHPOI
  377.          IFOPOI = ifos
  378.          IRET = MCHPOI
  379.       ENDIF
  380.       RETURN
  381. C
  382.  3001 CONTINUE
  383.       SEGINI MTRA,MTR1,MTR4
  384. C
  385. C     mise a zero de nopoin
  386. C
  387. *      DO 305 IA=1,nbpts
  388. *      NOPOIN(IA)=0
  389. *  305 CONTINUE
  390. C
  391.       MCHPOI=IP1
  392. c*      SEGACT MCHPOI
  393.       MSOUPO=IPCHP(1)
  394.       SEGACT MSOUPO
  395.       IPCOM(**)=NOCOMP(1)
  396.       IPHAR(**)=NOHARM(1)
  397.       NC=1
  398.       IK=0
  399. C
  400. C     creation de nopoin et de ipcom
  401. C
  402.       DO 360 ICH=1,2
  403.         MCHPOI=IPO(ICH)
  404. c*        SEGACT MCHPOI
  405.         NSOUPO=IPCHP(/1)
  406. C
  407. C     boucle sur les sous references d 1 chpoint
  408. C
  409.         DO 350 IA=1,NSOUPO
  410.           MSOUPO=IPCHP(IA)
  411.           SEGACT MSOUPO
  412.           MELEME=IGEOC
  413.           SEGACT MELEME
  414.           NBNN  =NUM(/1)
  415.           NBELEM=NUM(/2)
  416.           DO 310 IB=1,NBELEM
  417.             K=NUM(1,IB)
  418.             IF(NOPOIN(K).NE.0) GO TO 310
  419.             IK=IK+1
  420.             NOPOIN(K)=IK
  421.   310     CONTINUE
  422.           NCBBB=NOCOMP(/2)
  423.           DO 330 IB=1,NCBBB
  424.             NC=IPCOM(/2)
  425.             DO 320 IC=1,NC
  426.               IF(IPCOM(IC).NE.NOCOMP(IB)) GO TO 320
  427.               IF(IPHAR(IC).EQ.NOHARM(IB))  GO TO 330
  428.   320       CONTINUE
  429.             IPCOM(**)=NOCOMP(IB)
  430.             IPHAR(**)=NOHARM(IB)
  431.             NC=NC+1
  432.   330     CONTINUE
  433.   350   CONTINUE
  434.   360 CONTINUE
  435. C
  436.       NNIN=NC
  437.       NNNOE=IK
  438.       SEGINI MTRAV
  439. C
  440. C     initialisation a zero des tableaux
  441. C
  442.       SEGINI MTR5
  443. C
  444. *      DO 370 IB=1,NNNOE
  445. *      DO 370 IA=1,NNIN
  446. *      BB(IA,IB)=ZERO
  447. *      IBIN(IA,IB)=0
  448. *      IMOY(IA,IB) = 0
  449. *  370 CONTINUE
  450. C
  451. C     creation de inco
  452. C
  453.       DO 380 IA=1,NNIN
  454.         INCO(IA)=IPCOM(IA)
  455.         NHAR(IA)=IPHAR(IA)
  456.   380 CONTINUE
  457. C
  458. C     creation de bb,ibin,igeo
  459. C
  460.       DO 450 ICH=1,2
  461.         MCHPOI=IPO(ICH)
  462. c*        SEGACT MCHPOI
  463.         NSOUPO=IPCHP(/1)
  464.         DO 430 IA=1,NSOUPO
  465.           MSOUPO=IPCHP(IA)
  466.           SEGACT MSOUPO
  467.           MELEME=IGEOC
  468.           SEGACT MELEME
  469.           MPOVAL=IPOVAL
  470.           SEGACT MPOVAL
  471.           NBELEM=NUM(/2)
  472.           N=VPOCHA(/1)
  473.           NC=VPOCHA(/2)
  474.           NC1=NOCOMP(/2)
  475. C
  476.          DO 420 IB=1,NC1
  477.            DO 390 IC=1,NNIN
  478.              IF(NOCOMP(IB).NE.IPCOM(IC)) GO TO 390
  479.              IF(NOHARM(IB).EQ.IPHAR(IC)) GO TO 400
  480.   390      CONTINUE
  481.   400      CONTINUE
  482.            DO 411 ID=1,NBELEM
  483.              DMOY(IB)=DMOY(IB)+ABS(VPOCHA(ID,IB)/NBELEM)
  484.   411      CONTINUE
  485.            DO 410 ID=1,NBELEM
  486.               KI=NOPOIN(NUM(1,ID))
  487.               IGEO(KI)=NUM(1,ID)
  488.               IF (  IBIN(IC,KI) .EQ. 1) THEN
  489. C          la valeur au point est defini dans les deux champs
  490.                XX1 = BB(IC,KI)
  491.                XX2 = VPOCHA(ID,IB)
  492.                DXX = ABS ( XX2 - XX1 )
  493.                SXX = DMOY(IB)
  494.                IF ( DXX .LE. (1.D-4*SXX).or.noer.eq.1) THEN
  495.                   BB(IC,KI) = ( XX1 + XX2 ) / 2.D0
  496.                ELSE
  497. C          les valeurs des champs au meme point sont différentes
  498.                 IF (IIMPI.EQ.123)
  499.      &            write (IOIMP,*) xx1,xx2,sxx,DXX
  500.                 CALL ERREUR (651)
  501.                 RETURN
  502.                 SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  503.                 GOTO 666
  504.                ENDIF
  505.              ELSE
  506.                BB(IC,KI)=BB(IC,KI)+VPOCHA(ID,IB)
  507. *              DMOY(IC) = DMOY(IC) +ABS(BB(IC,KI)/NNNOE)
  508.              ENDIF
  509.              IBIN(IC,KI)=1
  510.   410      CONTINUE
  511.   420    CONTINUE
  512.   430  CONTINUE
  513.   450 CONTINUE
  514.       ITRAV=MTRAV
  515. C
  516. C     reconstuction de la partition
  517. C
  518.       CALL CRECHP(ITRAV,ICHPOI)
  519. C
  520.       SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  521.       IRET=ICHPOI
  522.       MCHPOI=ICHPOI
  523. c*      MCHPO1 = IP1
  524. c*      MCHPO2 = IP2
  525. c*      NAT1 = MCHPO1.JATTRI(/1)
  526. c*      NAT2 = MCHPO2.JATTRI(/1)
  527.       NAT=MAX(NAT1,NAT2,1)
  528.       NSOUPO = IPCHP(/1)
  529.       SEGADJ MCHPOI
  530.       IRET=MCHPOI
  531.       MTYPOI=MOT
  532.       IFOPOI = ifos
  533.       JATTRI(1) = 1
  534. C
  535.  666  CONTINUE
  536. c*      RETURN
  537.       END
  538.  
  539.  
  540.  
  541.  
  542.  

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