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

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