C CADRCL SOURCE CB215821 23/01/25 21:15:06 11573 C TRACE DE DEFORMES C CALCUL DU CADRE GLOBAL C PROJECTION SUCCESSIVE DE CHAQUE DEFORME (SELON ICLE) C 1995 option DIRE P.PEGON JRC-ISPRA C C PP option DIRE SUBROUTINE CADRCL(KABCOR,LABCO2,IOEIL,XPROJ, # ICLE,XMINT,YMINT,XMAXT,YMAXT,zmint,zmaxt,cgrav,diloc,LDIRE, > axez) IMPLICIT INTEGER(I-N) DIMENSION XO(3),XP(3),XN(3),XG(3),XV(3),UI(3),UJ(3),axez(3) dimension cgrav(*) C+PP option DIRE dimension diloc(*) C+PP COMMON /CCADRC/XN,XG,XO,UI,UJ C ATTENTION SEGMENT KABCOR(0) SEGMENT KABCO2(2,0) SEGMENT LABCO2(3,0) SEGMENT ICOR2(0) SEGMENT KXPRO2(NVEC) SEGMENT SXCORD REAL XCORD(IDIM,ITE) ENDSEGMENT SEGMENT SXCOR2 REAL XCOR2(IDIM,ITE) ENDSEGMENT -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC CCREEL SEGMENT XPROJ(3,ITE) SEGMENT XPRO2(3,ITE) C+PP option DIRE REAL*8 SSN LOGICAL LDIRE *dbg write(6,*) 'coucou cadrcl' C+PP C ICLE = 0 RECHERCHE DES MIN MAX SUR TOUTES LES DEFORMES C CALCUL DE LA PROJECTION XMINT=1E30 YMINT=XMINT ZMINT=XMINT XMAXT=-1E30 YMAXT=XMAXT ZMAXT=XMAXT NCORD=KABCOR(/1) IF (IDIM.EQ.3) GOTO 50 DO 10 ICORD=1,NCORD SXCORD=KABCOR(ICORD) * write(6,*) '1 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) DO 19 J=1,ITE XMINT=MIN(XMINT,XCORD(1,J)) XMAXT=MAX(XMAXT,XCORD(1,J)) YMINT=MIN(YMINT,XCORD(2,J)) YMAXT=MAX(YMAXT,XCORD(2,J)) ZMINT=0 ZMAXT=0 19 CONTINUE 10 CONTINUE IF (ICLE.NE.0) GOTO 1000 IF (LABCO2.EQ.0) GOTO 60 DO 61 ICORD=1,NCORD IF (LABCO2(3,ICORD).EQ.0) GOTO 61 KABCO2=LABCO2(1,ICORD) SXCORD=KABCOR(ICORD) * write(6,*) '2 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) NVEC=KABCO2(/2) IF (NVEC.EQ.0) GOTO 61 DO 7 IVEC=1,NVEC SXCOR2=KABCO2(1,IVEC) ICOR2=KABCO2(2,IVEC) ITE=XCOR2(/2) DO 18 J=1,ITE IF (ICOR2(J).EQ.0) GOTO 18 UX=XCOR2(1,J)-XCORD(1,J) UY=XCOR2(2,J)-XCORD(2,J) U1=XCOR2(1,J)-UX/3-UY/5 V1=XCOR2(2,J)-UY/3+UX/5 XMINT=MIN(XMINT,U1) XMAXT=MAX(XMAXT,U1) YMINT=MIN(YMINT,V1) YMAXT=MAX(YMAXT,V1) U1=XCOR2(1,J)-UX/3+UY/5 V1=XCOR2(2,J)-UY/3-UX/5 XMINT=MIN(XMINT,U1) XMAXT=MAX(XMAXT,U1) YMINT=MIN(YMINT,V1) YMAXT=MAX(YMAXT,V1) XMINT=MIN(XMINT,XCOR2(1,J)) XMAXT=MAX(XMAXT,XCOR2(1,J)) YMINT=MIN(YMINT,XCOR2(2,J)) YMAXT=MAX(YMAXT,XCOR2(2,J)) 18 CONTINUE 7 CONTINUE 61 CONTINUE 60 CONTINUE * pour eviter des calculs sur des quantites non significatives, * on agrandi un peu le cadre XDIFF=XMAXT-XMINT YDIFF=YMAXT-YMINT ZDIFF=ZMAXT-ZMINT XSOMM=XMAXT+XMINT YSOMM=YMAXT+YMINT ZSOMM=ZMAXT+ZMINT XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10) YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10) ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF) XMINT=(XSOMM-XDIFF)/2 XMAXT=(XSOMM+XDIFF)/2 YMINT=(YSOMM-YDIFF)/2 YMAXT=(YSOMM+YDIFF)/2 ZMINT=(ZSOMM-ZDIFF)/2 ZMAXT=(ZSOMM+ZDIFF)/2 * write(6,*) '1 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt IF (ICLE.NE.0) GOTO 1000 RETURN 50 CONTINUE XMINT=1E30 XMAXT=-XMINT YMINT=XMINT YMAXT=-YMINT ZMINT=XMINT ZMAXT=-YMINT IREF=(IOEIL-1)*4 SEGACT,MCOORD XO(1)=XCOOR(IREF+1) XO(2)=XCOOR(IREF+2) XO(3)=XCOOR(IREF+3) C POINT MOYEN DO 1 I=1,3 XG(I)=0.D0 1 CONTINUE NPTOT=0 DO 100 ICORD=1,NCORD SXCORD=KABCOR(ICORD) * write(6,*) '3 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) NPTOT=NPTOT+ITE DO 2 I=1,ITE DO 3 J=1,3 XG(J)=XG(J)+XCORD(J,I) 3 CONTINUE 2 CONTINUE 100 CONTINUE C+PP option DIRE IF (LDIRE)THEN DO J=1,3 XG(J)=cgrav(J) XN(J)=XO(J)-XG(J) ENDDO ELSE C+PP DO 4 J=1,3 XG(J)=XG(J)/NPTOT XN(J)=XO(J)-XG(J) cgrav(j)=xg(j) 4 CONTINUE C+PP option DIRE ENDIF C+PP C NORMALE AU PLAN SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2) SN=SN*SNI IF (SN.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN C+PP option DIRE SSN=SN C+PP DO 5 J=1,3 XN(J)=XN(J)/SN 5 CONTINUE C REPERE LOCAL SUR LE PLAN C+PP option DIRE IF (LDIRE)THEN DO J=1,3 UJ(J)=diloc(J) ENDDO ELSE C+PP UJ(1)=0.D0 UJ(2)=0.D0 UJ(3)=1.D0 C+PP option DIRE ENDIF C+PP 21 CONTINUE SV=UJ(1)*XN(1)+UJ(2)*XN(2)+UJ(3)*XN(3) DO 20 J=1,3 UJ(J)=UJ(J)-SV*XN(J) 20 CONTINUE SV=UJ(1)**2+UJ(2)**2+UJ(3)**2 IF (ABS(SV).LT.0.1) THEN UJ(1)=0.D0 UJ(2)=1.D0 UJ(3)=1.D0 GOTO 21 ENDIF SV=SQRT(SV) UJ(1)=UJ(1)/SV UJ(2)=UJ(2)/SV UJ(3)=UJ(3)/SV UI(1)=UJ(2)*XN(3)-UJ(3)*XN(2) UI(2)=UJ(3)*XN(1)-UJ(1)*XN(3) UI(3)=UJ(1)*XN(2)-UJ(2)*XN(1) C PROJECTION CONIQUE SUR LE PLAN DO 170 ICORD=1,NCORD SXCORD=KABCOR(ICORD) * write(6,*) '4 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) DO 12 I=1,ITE DO 13 J=1,3 XP(J)=XCORD(J,I) XV(J)=XP(J)-XO(J) 13 CONTINUE SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3) C+PP option DIRE IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN C+PP SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3)) $ *XN(3) DO 14 J=1,3 XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J) 14 CONTINUE XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3) YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3) ZPRO=-SN XMINT=MIN(XMINT,XPRO) XMAXT=MAX(XMAXT,XPRO) YMINT=MIN(YMINT,YPRO) YMAXT=MAX(YMAXT,YPRO) ZMINT=MIN(ZMINT,ZPRO) ZMAXT=MAX(ZMAXT,ZPRO) C+PP option DIRE ENDIF C+PP 12 CONTINUE 170 CONTINUE * write(6,*) 'cadrcl en 239 ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt IF (LABCO2.EQ.0) GOTO 11 DO 171 ICORD=1,NCORD IF (LABCO2(3,ICORD).EQ.0) GOTO 171 KABCO2=LABCO2(1,ICORD) SXCORD=KABCOR(ICORD) * write(6,*) '5 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) NVEC=KABCO2(/2) IF (NVEC.EQ.0) GOTO 171 DO 6 IVEC=1,NVEC SXCOR2=KABCO2(1,IVEC) ICOR2=KABCO2(2,IVEC) DO 17 I=1,ITE IF(I.GT.ICOR2(/1)) goto 17 IF (ICOR2(I).EQ.0) GOTO 17 C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO? DO J=1,3 XP(J)=XCORD(J,I) XV(J)=XP(J)-XO(J) ENDDO SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3) C+PP option DIRE IF (LDIRE.AND.-SV.LT.SSN) THEN ICOR2(I)=0 GOTO 17 ENDIF C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO? SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3)) $ *XN(3) DO J=1,3 XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J) ENDDO XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3) YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3) C+PP DO 15 J=1,3 XP(J)=XCOR2(J,I) XV(J)=XP(J)-XO(J) 15 CONTINUE SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3) C+PP option DIRE IF (LDIRE.AND.-SV.LT.SSN) THEN ICOR2(I)=0 GOTO 17 ENDIF C+PP SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3)) $ *XN(3) DO 16 J=1,3 XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J) 16 CONTINUE XPROO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3) YPROO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3) XMINT=MIN(XMINT,XPROO) XMAXT=MAX(XMAXT,XPROO) YMINT=MIN(YMINT,YPROO) YMAXT=MAX(YMAXT,YPROO) UX=XPROO-XPRO UY=YPROO-YPRO U1=XPROO-UX/3-UY/5 V1=YPROO-UY/3+UX/5 XMINT=MIN(XMINT,U1) XMAXT=MAX(XMAXT,V1) U1=XPROO-UX/3+UY/5 V1=YPROO-UY/3-UX/5 XMINT=MIN(XMINT,U1) XMAXT=MAX(XMAXT,V1) 17 CONTINUE 6 CONTINUE 171 CONTINUE 11 CONTINUE * pour eviter des calculs sur des quantites non significatives, * on agrandi un peu le cadre * write(6,*) '2 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt XDIFF=XMAXT-XMINT YDIFF=YMAXT-YMINT ZDIFF=ZMAXT-ZMINT XSOMM=XMAXT+XMINT YSOMM=YMAXT+YMINT ZSOMM=ZMAXT+ZMINT * write(6,*) '4 cadrcl ',xdiff,ydiff,zdiff,xsomm,ysomm,zsomm XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10) YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10) ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF) * write(6,*) '5 cadrcl ',xdiff,ydiff,zdiff XMINT=(XSOMM-XDIFF)/2 XMAXT=(XSOMM+XDIFF)/2 YMINT=(YSOMM-YDIFF)/2 YMAXT=(YSOMM+YDIFF)/2 ZMINT=(ZSOMM-ZDIFF)/2 ZMAXT=(ZSOMM+ZDIFF)/2 * write(6,*) '3 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt * axez pour tourner correctement avec opengl axez(1)=0 axez(2)=uj(3) axez(3)=sqrt(1-uj(3)**2) if (xn(3).lt.0) axez(3)=-axez(3) * write (6,*) ' axez ',axez(1),axez(2),axez(3) IF (ICLE.NE.0) GOTO 1000 RETURN 1000 CONTINUE C ON REMPLIT XPROJ POUR LA DEFORME CONCERNEE SXCORD=KABCOR(ICLE) * write(6,*) '6 cadrcl icle sxcord',icord,sxcord ITE=XCORD(/2) IF (IDIM.EQ.2) GOTO 1100 C+PP option DIRE SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2) SN=SN*SNI SSN=0.95*SSN C+PP C PROJECTION CONIQUE SUR LE PLAN DO 612 I=1,ITE DO 613 J=1,3 XP(J)=XCORD(J,I) XV(J)=XP(J)-XO(J) 613 CONTINUE SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3) SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))*XN(3) XPROJ(3,I)=-SN DO 614 J=1,3 XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J) 614 CONTINUE C+PP option DIRE IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN C+PP XPROJ(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3) XPROJ(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3) C+PP option DIRE ELSE XPROJ(1,I)=(XMINT+XMAXT)/2 XPROJ(2,I)=(YMINT+YMAXT)/2 ENDIF C+PP 612 CONTINUE IF (LABCO2.EQ.0) GOTO 618 IF (LABCO2(3,ICLE).EQ.0) GOTO 618 KABCO2=LABCO2(1,ICLE) NVEC=KABCO2(/2) SEGINI KXPRO2 LABCO2(2,ICLE)=KXPRO2 IF (NVEC.EQ.0) GOTO 618 DO 8 IVEC=1,NVEC SXCOR2=KABCO2(1,IVEC) ICOR2=KABCO2(2,IVEC) * ajout SG 2016/07/16 apparemment le ITE de XCOR2 n'est pas * forcement celui de XCORD en cas de coupe ITE=XCOR2(/2) SEGINI XPRO2 KXPRO2(IVEC)=XPRO2 DO 617 I=1,ITE IF (ICOR2(I).EQ.0) GOTO 617 DO 615 J=1,3 XP(J)=XCOR2(J,I) XV(J)=XP(J)-XO(J) 615 CONTINUE SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3) SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3)) $ *XN(3) XPRO2(3,I)=-SN DO 616 J=1,3 XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J) 616 CONTINUE * write(6,*) ' cadrcl on passe la 353 ',i XPRO2(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3) XPRO2(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3) 617 CONTINUE ** SEGSUP SXCOR2 8 CONTINUE 618 CONTINUE ** SEGSUP SXCORD * write(6,*) '6 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt RETURN 1100 CONTINUE DO 1110 I=1,IDIM DO 1110 J=1,ITE XPROJ(I,J)=XCORD(I,J) 1110 CONTINUE IF (LABCO2.EQ.0) GOTO 1111 IF (LABCO2(3,ICLE).EQ.0) GOTO 1111 KABCO2=LABCO2(1,ICLE) NVEC=KABCO2(/2) SEGINI KXPRO2 LABCO2(2,ICLE)=KXPRO2 IF (NVEC.EQ.0) GOTO 1111 DO 9 IVEC=1,NVEC SXCOR2=KABCO2(1,IVEC) ICOR2=KABCO2(2,IVEC) segact sxcor2 ITE=XCOR2(/2) SEGINI XPRO2 KXPRO2(IVEC)=XPRO2 * write(6,*) ' cadrcl on passe aussi la 381',j DO 1113 J=1,ITE IF (ICOR2(J).EQ.0) GOTO 1113 DO 1112 I=1,IDIM XPRO2(I,J)=XCOR2(I,J) 1112 CONTINUE * write(6,*) j,xcor2,(xcor2(i,j),i=1,idim) 1113 CONTINUE ** SEGSUP SXCOR2 9 CONTINUE 1111 CONTINUE ** SEGSUP SXCORD * axez pour tourner correctement avec opengl axez(1)=0 axez(2)=uj(3) axez(3)=sqrt(1-uj(3)**2) if (xn(3).lt.0) axez(3)=-axez(3) * write (6,*) ' axez ',axez(1),axez(2),axez(3) * write(6,*) '7 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt RETURN END