C ADDITE SOURCE CB215821 25/04/08 21:15:02 12227 C Sous-programme elementaire pour effectuer la translation, la rotation C ou l'affinite d'un objet C 10/2003 : Prise en compte cas IDIM=1. SUBROUTINE ADDITE(X,IPT1,IPT2,ICPR,IARG) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMCOORD -INC CCGEOME -INC CCTOURN DIMENSION X(*) SEGMENT ICPR(nbpts) idimp1=IDIM+1 XVECS=XVEC YVECS=YVEC ZVECS=ZVEC COSA=COS(ANGLE) SINA=SIN(ANGLE) ANG2=ANGLE*ANGLE C Initialisation du maillage IPT2 transforme de IPT1 NBSOUS=0 IF (IARG.EQ.0) THEN NBREF=0 ELSE NBREF=IPT1.LISREF(/1) ENDIF NBNN=IPT1.NUM(/1) NBELEM=IPT1.NUM(/2) SEGINI IPT2 IPT2.ITYPEL=IPT1.ITYPEL DO i=1,NBELEM IPT2.ICOLOR(i)=IPT1.ICOLOR(i) ENDDO C ICLE / KCLE : indique la transformation a realiser C KCLE = 1 : translation (operateurs PLUS,MOIN,DEDU 'TRAN') C KCLE = 2 : rotation (operateurs TOUR,DEDU 'ROTA') C KCLE = 3 : homothetie (operateur HOMO) C KCLE = 4 : affinite (operateur AFFI) C KCLE = 5 : symetrie point (operateur SYME 'POINT') C KCLE = 6 : symetrie droite (operateur SYME 'DROIT') C KCLE = 7 : symetrie plan (operateur SYME 'PLAN') C KCLE = 8 : projection plan (operateur PROJ 'PLAN') C KCLE = 9 : projection sphere (operateur PROJ 'SPHE') C KCLE = 10 : projection cylindre (operateur PROJ 'CYLI') C KCLE = 11 : projection conique (operateur PROJ 'CONI') C KCLE = 12 : projection torique (operateur PROJ 'TORI') C KCLE = 13 : projection droite (operateur PROJ 'DROI') C KCLE = 14 : projection cercle (operateur PROJ 'CERC') KCLE=ICLE IF (ICLE.GE.10000) KCLE=ICLE-10000 IF (ICLE.LT.0) KCLE=ICLE+10000 C Reservation de la place dans XCOOR segdes mcoord SEGACT MCOORD*mod NPACT=NBPTS NXCOUR=NPACT*idimp1 NBPTS=NPACT+NBELEM*NBNN SEGADJ MCOORD C Boucle sur les noeuds du maillage C Creation du noeud de IPT2 image de IPT1.NUM(i,j) et stocke dans ICPR DO 10 j=1,NBELEM DO 11 i=1,NBNN C L'image de IPT1.NUM(i,j) a-t-elle deja ete creee ? IF (ICPR(IPT1.NUM(i,j)).EQ.0) THEN IREF=(IPT1.NUM(i,j)-1)*idimp1 GOTO (51,52,53,54,55,56,57,58,59,60,61,62,63,64),KCLE C Translation de vecteur X : 51 DO k=1,IDIM XCOOR(NXCOUR+k)=X(k)+XCOOR(IREF+k) ENDDO XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Rotation (2D/3D) : 52 XD=XCOOR(IREF+1)-XPT1 YD=XCOOR(IREF+2)-YPT1 ZD=0. IF (IDIM.GE.3) ZD=XCOOR(IREF+3)-ZPT1 XE=XD*XV1+YD*YV1+ZD*ZV1 YE=XD*XV2+YD*YV2+ZD*ZV2 ZE=XD*XVEC+YD*YVEC+ZD*ZVEC XD=XE*COSA-YE*SINA YD=XE*SINA+YE*COSA ZD=ZE XCOOR(NXCOUR+1)=XD*XV1+YD*XV2+ZD*XVEC+XPT1 XCOOR(NXCOUR+2)=XD*YV1+YD*YV2+ZD*YVEC+YPT1 IF (IDIM.GE.3) XCOOR(NXCOUR+3)=XD*ZV1+YD*ZV2+ZD*ZVEC+ZPT1 XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Homothetie : 53 XCOOR(NXCOUR+1)=(XCOOR(IREF+1)-XPT1)*ANGLE+XPT1 IF (IDIM.GE.2) XCOOR(NXCOUR+2)=(XCOOR(IREF+2)-YPT1)*ANGLE+YPT1 IF (IDIM.GE.3) XCOOR(NXCOUR+3)=(XCOOR(IREF+3)-ZPT1)*ANGLE+ZPT1 XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1)*ANGLE GOTO 70 C Affinite (2D/3D) : 54 XD=XCOOR(IREF+1)-XPT1 YD=XCOOR(IREF+2)-YPT1 ZD=0. IF (IDIM.GE.3) ZD=XCOOR(IREF+3)-ZPT1 XE=XD*XV1+YD*YV1+ZD*ZV1 YE=XD*XV2+YD*YV2+ZD*ZV2 ZE=(XD*XVEC+YD*YVEC+ZD*ZVEC)*ANGLE XCOOR(NXCOUR+1)=XE*XV1+YE*XV2+ZE*XVEC+XPT1 XCOOR(NXCOUR+2)=XE*YV1+YE*YV2+ZE*YVEC+YPT1 IF (IDIM.GE.3) XCOOR(NXCOUR+3)=XE*ZV1+YE*ZV2+ZE*ZVEC+ZPT1 XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Symetrie par rapport a un point : 55 XCOOR(NXCOUR+1)=2*XPT1-XCOOR(IREF+1) IF (IDIM.GE.2) XCOOR(NXCOUR+2)=2*YPT1-XCOOR(IREF+2) IF (IDIM.GE.3) XCOOR(NXCOUR+3)=2*ZPT1-XCOOR(IREF+3) XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Symetrie par rapport a une droite (2D/3D) : 56 XD=XCOOR(IREF+1)-XPT1 YD=XCOOR(IREF+2)-YPT1 ZD=0. IF (IDIM.GE.3) ZD=XCOOR(IREF+3)-ZPT1 PVEC=2*(XD*XVEC+YD*YVEC+ZD*ZVEC) XCOOR(NXCOUR+1)=XPT1+XVEC*PVEC-XD XCOOR(NXCOUR+2)=YPT1+YVEC*PVEC-YD IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPT1+ZVEC*PVEC-ZD XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Symetrie par rapport a un plan (3D) : 57 XD=XCOOR(IREF+1)-XPT1 YD=XCOOR(IREF+2)-YPT1 ZD=XCOOR(IREF+3)-ZPT1 PVEC=2*(XD*XVEC+YD*YVEC+ZD*ZVEC) XCOOR(NXCOUR+1)=XPT1+XD-XVEC*PVEC XCOOR(NXCOUR+2)=YPT1+YD-YVEC*PVEC XCOOR(NXCOUR+3)=ZPT1+ZD-ZVEC*PVEC XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur un plan : 58 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=ZVEC-ZPOIN ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=ZPX-ZPOIN ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS DENUM=(XPOIN-XPT1)*XV1+(YPOIN-YPT1)*YV1+(ZPOIN-ZPT1)*ZV1 DENOM=XVECS*XV1+YVECS*YV1+ZVECS*ZV1 IF (DENOM.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN RAP=-DENUM/DENOM XCOOR(NXCOUR+1)=XPOIN+RAP*XVECS XCOOR(NXCOUR+2)=YPOIN+RAP*YVECS IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPOIN+RAP*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur une sphere (3D) : 59 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=ZVEC-ZPOIN ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=ZPX-ZPOIN ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS SCA=XVECS*(XPT1-XPOIN)+YVECS*(YPT1-YPOIN)+ZVECS*(ZPT1-ZPOIN) XV=XVECS*SCA YV=YVECS*SCA ZV=ZVECS*SCA S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2 IF (S2.GT.ANG2) CALL ERREUR(21) IF (IERR.NE.0) RETURN C=SQRT(ANG2-S2) IF (SCA.LT.0.) C=-C XCOOR(NXCOUR+1)=XPOIN+XV-C*XVECS XCOOR(NXCOUR+2)=YPOIN+YV-C*YVECS IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPOIN+ZV-C*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur un cylindre (3D) : 60 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=ZVEC-ZPOIN ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=ZPX-ZPOIN ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS XV2=YVECS*ZV1-ZVECS*YV1 YV2=ZVECS*XV1-XVECS*ZV1 ZV2=XVECS*YV1-YVECS*XV1 S2V2=XV2*XV2+YV2*YV2+ZV2*ZV2 IF (S2V2.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN C2=(XVECS*XV1+YVECS*YV1+ZVECS*ZV1)**2 IF (C2.EQ.1.) CALL ERREUR(21) IF (IERR.NE.0) RETURN 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*XVECS+YV4*YVECS+ZV4*ZVECS)**2 DIS2=DNUM/S2V2 IF (IIMPI.EQ.1) WRITE (IOIMP,1000) DIS2,ANGLE 1000 FORMAT (' DISTANCE AU CARRE DES AXES ',G12.5, . 'RAYON DU CYLINDRE ',G12.5) IF (DIS2.GT.ANG2) CALL ERREUR(21) IF (IERR.NE.0) RETURN DMU=SQRT((ANG2-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) XCOOR(NXCOUR+1)=XPOIN+XVECS*(DLA-DMU) XCOOR(NXCOUR+2)=YPOIN+YVECS*(DLA-DMU) IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPOIN+ZVECS*(DLA-DMU) XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur un cone (3D) : 61 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=ZVEC-ZPOIN ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=ZPX-ZPOIN ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS XV2=XPOIN-XPT1 YV2=YPOIN-YPT1 ZV2=ZPOIN-ZPT1 VECV1=XVECS*XV1+YVECS*YV1+ZVECS*ZV1 VEC2=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS V2V1=XV2*XV1+YV2*YV1+ZV2*ZV1 VECV2=XVECS*XV2+YVECS*YV2+ZVECS*ZV2 V22=XV2*XV2+YV2*YV2+ZV2*ZV2 A=VECV1*VECV1-ANGLE*VEC2 B=2*(V2V1*VECV1-ANGLE*VECV2) C=V2V1*V2V1-ANGLE*V22 DELTA=B*B-4*A*C IF (DELTA.LT.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN DEL=SQRT(DELTA) X1=(-B+DEL)/(2*A) XX=(-B-DEL)/(2*A) IF (ABS(X1).LT.ABS(XX)) XX=X1 XCOOR(NXCOUR+1)=XPOIN+XX*XVECS XCOOR(NXCOUR+2)=YPOIN+XX*YVECS XCOOR(NXCOUR+3)=ZPOIN+XX*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur un tore (3D) : 62 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=XCOOR(IREF+3) IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=ZVEC-ZPOIN ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=ZPX-ZPOIN ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS PR2=XV2 GR2=YV2 XOP=XPOIN-XPT1 YOP=YPOIN-YPT1 ZOP=ZPOIN-ZPT1 OPV=XOP*XVECS+YOP*YVECS+ZOP*ZVECS R2=XOP*XOP+YOP*YOP+ZOP*ZOP-GR2-PR2 VA=XVECS*XV1+YVECS*YV1+ZVECS*ZV1 QR2VA2=4*GR2*VA*VA OPA=XOP*XV1+YOP*YV1+ZOP*ZV1 HR2PV=8*GR2*OPA*VA R=4*GR2*OPA*OPA-4*PR2*GR2 RLMD=0 C Resolution iterative de l'equation du 4eme degre DO ITER=1,25 EXP1=RLMD*(RLMD+2*OPV)+R2 FDLM=EXP1*EXP1+QR2VA2*RLMD*RLMD+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.NE.0.) THEN IF (ABS(CORR/RLMD).LT.1E-5) GOTO 66 ENDIF ENDDO CALL ERREUR(40) RETURN 66 XCOOR(NXCOUR+1)=XPOIN+RLMD*XVECS XCOOR(NXCOUR+2)=YPOIN+RLMD*YVECS XCOOR(NXCOUR+3)=ZPOIN+RLMD*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur une droite (2D/3D) : 63 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=0. IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=0. ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=0. ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS DENUM=(XPOIN-XPT1)*XV1+(YPOIN-YPT1)*YV1 DENOM=XVECS*XV1+YVECS*YV1+ZVECS*ZV1 IF (DENOM.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN RAP=-DENUM/DENOM XCOOR(NXCOUR+1)=XPOIN+RAP*XVECS XCOOR(NXCOUR+2)=YPOIN+RAP*YVECS IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPOIN+RAP*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Projection sur un cercle (2D/3D) : 64 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=0. IF (ICLE.GE.10000) THEN XVECS=XVEC-XPOIN YVECS=YVEC-YPOIN ZVECS=0. ENDIF IF (ICLE.LE.0) THEN XPX=XP1-XPOIN YPX=YP1-YPOIN ZPX=ZP1-ZPOIN SPS=XPX*XVEC+YPX*YVEC+ZPX*ZVEC XPX=XP1-SPS*XVEC YPX=YP1-SPS*YVEC ZPX=ZP1-SPS*ZVEC XVECS=XPX-XPOIN YVECS=YPX-YPOIN ZVECS=0. ENDIF SVECS=XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS IF (SVECS.EQ.0.) CALL ERREUR(21) IF (IERR.NE.0) RETURN SVECS=SQRT(SVECS) XVECS=XVECS/SVECS YVECS=YVECS/SVECS ZVECS=ZVECS/SVECS SCA=XVECS*(XPT1-XPOIN)+YVECS*(YPT1-YPOIN) XV=XVECS*SCA YV=YVECS*SCA ZV=ZVECS*SCA S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2 IF (S2.GT.ANG2) CALL ERREUR(21) IF (IERR.NE.0) RETURN C=SQRT(ANG2-S2) IF (SCA.LT.0.) C=-C XCOOR(NXCOUR+1)=XPOIN+XV-C*XVECS XCOOR(NXCOUR+2)=YPOIN+YV-C*YVECS IF (IDIM.GE.3) XCOOR(NXCOUR+3)=ZPOIN+ZV-C*ZVECS XCOOR(NXCOUR+idimp1)=XCOOR(IREF+idimp1) GOTO 70 C Stockage du point image cree dans XCOOR et dans ICPR C On remplit le maillage IPT2 70 NXCOUR=NXCOUR+idimp1 IPT2.NUM(i,j)=NXCOUR/idimp1 ICPR(IPT1.NUM(i,j))=IPT2.NUM(i,j) ELSE C On recupere le noeud image dans ICPR et on remplit IPT2 IPT2.NUM(i,j)=ICPR(IPT1.NUM(i,j)) ENDIF 11 CONTINUE 10 CONTINUE C Fin de la boucle sur les noeuds de IPT1 C On a cree IPT2 maillage transforme de IPT2 (elements,noeuds) NBPTS=NXCOUR/idimp1 SEGADJ MCOORD RETURN END