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    GOUNAND   25/05/05    21:15:04     12259          
  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.         moterr(1:8)='CHPOINT'
  111.         interr(1)=ifo1
  112.         interr(2)=ifo2
  113.         interr(3)=IFOUR
  114. c-dbg      write(ioimp,*) '1132 FUCHPO',ip1,ip2
  115.         call erreur(1132)
  116.         ifos = IFOUR
  117.       END IF
  118.  
  119. C   les deux champs sont de nature diffuse
  120. C   on moyenne les composantes communes
  121. C
  122.       IF ( IP1 .NE. IP2) GOTO 60
  123. C
  124. C *** cas ou les 2 pointeurs ip1 et ip2 sont egaux
  125. C
  126. c*      SEGACT MCHPO1
  127.       NSOUPO=NSOUP1
  128.       NAT   =NAT1
  129.       SEGINI MCHPOI
  130.       DO 10 I=1,NAT
  131.         JATTRI(I)=MCHPO1.JATTRI(I)
  132.  10   CONTINUE
  133.       MOCHDE=MCHPO1.MOCHDE
  134.       MTYPOI=MCHPO1.MTYPOI
  135.       IFOPOI=ifos
  136.       DO 50 IA=1,NSOUPO
  137.         MSOUP1=MCHPO1.IPCHP(IA)
  138.         SEGACT MSOUP1
  139.         NC=MSOUP1.NOCOMP(/2)
  140.         SEGINI MSOUPO
  141.         IPCHP(IA)=MSOUPO
  142.         DO 20 IB=1,NC
  143.           NOCOMP(IB)=MSOUP1.NOCOMP(IB)
  144.           NOHARM(IB)=MSOUP1.NOHARM(IB)
  145.   20    CONTINUE
  146.         IGEOC=MSOUP1.IGEOC
  147.         MPOVA1=MSOUP1.IPOVAL
  148.         SEGACT MPOVA1
  149.         N  =MPOVA1.VPOCHA(/1)
  150.         NC1=MPOVA1.VPOCHA(/2)
  151. C
  152. C     erreur pb dimension tableau voir routine adchpo
  153.         IF (NC1.NE.NC) THEN
  154.           IRET=0
  155.           SEGSUP MSOUPO,MCHPOI
  156.           CALL ERREUR(114)
  157.           RETURN
  158.         ENDIF
  159. C
  160.         SEGINI MPOVAL
  161.         IPOVAL=MPOVAL
  162.         DO 40 IC=1,NC
  163.           DO 41 IB=1,N
  164.             VPOCHA(IB,IC)=MPOVA1.VPOCHA(IB,IC)
  165.    41     CONTINUE
  166.    40   CONTINUE
  167.    50 CONTINUE
  168. C
  169. C      on sort de la sous routine
  170.       IRET=MCHPOI
  171.       GOTO 666
  172. C
  173. C *** cas  ou les pointeurs ip1 et ip2 sont differents
  174. C
  175.   60  CONTINUE
  176.       IPO(1)=IP1
  177.       IPO(2)=IP2
  178.       MOT=MCHPO1.MTYPOI
  179.       IF(MOT.NE.MCHPO2.MTYPOI) THEN
  180.         MOT='ET OU +'
  181.       ENDIF
  182. C
  183. C     on verifie que les nbres de sous paquets sont egaux
  184. C
  185.       IF(NSOUP1.EQ.NSOUP2) GO TO 75
  186. C      traitement par la methode générale
  187.       GO TO 300
  188. C
  189. C     on regarde si les supports geometriques sont identiques a une
  190. C     permutation pres
  191. C
  192.    75 SEGINI MTR3
  193.       DO 100 IA=1,NSOUP1
  194.         MSOUP1=MCHPO1.IPCHP(IA)
  195.         SEGACT MSOUP1
  196.         DO 80 IB=1,NSOUP2
  197.           MSOUP2=MCHPO2.IPCHP(IB)
  198.           SEGACT MSOUP2
  199.           IF(MSOUP1.IGEOC.EQ.MSOUP2.IGEOC) GO TO 90
  200.    80   CONTINUE
  201. C
  202. C     il n y a pas egalite des supports geometriques a une permutation
  203. C     pres
  204. C
  205.         SEGSUP MTR3
  206. C      traitement par la methode générale
  207.         GO TO 300
  208. C
  209.    90   CONTINUE
  210. C    la permutation est rangée dans  index
  211.         INDEX(**)=IB
  212.   100 CONTINUE
  213. C
  214. C *** cas ou il y a egalite des supports geometriques a une permutation
  215. C     pres
  216. C
  217.       NSOUPO=NSOUP1
  218.       NAT=MAX(NAT1,NAT2,1)
  219.       SEGINI MCHPOI
  220.       JATTRI(1) = 1
  221.       IRET=MCHPOI
  222.       MTYPOI=MOT
  223.       MOCHDE=MCHPO1.MOCHDE
  224.       IFOPOI=ifos
  225. C
  226.       DO 250 IA=1,NSOUP1
  227.         MSOUP1=MCHPO1.IPCHP(IA)
  228.         MSOUP2=MCHPO2.IPCHP(INDEX(IA))
  229.         SEGACT MSOUP1,MSOUP2
  230. C
  231. C     comparaison des noms des composantes
  232. C
  233.         SEGINI MTR1,MTR4
  234.         NC1=MSOUP1.NOCOMP(/2)
  235.         NC2=MSOUP2.NOCOMP(/2)
  236.         DO 130 IB=1,NC1
  237.          IPCOM(**)=MSOUP1.NOCOMP(IB)
  238.          IPHAR(**)=MSOUP1.NOHARM(IB)
  239.   130   CONTINUE
  240.         SEGINI MTR2
  241.         DO 160 IB=1,NC2
  242.           DO 140 IC=1,NC1
  243.            IF(MSOUP2.NOCOMP(IB).NE.MSOUP1.NOCOMP(IC)) GOTO 140
  244.            IF(MSOUP2.NOHARM(IB).EQ.MSOUP1.NOHARM(IC)) GOTO 150
  245.   140    CONTINUE
  246. C         la composante du IB n'est pas commune
  247.          IPCOM(**)=MSOUP2.NOCOMP(IB)
  248.          IPHAR(**)=MSOUP2.NOHARM(IB)
  249.          IICOM(**)=IPCOM(/2)
  250.          GO TO 160
  251.   150    CONTINUE
  252. C          la composante est commune
  253.          IICOM(**)=IC
  254.   160   CONTINUE
  255. C
  256.         MPOVA1=MSOUP1.IPOVAL
  257.         MPOVA2=MSOUP2.IPOVAL
  258.         SEGACT MPOVA1,MPOVA2
  259.         N1=MPOVA1.VPOCHA(/1)
  260.         N2=MPOVA2.VPOCHA(/1)
  261.         NCX1=MPOVA1.VPOCHA(/2)
  262.         NCX2=MPOVA2.VPOCHA(/2)
  263.         IF(NCX1.NE.NC1) GOTO 170
  264.         IF(NCX2.NE.NC2) GOTO 170
  265.         IF(N1.NE.N2)    GOTO 170
  266.         GOTO 180
  267.   170   CONTINUE
  268.         SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  269. C
  270. C     pb avec les dimensions des chpoints
  271. C
  272.         CALL ERREUR(114)
  273.         RETURN
  274.         IRET=0
  275. C      on sort de la sous routine
  276.         GOTO 666
  277. C
  278.   180   CONTINUE
  279.         N=N1
  280.         NC=IPCOM(/2)
  281.         SEGINI MPOVAL
  282.         NNIN = NC
  283.         SEGINI MTR5
  284. C
  285. C     mise a 0 de vpocha
  286. C
  287. *        DO 190 IB=1,N
  288. *        DO 190 IC=1,NC
  289. *         VPOCHA(IB,IC)=ZERO
  290. *  190   CONTINUE
  291. C
  292. C     addition des chpoints
  293. C
  294. C     on place les valeurs de MCHPO1
  295.         DO 210 IC=1,NC1
  296.           DO 200 IB=1,N
  297.             VPOCHA(IB,IC) = MPOVA1.VPOCHA(IB,IC)+VPOCHA(IB,IC)
  298.             DMOY(IC) = DMOY(IC) + ABS(VPOCHA(IB,IC)/N)
  299.   200     CONTINUE
  300.           IF (IIMPI.EQ.123)
  301.      &      write (IOIMP,*) ' ic dmoy(ic) ',ic,dmoy(ic)
  302.   210   CONTINUE
  303. C
  304.         DO 230 IC=1,NC2
  305.           IIC=IICOM(IC)
  306.           DO 220 IB=1,N
  307.             IF (IIC .LE. NC1 ) THEN
  308. C      il s'agit d'ne composante commune
  309.               XX1 = MPOVA2.VPOCHA(IB,IC)
  310.               XX2 = VPOCHA(IB,IIC)
  311.               DXX = ABS ( XX2 - XX1)
  312. *             SXX = MIN(ABS ( XX1 + XX2 ) / 2.D0,1.D-50)
  313.               SXX = DMOY(IIC)
  314.               IF (DXX .LE. (1.D-4*SXX) .or.noer.eq.1) THEN
  315.                 VPOCHA(IB,IIC)= ( XX1 + XX2 ) / 2.D0
  316.               ELSE
  317. C        les valeurs des champ diffus au meme point sont différentes
  318.                IF (IIMPI.EQ.123)
  319.      &           write (IOIMP,*) xx1,xx2,SXX,DXX
  320.                  CALL ERREUR(651)
  321.                  RETURN
  322.                SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  323. C         on sort
  324.                GOTO 666
  325.               ENDIF
  326.             ELSE
  327. C    composantes non communes
  328.               VPOCHA(IB,IIC) = MPOVA2.VPOCHA(IB,IC)+VPOCHA(IB,IIC)
  329.             ENDIF
  330.   220     CONTINUE
  331.   230   CONTINUE
  332. C
  333.         SEGINI MSOUPO
  334.         DO 240 IB=1,NC
  335.          NOCOMP(IB)=IPCOM(IB)
  336.          NOHARM(IB)=IPHAR(IB)
  337.   240   CONTINUE
  338.         SEGSUP MTR1,MTR2,MTR4
  339.         IPOVAL=MPOVAL
  340.         IPT2=MSOUP1.IGEOC
  341. ****    SEGINI,IPT1=IPT2
  342.         IPT1 = IPT2
  343.         IGEOC=IPT1
  344.         IPCHP(IA)=MSOUPO
  345.         SEGSUP MTR5
  346.   250 CONTINUE
  347. C
  348.       SEGSUP MTR3
  349. C     on sort
  350.       GOTO 666
  351. C
  352. C *** cas ou les supports geometriques ne sont pas identiques
  353. C     a une permutation pres
  354. C
  355.   300 CONTINUE
  356. C
  357. C ****  a-t-on affaires a des champoints vides?
  358. C
  359.       MCHPOI=IP1
  360. c*      SEGACT MCHPOI
  361.       NS1=IPCHP(/1)
  362.       MCHPO2=IP2
  363. c*      SEGACT MCHPO2
  364.       NS2=MCHPO2.IPCHP(/1)
  365.       IF(NS1*NS2.NE.0) GO TO 3001
  366.       IF(NS1+NS2.NE.0) THEN
  367. C un seul des chpoints est vide
  368.           IF(NS1.EQ.0)  IP1=IP2
  369.           CALL ECRCHA('GEOM')
  370.           CALL ECROBJ('CHPOINT ',IP1)
  371.           CALL COPIER
  372.           CALL LIROBJ('CHPOINT',IRET,1,IRETOU)
  373.       ELSE
  374. C  les deux chpoints sont vides
  375.          NSOUPO=0
  376.          NAT=1
  377.          SEGINI MCHPOI
  378.          IFOPOI = ifos
  379.          IRET = MCHPOI
  380.       ENDIF
  381.       RETURN
  382. C
  383.  3001 CONTINUE
  384.       SEGINI MTRA,MTR1,MTR4
  385. C
  386. C     mise a zero de nopoin
  387. C
  388. *      DO 305 IA=1,nbpts
  389. *      NOPOIN(IA)=0
  390. *  305 CONTINUE
  391. C
  392.       MCHPOI=IP1
  393. c*      SEGACT MCHPOI
  394.       MSOUPO=IPCHP(1)
  395.       SEGACT MSOUPO
  396.       IPCOM(**)=NOCOMP(1)
  397.       IPHAR(**)=NOHARM(1)
  398.       NC=1
  399.       IK=0
  400. C
  401. C     creation de nopoin et de ipcom
  402. C
  403.       DO 360 ICH=1,2
  404.         MCHPOI=IPO(ICH)
  405. c*        SEGACT MCHPOI
  406.         NSOUPO=IPCHP(/1)
  407. C
  408. C     boucle sur les sous references d 1 chpoint
  409. C
  410.         DO 350 IA=1,NSOUPO
  411.           MSOUPO=IPCHP(IA)
  412.           SEGACT MSOUPO
  413.           MELEME=IGEOC
  414.           SEGACT MELEME
  415.           NBNN  =NUM(/1)
  416.           NBELEM=NUM(/2)
  417.           DO 310 IB=1,NBELEM
  418.             K=NUM(1,IB)
  419.             IF(NOPOIN(K).NE.0) GO TO 310
  420.             IK=IK+1
  421.             NOPOIN(K)=IK
  422.   310     CONTINUE
  423.           NCBBB=NOCOMP(/2)
  424.           DO 330 IB=1,NCBBB
  425.             NC=IPCOM(/2)
  426.             DO 320 IC=1,NC
  427.               IF(IPCOM(IC).NE.NOCOMP(IB)) GO TO 320
  428.               IF(IPHAR(IC).EQ.NOHARM(IB))  GO TO 330
  429.   320       CONTINUE
  430.             IPCOM(**)=NOCOMP(IB)
  431.             IPHAR(**)=NOHARM(IB)
  432.             NC=NC+1
  433.   330     CONTINUE
  434.   350   CONTINUE
  435.   360 CONTINUE
  436. C
  437.       NNIN=NC
  438.       NNNOE=IK
  439.       SEGINI MTRAV
  440. C
  441. C     initialisation a zero des tableaux
  442. C
  443.       SEGINI MTR5
  444. C
  445. *      DO 370 IB=1,NNNOE
  446. *      DO 370 IA=1,NNIN
  447. *      BB(IA,IB)=ZERO
  448. *      IBIN(IA,IB)=0
  449. *      IMOY(IA,IB) = 0
  450. *  370 CONTINUE
  451. C
  452. C     creation de inco
  453. C
  454.       DO 380 IA=1,NNIN
  455.         INCO(IA)=IPCOM(IA)
  456.         NHAR(IA)=IPHAR(IA)
  457.   380 CONTINUE
  458. C
  459. C     creation de bb,ibin,igeo
  460. C
  461.       DO 450 ICH=1,2
  462.         MCHPOI=IPO(ICH)
  463. c*        SEGACT MCHPOI
  464.         NSOUPO=IPCHP(/1)
  465.         DO 430 IA=1,NSOUPO
  466.           MSOUPO=IPCHP(IA)
  467.           SEGACT MSOUPO
  468.           MELEME=IGEOC
  469.           SEGACT MELEME
  470.           MPOVAL=IPOVAL
  471.           SEGACT MPOVAL
  472.           NBELEM=NUM(/2)
  473.           N=VPOCHA(/1)
  474.           NC=VPOCHA(/2)
  475.           NC1=NOCOMP(/2)
  476. C
  477.          DO 420 IB=1,NC1
  478.            DO 390 IC=1,NNIN
  479.              IF(NOCOMP(IB).NE.IPCOM(IC)) GO TO 390
  480.              IF(NOHARM(IB).EQ.IPHAR(IC)) GO TO 400
  481.   390      CONTINUE
  482.   400      CONTINUE
  483.            DO 411 ID=1,NBELEM
  484.              DMOY(IB)=DMOY(IB)+ABS(VPOCHA(ID,IB)/NBELEM)
  485.   411      CONTINUE
  486.            DO 410 ID=1,NBELEM
  487.               KI=NOPOIN(NUM(1,ID))
  488.               IGEO(KI)=NUM(1,ID)
  489.               IF (  IBIN(IC,KI) .EQ. 1) THEN
  490. C          la valeur au point est defini dans les deux champs
  491.                XX1 = BB(IC,KI)
  492.                XX2 = VPOCHA(ID,IB)
  493.                DXX = ABS ( XX2 - XX1 )
  494.                SXX = DMOY(IB)
  495.                IF ( DXX .LE. (1.D-4*SXX).or.noer.eq.1) THEN
  496.                   BB(IC,KI) = ( XX1 + XX2 ) / 2.D0
  497.                ELSE
  498. C          les valeurs des champs au meme point sont différentes
  499.                 IF (IIMPI.EQ.123)
  500.      &            write (IOIMP,*) xx1,xx2,sxx,DXX
  501.                 CALL ERREUR (651)
  502.                 RETURN
  503.                 SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  504.                 GOTO 666
  505.                ENDIF
  506.              ELSE
  507.                BB(IC,KI)=BB(IC,KI)+VPOCHA(ID,IB)
  508. *              DMOY(IC) = DMOY(IC) +ABS(BB(IC,KI)/NNNOE)
  509.              ENDIF
  510.              IBIN(IC,KI)=1
  511.   410      CONTINUE
  512.   420    CONTINUE
  513.   430  CONTINUE
  514.   450 CONTINUE
  515.       ITRAV=MTRAV
  516. C
  517. C     reconstuction de la partition
  518. C
  519.       CALL CRECHP(ITRAV,ICHPOI)
  520. C
  521.       SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  522.       IRET=ICHPOI
  523.       MCHPOI=ICHPOI
  524. c*      MCHPO1 = IP1
  525. c*      MCHPO2 = IP2
  526. c*      NAT1 = MCHPO1.JATTRI(/1)
  527. c*      NAT2 = MCHPO2.JATTRI(/1)
  528.       NAT=MAX(NAT1,NAT2,1)
  529.       NSOUPO = IPCHP(/1)
  530.       SEGADJ MCHPOI
  531.       IRET=MCHPOI
  532.       MTYPOI=MOT
  533.       IFOPOI = ifos
  534.       JATTRI(1) = 1
  535. C
  536.  666  CONTINUE
  537. c*      RETURN
  538.       END
  539.  
  540.  

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