C PROOBJ SOURCE SP204843 25/03/14 21:15:08 12201 C CE SOUS-PROGRAMME PREPARE LA PROJECTION D'UN OBJET SUR UNE C SURFACE : C PLAN P1 P2 P3 C SPHE C P C CYLI C1 C2 P C CONI C C1 P C TORI C A Cs P C NOUVELLE POSSIBILITE EN 2D NOVEMBRE 1986 C DROI P1 P2 C CERCLE P1 C C JANVIER 1987 INTRODUCTION DE LA PROJECTION CONIQUE C C SEPTEMBRE 1996 PROJECTION D'UN OBJET MAILLAGE SUR UN OBJET C MAILLAGE. SUBROUTINE PROOBJ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCREEL -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC CCGEOME -INC CCTOURN CHARACTER*4 MCLE(7),MTYPP(3) DATA MTYPP/'DIRE','CONI','CYLI'/ DATA MCLE/'PLAN','SPHE','CYLI','CONI','TORI','DROI','CERC'/ SEGACT MCOORD MELEME=0 CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU) IF (IRETOU.EQ.1) GOTO 1 C C'EST UN POINT QU'ON PROJETTE CALL LIROBJ('POINT ',IPOIN,1,IRETOU) IF (IERR.NE.0) RETURN 1 CONTINUE * CONIQUE OU CYLINDRIQUE CALL LIRMOT(MTYPP,3,IYYT,0) CALL LIROBJ('POINT ',IP1,1,IRETOU) IF (IERR.NE.0) RETURN IF (IYYT.EQ.3) THEN IP2 = 0 CALL LIROBJ('POINT ',IP2,1,IRETOU) IF (IERR.NE.0) RETURN ENDIF CALL LIROBJ('MAILLAGE',IPPP,0,IRETO) IF (IYYT.EQ.1) THEN IF(IRETO.NE.0) THEN if (meleme.ne.0) then CALL PROOB1(MELEME,IPPP,IP1) else call erreur(821) return endif RETURN ENDIF ENDIF IF (IYYT.EQ.2) THEN IF(IRETO.NE.0) THEN CALL PROOB4(MELEME,IPPP,IP1) RETURN ENDIF ENDIF SEGACT MCOORD IREF=(IDIM+1)*(IP1-1) XVEC=XCOOR(IREF+1) YVEC=XCOOR(IREF+2) ZVEC=XCOOR(IREF+3) IF (IDIM.EQ.2) ZVEC=0. IF (IYYT.EQ.3) THEN IF (IDIM.EQ.2) THEN INTERR(1)=2 CALL ERREUR(709) RETURN ENDIF C write(6,*) 'Option CYL2' XP1 = XVEC YP1 = YVEC ZP1 = ZVEC IREF=(IDIM+1)*(IP2-1) XP2=XCOOR(IREF+1) YP2=XCOOR(IREF+2) ZP2=XCOOR(IREF+3) IF (IDIM.EQ.2) ZP2=0. XVEC=XP2-XP1 YVEC=YP2-YP1 ZVEC=ZP2-ZP1 SVEC=XVEC**2+YVEC**2+ZVEC**2 UMAX=max(ABS(XP1),ABS(XP2)) VMAX=max(ABS(YP1),ABS(YP2)) WMAX=max(ABS(ZP1),ABS(ZP2)) IF(VMAX.GT.UMAX) UMAX=VMAX IF(WMAX.GT.UMAX) UMAX=WMAX IF (SVEC.LT.(UMAX*XZPREC)) THEN CALL ERREUR(237) RETURN ENDIF SVEC=SQRT(SVEC) XVEC=XVEC/SVEC YVEC=YVEC/SVEC ZVEC=ZVEC/SVEC C write(6,*) 'XP1,YP1,ZP1',XP1,YP1,ZP1 C write(6,*) 'XP2,YP2,ZP2',XP2,YP2,ZP2 C write(6,*) 'XVEC,YVEC,ZVEC,SVEC',XVEC,YVEC,ZVEC,SVEC ENDIF IALIR=0 IF (IYYT.EQ.0) IALIR=1 CALL LIRMOT(MCLE,7,IVAL,IALIR) IF (IVAL.EQ.0) THEN IVAL=IYYT+2 IYYT=0 ENDIF IF (IERR.NE.0) RETURN IF (IDIM.EQ.2.AND.IVAL.LE.5) CALL ERREUR(34) IF (IDIM.EQ.3.AND.IVAL.GE.6) CALL ERREUR(34) IF (IERR.NE.0) RETURN ICLE=IVAL+7 GOTO (10,20,30,40,50,60,70),IVAL 10 CONTINUE C PLAN ON LIT TROIS POINTS CALL LIROBJ('POINT ',IP1,1,IRETOU) CALL LIROBJ('POINT ',IP2,1,IRETOU) CALL LIROBJ('POINT ',IP3,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IP1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IDIM+1)*(IP2-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 IREF=(IDIM+1)*(IP3-1) XV3=XCOOR(IREF+1)-XPT1 YV3=XCOOR(IREF+2)-YPT1 ZV3=XCOOR(IREF+3)-ZPT1 C ON GARDE LE VECTEUR NORMAL NORMALISE XV1=YV2*ZV3-ZV2*YV3 YV1=ZV2*XV3-XV2*ZV3 ZV1=XV2*YV3-YV2*XV3 SV1=SQRT(XV1**2+YV1**2+ZV1**2) IF (SV1.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 GOTO 100 20 CONTINUE C SPHERE ON LIT LE CENTRE ET UN POINT CALL LIROBJ('POINT ',IPCEN,1,IRETOU) CALL LIROBJ('POINT ',IP1,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IPCEN-1)*(IDIM+1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IP1-1)*(IDIM+1) XV=XCOOR(IREF+1)-XPT1 YV=XCOOR(IREF+2)-YPT1 ZV=XCOOR(IREF+3)-ZPT1 ANGLE=SQRT(XV**2+YV**2+ZV**2) IF (ANGLE.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN GOTO 100 30 CONTINUE C CYLINDRE ON LIT DEUX POINTS DE L'AXE ET UN POINT DU CYLINDRE CALL LIROBJ('POINT ',IPOIN1,1,IRETOU) CALL LIROBJ('POINT ',IPOIN2,1,IRETOU) CALL LIROBJ('POINT ',IP1,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IPOIN1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IDIM+1)*(IPOIN2-1) XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 S=SQRT(XV1**2+YV1**2+ZV1**2) IF (S.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV1=XV1/S YV1=YV1/S ZV1=ZV1/S IREF=(IDIM+1)*(IP1-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 XV3=YV1*ZV2-ZV1*YV2 YV3=ZV1*XV2-XV1*ZV2 ZV3=XV1*YV2-YV1*XV2 ANGLE=SQRT(XV3**2+YV3**2+ZV3**2) C RAYON DU CYLINDRE IF (ANGLE.EQ.0) CALL ERREUR(21) IF (IERR.NE.0) RETURN GOTO 100 40 CONTINUE C CONE ON LIT LE SOMMET UN POINT DE L'AXE ET UN POINT DU CONE CALL LIROBJ('POINT ',IPT1,1,IRETOU) CALL LIROBJ('POINT ',IP1,1,IRETOU) CALL LIROBJ('POINT ',IP2,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IPT1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IDIM+1)*(IP1-1) XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 SV1=SQRT(XV1**2+YV1**2+ZV1**2) IF (SV1.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 IREF=(IDIM+1)*(IP2-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 SV2=SQRT(XV2**2+YV2**2+ZV2**2) IF (SV2.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV2=XV2/SV2 YV2=YV2/SV2 ZV2=ZV2/SV2 ANGLE=(XV1*XV2+YV1*YV2+ZV1*ZV2)**2 IF (IIMPI.NE.0) WRITE (IOIMP,1012) ANGLE 1012 FORMAT(' COS **2 DE ANGLE AU SOMMET ',G12.5) GOTO 100 50 CONTINUE C TORE ON LIT LE CENTRE UN POINT DE L'AXE UN CENTRE SECONDAIRE ET C UN POINT CALL LIROBJ('POINT ',IPT1,1,IRETOU) CALL LIROBJ('POINT ',IPAX,1,IRETOU) CALL LIROBJ('POINT ',IPCS,1,IRETOU) CALL LIROBJ('POINT ',IP1 ,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IPT1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IDIM+1)*(IPAX-1) XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 SV1=XV1**2+YV1**2+ZV1**2 IF (SV1.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SV1=SQRT(SV1) XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 IREF=(IDIM+1)*(IPCS-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 SCA=XV2*XV1+YV2*YV1+ZV2*ZV1 XPT1=XPT1+SCA*XV1 YPT1=YPT1+SCA*YV1 ZPT1=ZPT1+SCA*ZV1 XV2=XV2-SCA*XV1 YV2=YV2-SCA*YV1 ZV2=ZV2-SCA*ZV1 GR2=XV2**2+YV2**2+ZV2**2 IF (GR2.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IP1-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 SCA=XV2*XV1+YV2*YV1+ZV2*ZV1 XC=XV2-SCA*XV1 YC=YV2-SCA*YV1 ZC=ZV2-SCA*ZV1 SC=XC**2+YC**2+ZC**2 IF (SC.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN RAP=SQRT(GR2/SC) XC=XC*RAP YC=YC*RAP ZC=ZC*RAP PR2=(XV2-XC)**2+(YV2-YC)**2+(ZV2-ZC)**2 IF (PR2.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV2=PR2 YV2=GR2 IF (IIMPI.NE.0) WRITE (IOIMP,1015) XV2,YV2,XPT1,YPT1,ZPT1, # XV1,YV1,ZV1 1015 FORMAT(' CARRE PETIT RAYON ',G12.5, # ' CARRE GRAND RAYON ',G12.5,/,' CENTRE ',3G13.5,' AXE ', # 3G12.5) GOTO 100 60 CONTINUE C DROITE ON LIT DEUX POINTS CALL LIROBJ('POINT ',IP1,1,IRETOU) CALL LIROBJ('POINT ',IP2,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IDIM+1)*(IP1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=0. IREF=(IDIM+1)*(IP2-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=0. XV3=0. YV3=0. ZV3=1. C ON GARDE LE VECTEUR NORMAL NORMALISE XV1=YV2*ZV3-ZV2*YV3 YV1=ZV2*XV3-XV2*ZV3 ZV1=XV2*YV3-YV2*XV3 SV1=SQRT(XV1**2+YV1**2+ZV1**2) IF (SV1.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 GOTO 100 70 CONTINUE C CERCLE ON LIT LE CENTRE ET UN POINT CALL LIROBJ('POINT ',IPCEN,1,IRETOU) CALL LIROBJ('POINT ',IP1,1,IRETOU) IF (IERR.NE.0) RETURN IREF=(IPCEN-1)*(IDIM+1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=0. IREF=(IP1-1)*(IDIM+1) XV=XCOOR(IREF+1)-XPT1 YV=XCOOR(IREF+2)-YPT1 ZV=0. ANGLE=SQRT(XV**2+YV**2+ZV**2) IF (ANGLE.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN GOTO 100 100 CONTINUE IF (MELEME.EQ.0) GOTO 101 IF (IYYT.EQ.2) THEN IVAL=IVAL+10000 ICLE=ICLE+10000 ENDIF IF (IYYT.EQ.3) THEN IVAL=IVAL-10000 ICLE=ICLE-10000 ENDIF CALL INTOPE(MELEME) RETURN 101 CONTINUE IREF=(IDIM+1)*(IPOIN-1) XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) TPOIN=XCOOR(IREF+4) IF (IDIM.EQ.2) ZPOIN=0 IF (IDIM.EQ.2) TPOIN=XCOOR(IREF+3) IF (IYYT.EQ.2) THEN XVEC=XVEC-XPOIN YVEC=YVEC-YPOIN ZVEC=ZVEC-ZPOIN IF (IDIM.EQ.2) ZVEC=0. ENDIF IF (IYYT.EQ.3) THEN CALL ERREUR(251) RETURN ENDIF SVEC=SQRT(XVEC**2+YVEC**2+ZVEC**2) IF (SVEC.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN XVEC=XVEC/SVEC YVEC=YVEC/SVEC ZVEC=ZVEC/SVEC GOTO (110,120,130,140,150,160,170),IVAL 110 CONTINUE XV=XPOIN-XPT1 YV=YPOIN-YPT1 ZV=ZPOIN-ZPT1 DENUM=XV*XV1+YV*YV1+ZV*ZV1 DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 IF (DENOM.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN RAP=-DENUM/DENOM XPOIN=XPOIN+RAP*XVEC YPOIN=YPOIN+RAP*YVEC ZPOIN=ZPOIN+RAP*ZVEC GOTO 200 120 CONTINUE XV=XPT1-XPOIN YV=YPT1-YPOIN ZV=ZPT1-ZPOIN SCA=XVEC*XV+YVEC*YV+ZVEC*ZV XV=XVEC*SCA YV=YVEC*SCA ZV=ZVEC*SCA S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2 IF (S2.GT.ANGLE**2) CALL ERREUR(21) IF (IERR.NE.0) RETURN C=SQRT(ANGLE**2-S2) IF (SCA.LT.0.) C=-C XPOIN=XPOIN+XV-C*XVEC YPOIN=YPOIN+YV-C*YVEC ZPOIN=ZPOIN+ZV-C*ZVEC GOTO 200 130 CONTINUE C V2 = VEC VECT V1 XV2=YVEC*ZV1-ZVEC*YV1 YV2=ZVEC*XV1-XVEC*ZV1 ZV2=XVEC*YV1-YVEC*XV1 S2V2=XV2**2+YV2**2+ZV2**2 IF (S2V2.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN C2=(XVEC*XV1+YVEC*YV1+ZVEC*ZV1)**2 IF (C2.EQ.1.) CALL ERREUR(21) IF (IERR.NE.0) RETURN C V3= POIN PT1 XV3=XPT1-XPOIN YV3=YPT1-YPOIN ZV3=ZPT1-ZPOIN XV4=YV3*ZV1-ZV3*YV1 YV4=ZV3*XV1-XV3*ZV1 ZV4=XV3*YV1-YV3*XV1 DNUM=(XV4*XVEC+YV4*YVEC+ZV4*ZVEC)**2 DIS2=DNUM/S2V2 IF (IIMPI.NE.0) WRITE (IOIMP,1000) DIS2,ANGLE 1000 FORMAT (' DISTANCE AU CARRE DES AXES ',G12.5, # 'RAYON DU CYLINDRE ',G12.5) IF (DIS2.GT.ANGLE**2) CALL ERREUR(21) IF (IERR.NE.0) RETURN DMU=SQRT((ANGLE**2-DIS2)/(1.-C2)) DNUM=XV2*XV4+YV2*YV4+ZV2*ZV4 DLA=DNUM/S2V2 DMU=SIGN(DMU,DLA) IF (IIMPI.NE.0) WRITE (IOIMP,1001) DLA,DMU 1001 FORMAT(' LAMBDA,MU ',2G13.5) XPOIN=XPOIN+XVEC*(DLA-DMU) YPOIN=YPOIN+YVEC*(DLA-DMU) ZPOIN=ZPOIN+ZVEC*(DLA-DMU) GOTO 200 140 CONTINUE XV2=XPOIN-XPT1 YV2=YPOIN-YPT1 ZV2=ZPOIN-ZPT1 VECV1=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 VEC2=XVEC**2+YVEC**2+ZVEC**2 V2V1=XV2*XV1+YV2*YV1+ZV2*ZV1 VECV2=XVEC*XV2+YVEC*YV2+ZVEC*ZV2 V22=XV2**2+YV2**2+ZV2**2 A=VECV1**2-ANGLE*VEC2 B=2*(V2V1*VECV1-ANGLE*VECV2) C=V2V1**2-ANGLE*V22 DELTA=B**2-4*A*C IF (DELTA.LT.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN DEL=SQRT(DELTA) X1=(-B+DEL)/(2*A) X2=(-B-DEL)/(2*A) X=X2 IF (ABS(X1).LT.ABS(X2)) X=X1 XPOIN=XPOIN+X*XVEC YPOIN=YPOIN+X*YVEC ZPOIN=ZPOIN+X*ZVEC GOTO 200 150 CONTINUE PR2=XV2 GR2=YV2 XOP=XPOIN-XPT1 YOP=YPOIN-YPT1 ZOP=ZPOIN-ZPT1 OPV=XOP*XVEC+YOP*YVEC+ZOP*ZVEC R2=XOP**2+YOP**2+ZOP**2-GR2-PR2 VA=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 QR2VA2=4*GR2*VA**2 OPA=XOP*XV1+YOP*YV1+ZOP*ZV1 HR2PV=8*GR2*OPA*VA R=4*GR2*OPA**2-4*PR2*GR2 RLMD=0 C RESOLUTION DE L'EQUATION DU 4EME DEGRE PAR ITERATION DO 151 ITER=1,25 EXP1=RLMD*(RLMD+2*OPV)+R2 FDLM=EXP1**2+QR2VA2*RLMD**2+HR2PV*RLMD+R FPDLM=4*EXP1*(RLMD+OPV)+QR2VA2*2*RLMD+HR2PV IF (FPDLM.EQ.0.) CALL ERREUR(40) IF (IERR.NE.0) RETURN CORR=FDLM/FPDLM IF (IIMPI.NE.0) WRITE(IOIMP,1016) ITER,RLMD,CORR 1016 FORMAT(' ITER ',I6,' LAMBDA ',G12.5,' CORR ',G12.5) RLMD=RLMD-CORR IF (RLMD.EQ.0.) GOTO 151 IF (ABS(CORR/RLMD).LT.1E-5) GOTO 152 151 CONTINUE CALL ERREUR(40) RETURN 152 CONTINUE XPOIN=XPOIN+RLMD*XVEC YPOIN=YPOIN+RLMD*YVEC ZPOIN=ZPOIN+RLMD*ZVEC GOTO 200 160 CONTINUE XV=XPOIN-XPT1 YV=YPOIN-YPT1 ZV=ZPOIN-ZPT1 DENUM=XV*XV1+YV*YV1+ZV*ZV1 DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 IF (DENOM.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN RAP=-DENUM/DENOM XPOIN=XPOIN+RAP*XVEC YPOIN=YPOIN+RAP*YVEC ZPOIN=ZPOIN+RAP*ZVEC GOTO 200 170 CONTINUE XV=XPT1-XPOIN YV=YPT1-YPOIN ZV=ZPT1-ZPOIN SCA=XVEC*XV+YVEC*YV+ZVEC*ZV XV=XVEC*SCA YV=YVEC*SCA ZV=ZVEC*SCA S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2 IF (S2.GT.ANGLE**2) CALL ERREUR(21) IF (IERR.NE.0) RETURN C=SQRT(ANGLE**2-S2) IF (SCA.LT.0.) C=-C XPOIN=XPOIN+XV-C*XVEC YPOIN=YPOIN+YV-C*YVEC ZPOIN=ZPOIN+ZV-C*ZVEC GOTO 200 200 CONTINUE SEGACT MCOORD*mod NBPTS=nbpts+1 SEGADJ MCOORD XCOOR((NBPTS-1)*(IDIM+1)+1)=XPOIN XCOOR((NBPTS-1)*(IDIM+1)+2)=YPOIN IF (IDIM.EQ.3) XCOOR((NBPTS-1)*(IDIM+1)+3)=ZPOIN XCOOR(NBPTS*(IDIM+1))=TPOIN CALL ECROBJ('POINT ',NBPTS) RETURN END