cadrcl
C CADRCL SOURCE PV090527 24/06/12 21:15:01 11940 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 # ICLE,XMINT,YMINT,XMAXT,YMAXT,zmint,zmaxt,cgrav,diloc,LDIRE, > axez) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMCOORD dimension cgrav(*),axez(*) C+PP option DIRE dimension diloc(*) LOGICAL LDIRE C+PP DIMENSION XO(3),XP(3),XN(3),XG(3),XV(3),UI(3),UJ(3) 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 SEGMENT XPROJ(3,ITE) SEGMENT XPRO2(3,ITE) C+PP option DIRE REAL*8 SSN *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 ZMINT=0. ZMAXT=0. DO ICORD=1,NCORD SXCORD=KABCOR(ICORD) * write(6,*) '1 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) DO 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)) ENDDO ENDDO IF (ICLE.NE.0) GOTO 1000 IF (LABCO2.EQ.0) GOTO 60 DO 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 IVEC=1,NVEC SXCOR2=KABCO2(1,IVEC) ICOR2=KABCO2(2,IVEC) ITE=XCOR2(/2) DO 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)) ENDDO 18 CONTINUE ENDDO 61 CONTINUE ENDDO 60 CONTINUE * pour eviter des calculs sur des quantites non significatives, * on agrandit un peu le cadre XDIFF=XMAXT-XMINT YDIFF=YMAXT-YMINT c* ZDIFF=ZMAXT-ZMINT ZDIFF=0. XSOMM=XMAXT+XMINT YSOMM=YMAXT+YMINT c* ZSOMM=ZMAXT+ZMINT ZSOMM=0. 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 C* Cas DIME = 3 : 50 CONTINUE SEGACT,MCOORD IREF=(IOEIL-1)*4 XO(1)=XCOOR(IREF+1) XO(2)=XCOOR(IREF+2) XO(3)=XCOOR(IREF+3) C+PP option DIRE IF (LDIRE)THEN DO J=1,3 XG(J)=cgrav(J) XN(J)=XO(J)-XG(J) ENDDO C+PP option DIRE C+PP ELSE C POINT MOYEN XG(1)=0.D0 XG(2)=0.D0 XG(3)=0.D0 NPTOT=0 DO ICORD=1,NCORD SXCORD=KABCOR(ICORD) * write(6,*) '3 cadrcl icord sxcord',icord,sxcord ITE=XCORD(/2) NPTOT=NPTOT+ITE DO I=1,ITE XG(1)=XG(1)+XCORD(1,I) XG(2)=XG(2)+XCORD(2,I) XG(3)=XG(3)+XCORD(3,I) ENDDO ENDDO DO J=1,3 XG(J)=XG(J)/NPTOT XN(J)=XO(J)-XG(J) cgrav(j)=XG(j) ENDDO 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.) THEN RETURN ENDIF C+PP option DIRE SSN=SN C+PP XN(1)=XN(1)/SN XN(2)=XN(2)/SN XN(3)=XN(3)/SN C REPERE LOCAL SUR LE PLAN C+PP option DIRE IF (LDIRE)THEN UJ(1)=diloc(1) UJ(2)=diloc(2) UJ(3)=diloc(3) 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 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.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) if (icor2.le.0) goto 6 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 J=1,ITE DO I=1,IDIM XPROJ(I,J)=XCORD(I,J) ENDDO ENDDO 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 I=1,IDIM XPRO2(I,J)=XCOR2(I,J) ENDDO * 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 c RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales