Télécharger la plaquette Cast3M


fr en

Télécharger tciso.eso

Retour à la liste

Numérotation des lignes : 

tciso

  1. C TCISO     SOURCE    PV090527  26/09/08    21:15:11     12640          
  2. C
  3.       SUBROUTINE TCISO(VCHC,XX,YY,zz,VV,NPT,NISO)
  4. C
  5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  6. C                                                                 C
  7. C         TRACER DES ISOVALEURS D UN CHAMPOINT                    C
  8. C                PAR COLORIAGE DE ZONE                            C
  9. C                LE NOEUD 1 EST CENTRAL A L'ELEMENT               C
  10. C         TRAITE LE CAS OU TOUT L'ELEMENT EST D'UNE COULEUR       C
  11. C         APPELE TRISO SINON                                      C
  12. C                                                                 C
  13. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  14. C
  15.       REAL VCHC
  16. C
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCREEL
  21. -INC CCGEOME
  22. -INC CCTRACE
  23. C
  24.       PARAMETER (NTR=64)
  25.       LOGICAL RANGE,IDEP
  26.       DIMENSION VCHC(*),XTR(NTR),YTR(NTR),ZTR(NTR),VVN(NTR)
  27.       DIMENSION ROUG(NTR),VERT(NTR),BLEU(NTR)
  28.       dimension vvo(3)
  29.       dimension xx(*),yy(*),zz(*),vv(*)
  30.       real*8 vdiff,up,upos,xxx
  31. *      RANGE(XXX)= XXX.GE.-0.000001.AND.XXX.LE.1.000001
  32.       RANGE(XXX)= XXX.GE.(0.d0-xszpre).AND.XXX.LE.(1.d0+xszpre)
  33.  
  34. *      write(ioimp,*) 'coucou tciso, npt,niso=',npt,niso
  35. *      WRITE (IOIMP,9111) (XX(I),YY(I),ZZ(I),VV(I),I=1,NPT)
  36. * 9111 FORMAT(5(2X,4E12.5))
  37.  
  38.  
  39.       VSTART= -xsgran
  40.       VFINAL= xsgran
  41.       VALHAU=VSTART
  42.       if (iogra.eq.6) then
  43.          valbas=vchc(1)
  44. *       valhau=vchc(niso)
  45.          valhau=vchc(max(niso-1,1))
  46.          do 300 i=1,npt
  47.             vvn(i)=(vv(i)-valbas)/(valhau-valbas)
  48.  300     continue
  49.          xtr(1)=xx(1)
  50.          ytr(1)=yy(1)
  51.          ztr(1)=zz(1)
  52.          vvo(1)=vvn(1)
  53.          do 310 ipt=2,npt
  54.             ipn=ipt+1
  55.             if (ipn.gt.npt) ipn=2
  56.             xtr(2)=xx(ipt)
  57.             ytr(2)=yy(ipt)
  58.             ztr(2)=zz(ipt)
  59.             vvo(2)=vvn(ipt)
  60.             xtr(3)=xx(ipn)
  61.             ytr(3)=yy(ipn)
  62.             ztr(3)=zz(ipn)
  63.             vvo(3)=vvn(ipn)
  64. C           Calcul des couleurs depuis la palette
  65.             DO I=1,3
  66.               XALFA=VVO(I)
  67.               CALL PALET3(XALFA,IPALET,R1,V1,B1)
  68.               ROUG(I)=R1
  69.               VERT(I)=V1
  70.               BLEU(I)=B1
  71.             ENDDO
  72.             CALL OGLTRISO(XTR,YTR,ZTR,VVO,3,ROUG,VERT,BLEU)
  73. *       WRITE(IOIMP,*) ' coul ',vvn(1),vvn(2),vvn(3)
  74.  310     continue
  75.       endif
  76.       IF (ISOTYP.GT.0.and.iogra.ne.6) THEN
  77.          DO 50 KK=1,NISO
  78.             VALBAS=VALHAU
  79.             VALHAU=VFINAL
  80. *      IF (KK.NE.NISO) VALHAU=(VCHC(KK)+VCHC(KK+1))/2
  81. *      TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+5
  82.             IF (KK.NE.NISO) VALHAU=VCHC(KK)
  83. *      TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))/1e+5
  84.             TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))*xszpre
  85. *      toll=max(REAL(XPETIT),toll)
  86.             toll=max(xspeti,toll)
  87.             NP=0
  88. C  VALBAS ET VALHAU SONT LES FRONTIERES DE LA ZONE A COLORIER
  89. C  JE CRAINS QU'IL FAILLE RECENSER LES CAS POSSIBLES
  90. C  LES POINT EXTERIEURS SONT ILS TOUS DANS LA ZONE ?
  91.             np=0
  92. *     SG : 20160727 : je me suis permis de rajouter ce branchement au
  93. *     label 11 directement car meme si les points exterieurs sont tous
  94. *     en zone, cela n'est pas forcement le cas du noeud central car :
  95. *     - soit la valeur du noeud central est independante des valeurs
  96. *      exterieures (cas des quafs TRI7,QUA9)
  97. *     - soit la valeur du noeud central est une combinaison non convexe 
  98. *       des valeurs exterieures (cas du QUA8).
  99.             goto 11
  100.             do 10 ipt=2,npt
  101.                IF ((VALBAS-toll).LE.VV(IPT).AND.(VALHAU+toll).GE.VV(IPT)
  102.      $              ) THEN
  103.                   NP=NP+1
  104.                   XTR(NP)=XX(IPT)
  105.                   YTR(NP)=YY(IPT)
  106.                ELSE
  107.                   np=0
  108.                   goto 11
  109.                ENDIF
  110.  10         continue
  111.             if (niso.lt.16) then
  112. *          CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8)))
  113.                CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO)))
  114.             else
  115.                CALL TRAISO(NP,XTR,YTR,KK)
  116.             endif
  117.             goto 51
  118.  11         CONTINUE
  119. *            write(ioimp,*) 'un des points nest pas dans la zone'
  120. C un des points n'est pas dans la zone
  121. C 1 est il dedans  ?
  122.             IF ((VALBAS-TOLL).LE.VV(1).AND.(VALHAU+TOLL).GE.VV(1)) THEN
  123.                DO 20 ipt=2,npt
  124.                   iptn=ipt+1
  125.                   if (iptn.gt.npt) iptn=2
  126. C IPT est il dedans
  127.                   IF ((VALBAS-toll).LE.VV(IPT).AND.(VALHAU+toll).GE
  128.      $                 .VV(IPT)) THEN
  129.                      NP=NP+1
  130.                      XTR(NP)=XX(IPT)
  131.                      YTR(NP)=YY(IPT)
  132.                   ELSE
  133. C SI IPT EST DEDANS INUTILE DE TESTER LE RAYON
  134.                      vdiff=sign(max(toll,abs(vv(ipt)-vv(1))),vv(ipt)
  135.      $                    -vv(1))
  136.                      UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff)
  137.                      UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff)
  138.                      UP=MAX(UPOSB,UPOSH)
  139.                   up=max(-2*abs(vdiff),up)
  140.                   up=min(2*abs(vdiff),up)
  141.                   UP=UP/abs(VDIFF)
  142.                      IF (RANGE(UP)) THEN
  143.                         NP=NP+1
  144.                         XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1))
  145.                         YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1))
  146.                      ELSE
  147.                         UP=MIN(UPOSB,UPOSH)
  148.                   up=max(-2*abs(vdiff),up)
  149.                   up=min(2*abs(vdiff),up)
  150.                   UP=UP/abs(VDIFF)
  151.                         IF (RANGE(UP)) THEN
  152.                            NP=NP+1
  153.                            XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1))
  154.                            YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1))
  155.                         ENDIF
  156.                      ENDIF
  157.                   ENDIF
  158.                   vdiff=sign(max(toll,abs(vv(iptn)-vv(ipt))),vv(iptn)
  159.      $                 -vv(ipt))
  160.                   UPOSH=(VALHAU+TOLL-VV(ipt))*sign(1.d0,vdiff)
  161.                   UPOSB=(VALBAS-TOLL-VV(ipt))*sign(1.d0,vdiff)
  162.                   UP=MIN(UPOSB,UPOSH)
  163.                   up=max(-2*abs(vdiff),up)
  164.                   up=min(2*abs(vdiff),up)
  165.                   UP=UP/abs(VDIFF)
  166.                   IF (RANGE(UP)) THEN
  167.                      NP=NP+1
  168.                      XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt))
  169.                      YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt))
  170.                   ENDIF
  171.  
  172.                   UP=MAX(UPOSB,UPOSH)
  173.                   up=max(-2*abs(vdiff),up)
  174.                   up=min(2*abs(vdiff),up)
  175.                   UP=UP/abs(VDIFF)
  176.                   IF (RANGE(UP)) THEN
  177.                      NP=NP+1
  178.                      XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt))
  179.                      YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt))
  180.                   ENDIF
  181.  20            continue
  182.                if (niso.lt.16) then
  183. *          CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8)))
  184.                   CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO)))
  185.                else
  186.                   CALL TRAISO(NP,XTR,YTR,KK)
  187.                endif
  188.                goto 50
  189.             endif
  190. C  1 n'est pas dans la zone !
  191. C on tourne autour de 1 en cherchant un point de depart !
  192. *            write(ioimp,*) '1 nest pas dans la zone'
  193.             iq=0
  194.             if (vv(1).lt.(valbas-toll)) iq=1
  195.             if (vv(1).gt.(valhau+toll)) iq=2
  196. **          if (iq.eq.0) write (ioimp,*) ' prob iq '
  197. * sg : iptf n'etait pas initialisee. On l'initialise a 1 mais on nest
  198. *      pas sur
  199.             iptf=1
  200.             iptft=2
  201.             iptd=npt
  202.  44         continue
  203.             np=0
  204.             do 30 ipt=iptft,npt
  205.                if (iptd.ge.iptf.and.ipt.gt.iptd) goto 50
  206.                if (iq.eq.1.and.vv(ipt).lt.(valbas-toll)) goto 30
  207.                if (iq.eq.2.and.vv(ipt).gt.(valhau+toll)) goto 30
  208.                goto 31
  209.  30         continue
  210. C  pas de point de depart il n'y a rien a faire
  211. *            write(ioimp,*) 'pas de point de depart'
  212.             goto 50
  213.  31         continue
  214.             do 32 irec=1,npt-2
  215.                ipt2=ipt-irec
  216.                if (ipt2.le.1) ipt2=ipt2+npt-1
  217.                if (iq.eq.1.and.vv(ipt2).lt.(valbas-toll)) goto 33
  218.                if (iq.eq.2.and.vv(ipt2).gt.(valhau+toll)) goto 33
  219.  32         continue
  220. **     WRITE(IOIMP,*) ' prob toujours pas de point de depart '
  221.  33         continue
  222.             iptd=ipt2
  223. C  IPTD ne traverse pas et iptd+1 traverse
  224.             do 40 iptb=iptd,iptd+npt-2
  225.                ipt=iptb
  226.                if (ipt.gt.npt) ipt=ipt-npt+1
  227.                iptn=ipt+1
  228.                if (iptn.gt.npt) iptn=iptn-npt+1
  229.                vdiff=sign(max(toll,abs(vv(iptn)-vv(ipt))),vv(iptn)
  230.      $              -vv(ipt))
  231.                UPOSH=(VALHAU+TOLL-VV(ipt))*sign(1.d0,vdiff)
  232.                UPOSB=(VALBAS-TOLL-VV(ipt))*sign(1.d0,vdiff)
  233.                UP=MIN(UPOSB,UPOSH)
  234.                   up=max(-2*abs(vdiff),up)
  235.                   up=min(2*abs(vdiff),up)
  236.                   UP=UP/abs(VDIFF)
  237.                IF (RANGE(UP)) THEN
  238.                   NP=NP+1
  239.                   XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt))
  240.                   YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt))
  241.                ENDIF
  242.                UP=MAX(UPOSB,UPOSH)
  243.                   up=max(-2*abs(vdiff),up)
  244.                   up=min(2*abs(vdiff),up)
  245.                   UP=UP/abs(VDIFF)
  246.                IF (RANGE(UP)) THEN
  247.                   NP=NP+1
  248.                   XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt))
  249.                   YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt))
  250.                ENDIF
  251. C IPTN Y est il ?
  252.                npin=np
  253.                IF ((VALBAS-toll).LE.VV(IPTN).AND.(VALHAU+toll).GE
  254.      $              .VV(IPTN)) THEN
  255.                   NP=NP+1
  256.                   XTR(NP)=XX(IPTN)
  257.                   YTR(NP)=YY(IPTN)
  258.                ELSE
  259. C SI IPTN EST DEDANS INUTILE DE TESTER LE RAYON
  260.                   vdiff=sign(max(toll,abs(vv(iptn)-vv(1))),vv(iptn)-vv(1
  261.      $                 ))
  262.                   UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff)
  263.                   UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff)
  264.                   UP=MAX(UPOSB,UPOSH)
  265.                   up=max(-2*abs(vdiff),up)
  266.                   up=min(2*abs(vdiff),up)
  267.                   UP=UP/abs(VDIFF)
  268.                   IF (RANGE(UP)) THEN
  269.                      NP=NP+1
  270.                      XTR(NP)=XX(1)+UP*(XX(iptn)-XX(1))
  271.                      YTR(NP)=YY(1)+UP*(YY(iptn)-YY(1))
  272.                   ELSE
  273.                      UP=MIN(UPOSB,UPOSH)
  274.                   up=max(-2*abs(vdiff),up)
  275.                   up=min(2*abs(vdiff),up)
  276.                   UP=UP/abs(VDIFF)
  277.                      IF (RANGE(UP)) THEN
  278.                         NP=NP+1
  279.                         XTR(NP)=XX(1)+UP*(XX(iptn)-XX(1))
  280.                         YTR(NP)=YY(1)+UP*(YY(iptn)-YY(1))
  281.                      ENDIF
  282.                   ENDIF
  283.                ENDIF
  284.                if (np.eq.npin) goto 41
  285. C  on est arrive au bout dans ce sens
  286.  40         continue
  287. **      WRITE(IOIMP,*) ' on ne devrai pas passer la', valhau, valbas,toll
  288.  41         continue
  289.             iptf=ipt
  290. C on revient en arriere de iptf a iptd
  291. C  apres iptf, il ne se passe plus rien
  292. C  avant iptd, il ne se passe plus rien
  293.             do 42 iptb=iptf,iptf-npt+1,-1
  294.                ipt=iptb
  295.                if (ipt.le.1) ipt=ipt+npt-1
  296.                if (ipt.eq.iptd) goto 48
  297.                vdiff=sign(max(toll,abs(vv(ipt)-vv(1))),vv(ipt)-vv(1))
  298.                UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff)
  299.                UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff)
  300.                UP=MIN(UPOSB,UPOSH)
  301.                   up=max(-2*abs(vdiff),up)
  302.                   up=min(2*abs(vdiff),up)
  303.                   UP=UP/abs(VDIFF)
  304.                IF (RANGE(UP)) THEN
  305.                   NP=NP+1
  306.                   XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1))
  307.                   YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1))
  308.                ELSE
  309.                   UP=MAX(UPOSB,UPOSH)
  310.                   up=max(-2*abs(vdiff),up)
  311.                   up=min(2*abs(vdiff),up)
  312.                   UP=UP/abs(VDIFF)
  313.                   IF (RANGE(UP)) THEN
  314.                      NP=NP+1
  315.                      XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1))
  316.                      YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1))
  317.                   ENDIF
  318.                ENDIF
  319.  42         continue
  320.  48         continue
  321.             if (niso.lt.16) then
  322. *          CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8)))
  323.                CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO)))
  324.             else
  325.                CALL TRAISO(NP,XTR,YTR,KK)
  326.             endif
  327.             iptft=max(iptft,iptf+1)
  328.             if (iptft.lt.npt) goto 44
  329.  50      CONTINUE
  330.  51      CONTINUE
  331.          IF (ISOTYP.EQ.2) THEN
  332.             call chcoul(IDNOIR)
  333.             DO 250 KK=1,NISO-1
  334. *      TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+5
  335. *       VALDES = (VCHC(KK)+VCHC(KK+1))/2
  336. *      TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))/1e+5
  337.                TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))*xszpre
  338.                VALDES = VCHC(KK)
  339.                do 220 ipt=2,npt
  340.                   NP=0
  341.                   iptn=ipt+1
  342.                   if (iptn.gt.npt) iptn=2
  343.                   UPOS=-1.
  344.                   IF (ABS(VV(iptn)-VV(ipt)).GT.TOLL)
  345.      *                 UPOS=(VALDES-VV(ipt))/(VV(iptn)-VV(ipt))
  346.                   IF (RANGE(UPOS)) THEN
  347.                      NP=NP+1
  348.                      XTR(NP)=XX(ipt)+UPOS*(XX(iptn)-XX(ipt))
  349.                      YTR(NP)=YY(ipt)+UPOS*(YY(iptn)-YY(ipt))
  350.                   ENDIF
  351.                   UPOS=-1.
  352.                   IF (ABS(VV(ipt)-VV(1)).GT.TOLL)
  353.      *                 UPOS=(VALDES-VV(1))/(VV(ipt)-VV(1))
  354.                   IF (RANGE(UPOS)) THEN
  355.                      NP=NP+1
  356.                      XTR(NP)=XX(1)+UPOS*(XX(ipt)-XX(1))
  357.                      YTR(NP)=YY(1)+UPOS*(YY(ipt)-YY(1))
  358.                   ENDIF
  359.                   UPOS=-1.
  360.                   IF (ABS(VV(iptn)-VV(1)).GT.TOLL)
  361.      *                 UPOS=(VALDES-VV(1))/(VV(iptn)-VV(1))
  362.                   IF (RANGE(UPOS)) THEN
  363.                      NP=NP+1
  364.                      XTR(NP)=XX(1)+UPOS*(XX(iptn)-XX(1))
  365.                      YTR(NP)=YY(1)+UPOS*(YY(iptn)-YY(1))
  366.                   ENDIF
  367. *  il convient de fermer la ligne
  368.                   if (np.gt.2) then
  369.                      np=np+1
  370.                      xtr(np)=xtr(1)
  371.                      ytr(np)=ytr(1)
  372.                   endif
  373.                   if (np.gt.1) call polrl(np,xtr,ytr,ztr)
  374.  220           continue
  375.  250        CONTINUE
  376.          ENDIF
  377.       ELSEIF (iogra.ne.6) THEN
  378.          DO 150 KK=1,NISO-1
  379. *      TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+3
  380.             TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))*xszpre
  381.             VALDES = VCHC(KK)
  382.             NP=0
  383.             do 270 ipt=2,npt
  384.                iptn=ipt+1
  385.                if (iptn.gt.npt) iptn=2
  386.                UPOS=-1.
  387.                IF (ABS(VV(iptn)-VV(ipt)).GT.TOLL)
  388.      *              UPOS=(VALDES-VV(ipt))/(VV(iptn)-VV(ipt))
  389.                IF (RANGE(UPOS)) THEN
  390.                   NP=NP+1
  391.                   XTR(NP)=XX(ipt)+UPOS*(XX(iptn)-XX(ipt))
  392.                   YTR(NP)=YY(ipt)+UPOS*(YY(iptn)-YY(ipt))
  393.                ENDIF
  394.                UPOS=-1.
  395.                IF (ABS(VV(ipt)-VV(1)).GT.TOLL)
  396.      *              UPOS=(VALDES-VV(1))/(VV(ipt)-VV(1))
  397.                IF (RANGE(UPOS)) THEN
  398.                   NP=NP+1
  399.                   XTR(NP)=XX(1)+UPOS*(XX(ipt)-XX(1))
  400.                   YTR(NP)=YY(1)+UPOS*(YY(ipt)-YY(1))
  401.                ENDIF
  402.                UPOS=-1.
  403.                IF (ABS(VV(iptn)-VV(1)).GT.TOLL)
  404.      *              UPOS=(VALDES-VV(1))/(VV(iptn)-VV(1))
  405.                IF (RANGE(UPOS)) THEN
  406.                   NP=NP+1
  407.                   XTR(NP)=XX(1)+UPOS*(XX(iptn)-XX(1))
  408.                   YTR(NP)=YY(1)+UPOS*(YY(iptn)-YY(1))
  409.                ENDIF
  410. *  il convient de fermer la ligne
  411.                if (np.gt.2) then
  412.                   np=np+1
  413.                   xtr(np)=xtr(1)
  414.                   ytr(np)=ytr(1)
  415.                endif
  416.                if (np.gt.1) then
  417. *sg                  if (niso.lt.16) then
  418.                   if (niso.lt.13) then
  419. *sg                     CALL CHCOUL(ICOTAB(KK*(2-NISO/8)))
  420.                      CALL CHCOUL(ICOTAB(ISOTA0(KK,NISO)))
  421.                   else
  422.                      CALL CHCOUL(ICOTAB(MOD(KK,12)+1))
  423. *sg                     CALL CHCOUL(KK)
  424.                   endif
  425.                   call polrl(np,xtr,ytr,ztr)
  426.                endif
  427.  270        continue
  428.  150     CONTINUE
  429.       ENDIF
  430.       RETURN
  431.       END
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  

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