C TCISO SOURCE PV 21/11/02 21:15:27 11158 C SUBROUTINE TCISO(VCHC,XX,YY,zz,VV,NPT,NISO) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C TRACER DES ISOVALEURS D UN CHAMPOINT C C PAR COLORIAGE DE ZONE C C LE NOEUD 1 EST CENTRAL A L'ELEMENT C C TRAITE LE CAS OU TOUT L'ELEMENT EST D'UNE COULEUR C C APPELE TRISO SINON C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C REAL VCHC C -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCGEOME -INC CCTRACE C PARAMETER (NTR=64) LOGICAL RANGE,IDEP DIMENSION VCHC(*),XTR(NTR),YTR(NTR),ZTR(NTR),VVN(NTR) dimension vvo(3) dimension xx(*),yy(*),zz(*),vv(*) real*8 vdiff,up,upos,xxx * RANGE(XXX)= XXX.GE.-0.000001.AND.XXX.LE.1.000001 RANGE(XXX)= XXX.GE.(0.d0-xszpre).AND.XXX.LE.(1.d0+xszpre) * write(ioimp,*) 'coucou tciso, npt,niso=',npt,niso * WRITE (IOIMP,9111) (XX(I),YY(I),ZZ(I),VV(I),I=1,NPT) * 9111 FORMAT(5(2X,4E12.5)) VSTART= -xsgran VFINAL= xsgran VALHAU=VSTART if (iogra.eq.6) then valbas=vchc(1) * valhau=vchc(niso) valhau=vchc(max(niso-1,1)) do 300 i=1,npt vvn(i)=(vv(i)-valbas)/(valhau-valbas) 300 continue xtr(1)=xx(1) ytr(1)=yy(1) ztr(1)=zz(1) vvo(1)=vvn(1) do 310 ipt=2,npt ipn=ipt+1 if (ipn.gt.npt) ipn=2 xtr(2)=xx(ipt) ytr(2)=yy(ipt) ztr(2)=zz(ipt) vvo(2)=vvn(ipt) xtr(3)=xx(ipn) ytr(3)=yy(ipn) ztr(3)=zz(ipn) vvo(3)=vvn(ipn) call ogltriso(xtr,ytr,ztr,vvo,3) * WRITE(IOIMP,*) ' coul ',vvn(1),vvn(2),vvn(3) 310 continue endif IF (ISOTYP.GT.0.and.iogra.ne.6) THEN DO 50 KK=1,NISO VALBAS=VALHAU VALHAU=VFINAL * IF (KK.NE.NISO) VALHAU=(VCHC(KK)+VCHC(KK+1))/2 * TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+5 IF (KK.NE.NISO) VALHAU=VCHC(KK) * TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))/1e+5 TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))*xszpre * toll=max(REAL(XPETIT),toll) toll=max(xspeti,toll) NP=0 C VALBAS ET VALHAU SONT LES FRONTIERES DE LA ZONE A COLORIER C JE CRAINS QU'IL FAILLE RECENSER LES CAS POSSIBLES C LES POINT EXTERIEURS SONT ILS TOUS DANS LA ZONE ? np=0 * SG : 20160727 : je me suis permis de rajouter ce branchement au * label 11 directement car meme si les points exterieurs sont tous * en zone, cela n'est pas forcement le cas du noeud central car : * - soit la valeur du noeud central est independante des valeurs * exterieures (cas des quafs TRI7,QUA9) * - soit la valeur du noeud central est une combinaison non convexe * des valeurs exterieures (cas du QUA8). goto 11 do 10 ipt=2,npt IF ((VALBAS-toll).LE.VV(IPT).AND.(VALHAU+toll).GE.VV(IPT) $ ) THEN NP=NP+1 XTR(NP)=XX(IPT) YTR(NP)=YY(IPT) ELSE np=0 goto 11 ENDIF 10 continue if (niso.lt.16) then * CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8))) CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO))) else CALL TRAISO(NP,XTR,YTR,KK) endif goto 51 11 CONTINUE * write(ioimp,*) 'un des points nest pas dans la zone' C un des points n'est pas dans la zone C 1 est il dedans ? IF ((VALBAS-TOLL).LE.VV(1).AND.(VALHAU+TOLL).GE.VV(1)) THEN DO 20 ipt=2,npt iptn=ipt+1 if (iptn.gt.npt) iptn=2 C IPT est il dedans IF ((VALBAS-toll).LE.VV(IPT).AND.(VALHAU+toll).GE $ .VV(IPT)) THEN NP=NP+1 XTR(NP)=XX(IPT) YTR(NP)=YY(IPT) ELSE C SI IPT EST DEDANS INUTILE DE TESTER LE RAYON vdiff=sign(max(toll,abs(vv(ipt)-vv(1))),vv(ipt) $ -vv(1)) UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff) UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff) UP=MAX(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1)) ELSE UP=MIN(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1)) ENDIF ENDIF ENDIF vdiff=sign(max(toll,abs(vv(iptn)-vv(ipt))),vv(iptn) $ -vv(ipt)) UPOSH=(VALHAU+TOLL-VV(ipt))*sign(1.d0,vdiff) UPOSB=(VALBAS-TOLL-VV(ipt))*sign(1.d0,vdiff) UP=MIN(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt)) ENDIF UP=MAX(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt)) ENDIF 20 continue if (niso.lt.16) then * CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8))) CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO))) else CALL TRAISO(NP,XTR,YTR,KK) endif goto 50 endif C 1 n'est pas dans la zone ! C on tourne autour de 1 en cherchant un point de depart ! * write(ioimp,*) '1 nest pas dans la zone' iq=0 if (vv(1).lt.(valbas-toll)) iq=1 if (vv(1).gt.(valhau+toll)) iq=2 ** if (iq.eq.0) write (ioimp,*) ' prob iq ' * sg : iptf n'etait pas initialisee. On l'initialise a 1 mais on nest * pas sur iptf=1 iptft=2 iptd=npt 44 continue np=0 do 30 ipt=iptft,npt if (iptd.ge.iptf.and.ipt.gt.iptd) goto 50 if (iq.eq.1.and.vv(ipt).lt.(valbas-toll)) goto 30 if (iq.eq.2.and.vv(ipt).gt.(valhau+toll)) goto 30 goto 31 30 continue C pas de point de depart il n'y a rien a faire * write(ioimp,*) 'pas de point de depart' goto 50 31 continue do 32 irec=1,npt-2 ipt2=ipt-irec if (ipt2.le.1) ipt2=ipt2+npt-1 if (iq.eq.1.and.vv(ipt2).lt.(valbas-toll)) goto 33 if (iq.eq.2.and.vv(ipt2).gt.(valhau+toll)) goto 33 32 continue ** WRITE(IOIMP,*) ' prob toujours pas de point de depart ' 33 continue iptd=ipt2 C IPTD ne traverse pas et iptd+1 traverse do 40 iptb=iptd,iptd+npt-2 ipt=iptb if (ipt.gt.npt) ipt=ipt-npt+1 iptn=ipt+1 if (iptn.gt.npt) iptn=iptn-npt+1 vdiff=sign(max(toll,abs(vv(iptn)-vv(ipt))),vv(iptn) $ -vv(ipt)) UPOSH=(VALHAU+TOLL-VV(ipt))*sign(1.d0,vdiff) UPOSB=(VALBAS-TOLL-VV(ipt))*sign(1.d0,vdiff) UP=MIN(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt)) ENDIF UP=MAX(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UP*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UP*(YY(iptn)-YY(ipt)) ENDIF C IPTN Y est il ? npin=np IF ((VALBAS-toll).LE.VV(IPTN).AND.(VALHAU+toll).GE $ .VV(IPTN)) THEN NP=NP+1 XTR(NP)=XX(IPTN) YTR(NP)=YY(IPTN) ELSE C SI IPTN EST DEDANS INUTILE DE TESTER LE RAYON vdiff=sign(max(toll,abs(vv(iptn)-vv(1))),vv(iptn)-vv(1 $ )) UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff) UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff) UP=MAX(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(iptn)-XX(1)) YTR(NP)=YY(1)+UP*(YY(iptn)-YY(1)) ELSE UP=MIN(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(iptn)-XX(1)) YTR(NP)=YY(1)+UP*(YY(iptn)-YY(1)) ENDIF ENDIF ENDIF if (np.eq.npin) goto 41 C on est arrive au bout dans ce sens 40 continue ** WRITE(IOIMP,*) ' on ne devrai pas passer la', valhau, valbas,toll 41 continue iptf=ipt C on revient en arriere de iptf a iptd C apres iptf, il ne se passe plus rien C avant iptd, il ne se passe plus rien do 42 iptb=iptf,iptf-npt+1,-1 ipt=iptb if (ipt.le.1) ipt=ipt+npt-1 if (ipt.eq.iptd) goto 48 vdiff=sign(max(toll,abs(vv(ipt)-vv(1))),vv(ipt)-vv(1)) UPOSH=(VALHAU+TOLL-VV(1))*sign(1.d0,vdiff) UPOSB=(VALBAS-TOLL-VV(1))*sign(1.d0,vdiff) UP=MIN(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1)) ELSE UP=MAX(UPOSB,UPOSH) up=max(-2*abs(vdiff),up) up=min(2*abs(vdiff),up) UP=UP/abs(VDIFF) IF (RANGE(UP)) THEN NP=NP+1 XTR(NP)=XX(1)+UP*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UP*(YY(ipt)-YY(1)) ENDIF ENDIF 42 continue 48 continue if (niso.lt.16) then * CALL TRAISO(NP,XTR,YTR,ICOTAB(KK*(2-NISO/8))) CALL TRAISO(NP,XTR,YTR,ICOTAB(ISOTAB(KK,NISO))) else CALL TRAISO(NP,XTR,YTR,KK) endif iptft=max(iptft,iptf+1) if (iptft.lt.npt) goto 44 50 CONTINUE 51 CONTINUE IF (ISOTYP.EQ.2) THEN call chcoul(IDNOIR) DO 250 KK=1,NISO-1 * TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+5 * VALDES = (VCHC(KK)+VCHC(KK+1))/2 * TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))/1e+5 TOLL=ABS(VCHC(min(niso-1,KK+1))-VCHC(max(1,KK-1)))*xszpre VALDES = VCHC(KK) do 220 ipt=2,npt NP=0 iptn=ipt+1 if (iptn.gt.npt) iptn=2 UPOS=-1. IF (ABS(VV(iptn)-VV(ipt)).GT.TOLL) * UPOS=(VALDES-VV(ipt))/(VV(iptn)-VV(ipt)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UPOS*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UPOS*(YY(iptn)-YY(ipt)) ENDIF UPOS=-1. IF (ABS(VV(ipt)-VV(1)).GT.TOLL) * UPOS=(VALDES-VV(1))/(VV(ipt)-VV(1)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(1)+UPOS*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UPOS*(YY(ipt)-YY(1)) ENDIF UPOS=-1. IF (ABS(VV(iptn)-VV(1)).GT.TOLL) * UPOS=(VALDES-VV(1))/(VV(iptn)-VV(1)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(1)+UPOS*(XX(iptn)-XX(1)) YTR(NP)=YY(1)+UPOS*(YY(iptn)-YY(1)) ENDIF * il convient de fermer la ligne if (np.gt.2) then np=np+1 xtr(np)=xtr(1) ytr(np)=ytr(1) endif if (np.gt.1) call polrl(np,xtr,ytr,ztr) 220 continue 250 CONTINUE ENDIF ELSEIF (iogra.ne.6) THEN DO 150 KK=1,NISO-1 * TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))/1e+3 TOLL=ABS(VCHC(min(niso,KK+1))-VCHC(max(1,KK-1)))*xszpre VALDES = VCHC(KK) NP=0 do 270 ipt=2,npt iptn=ipt+1 if (iptn.gt.npt) iptn=2 UPOS=-1. IF (ABS(VV(iptn)-VV(ipt)).GT.TOLL) * UPOS=(VALDES-VV(ipt))/(VV(iptn)-VV(ipt)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(ipt)+UPOS*(XX(iptn)-XX(ipt)) YTR(NP)=YY(ipt)+UPOS*(YY(iptn)-YY(ipt)) ENDIF UPOS=-1. IF (ABS(VV(ipt)-VV(1)).GT.TOLL) * UPOS=(VALDES-VV(1))/(VV(ipt)-VV(1)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(1)+UPOS*(XX(ipt)-XX(1)) YTR(NP)=YY(1)+UPOS*(YY(ipt)-YY(1)) ENDIF UPOS=-1. IF (ABS(VV(iptn)-VV(1)).GT.TOLL) * UPOS=(VALDES-VV(1))/(VV(iptn)-VV(1)) IF (RANGE(UPOS)) THEN NP=NP+1 XTR(NP)=XX(1)+UPOS*(XX(iptn)-XX(1)) YTR(NP)=YY(1)+UPOS*(YY(iptn)-YY(1)) ENDIF * il convient de fermer la ligne if (np.gt.2) then np=np+1 xtr(np)=xtr(1) ytr(np)=ytr(1) endif if (np.gt.1) then *sg if (niso.lt.16) then if (niso.lt.13) then *sg CALL CHCOUL(ICOTAB(KK*(2-NISO/8))) CALL CHCOUL(ICOTAB(ISOTA0(KK,NISO))) else CALL CHCOUL(ICOTAB(MOD(KK,12)+1)) *sg CALL CHCOUL(KK) endif call polrl(np,xtr,ytr,ztr) endif 270 continue 150 CONTINUE ENDIF RETURN END