C DEPLAC SOURCE SP204843 24/03/15 21:15:03 11871 C CE SOUS-PROGRAMME A POUR PROPOS DE DEPLACER L'ENSEMBLE DES POINTS C CONTENUS DANS UN OBJET C 09/2003 : Modifications dans le cas IDIM=1. SUBROUTINE DEPLAC IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCREEL -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMCHPOI CHARACTER*4 MCLE(11) CHARACTER*4 MCL2(3) CHARACTER*4 MCL3(7) CHARACTER*4 MCL4(2) CHARACTER*4 MCL5(2) DIMENSION XU(3),XV(3),XW(3),XP(3) REAL*8 XXX,RAP SEGMENT ICPR(nbpts) DATA NCLE / 11 / DATA MCLE / 'PLUS','MOIN','TOUR','HOMO','AFFI','SYME','PROJ', . 'COOR','MILI','BARS','DEDU' / DATA MCL2 / 'POIN','DROI','PLAN' / DATA MCL3 / 'PLAN','SPHE','CYLI','CONI','TORI','DROI','CERC' / DATA MCL4 / 'DIRE','CONI' / DATA MCL5 / 'CYLI','CART' / C Lecture de l'option de la directive DEPLACER CALL LIRMOT(MCLE,NCLE,ICLE,1) IF (IERR.NE.0) RETURN C Cas IDIM=1, seules les options suivantes sont disponibles : C 'PLUS','MOIN','HOMO','SYME','MILI','DEDU' IF (IDIM.EQ.1) THEN IF ((ICLE.NE.1).AND.(ICLE.NE.2).AND.(ICLE.NE.4).AND. . (ICLE.NE.6).AND.(ICLE.NE.9).AND.(ICLE.NE.11)) THEN MOTERR(1:4)=MCLE(ICLE) INTERR(1)=IDIM CALL ERREUR(971) RETURN ENDIF ENDIF idimp1=IDIM+1 ITYPLU=2 C Option BARS : Elements de BARSOUM C --------------- C Deplacement des noeuds milieu autour de la pointe de fissure IF (ICLE.EQ.10) THEN CALL BARSOU RETURN ENDIF C Option DEDU : C --------------- C Transformation affine definie par au plus (IDIM+1) points IF (ICLE.EQ.11) THEN CALL DEDU(1) RETURN ENDIF CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU) ITYPSA=MELEME IF (IRETOU.NE.1) THEN ITYPLU=1 CALL LIROBJ('POINT',IP,1,IRETOU) IF (IERR.NE.0) RETURN ITYPSA=IP CALL CRELEM(IP) MELEME=IP ENDIF SEGACT MCOORD*MOD C Option MILI : C --------------- IF (ICLE.EQ.9) GOTO 1900 C Remplissage du tableau ICPR des points a deplacer SEGINI ICPR C Est-ce bien utile apres un SEGINI ? c DO i=1,ICPR(/1) c ICPR(i)=0 c ENDDO c petite numerotation locale dans ICPR NICPR=0 SEGACT MELEME IPT1=MELEME DO io=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) THEN IPT1=LISOUS(io) SEGACT IPT1 ENDIF DO i=1,IPT1.NUM(/1) DO j=1,IPT1.NUM(/2) I1=IPT1.NUM(i,j) if(ICPR(I1).eq.0) then NICPR=NICPR+1 ICPR(I1)=NICPR endif ENDDO ENDDO IF (LISOUS(/1).NE.0) SEGDES IPT1 ENDDO SEGDES MELEME C Lecture des donnees supplementaires (selon ICLE) et action C ICLE=1 : PLUS 1 POINT (VECTEUR) C ICLE=2 : MOIN 1 POINT (VECTEUR) C ICLE=3 : TOUR 1 NOMBRE (ANGLE) 1 OU 2 POINTS (AXE 2 OU 3 D) C ICLE=4 : HOMO 1 NOMBRE (RAPPORT) 1 POINT (CENTRE) C ICLE=5 : AFFI 1 NOMBRE (RAPPORT) 1 POINT (INVARIANT) C 1 POINT(VECTEUR) C ICLE=6 : SYME 1 MOT (POIN DROI OU PLAN) 1,2 OU 3 POINTS C ICLE=7 : PROJ 1 MOT (PLAN SPHE CYLI CONI TORI DROI CERC) C DES DONNEES EN RAPPORT C ICLE=8 : COOR 1 MOT (CART CYLI) GOTO (100,200,300,400,500,600,700,1800),ICLE C Option PLUS : C --------------- 100 XSENS=1.D0 GOTO 201 C Option MOIN : C --------------- 200 XSENS=-1.D0 201 CALL LIROBJ('POINT',IP,0,IRETOU) IF (IERR.NE.0) GOTO 10000 IF (IRETOU.EQ.0) THEN CALL LIROBJ('CHPOINT ',IPCH,1,IRETOU) IF (IERR.NE.0) GOTO 10000 CALL DEPCHP(ICPR,IPCH,XSENS) GOTO 1000 ENDIF IREF=idimp1*(IP-1) DO j=1,IDIM XV(j)=XCOOR(IREF+j)*XSENS ENDDO DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) DO j=1,IDIM XCOOR(IREF+j)=XCOOR(IREF+j)+XV(j) ENDDO ENDIF ENDDO GOTO 1000 C Option TOUR : C --------------- c lecture de l'angle 300 CALL LIRREE(ANG,0,IRETOU) IF (IERR.NE.0) GOTO 10000 IPCH1=0 C -- bp (nouvelle option 06/2013) : TOUR CHPOint -- IF (IRETOU.EQ.0) THEN CALL LIROBJ('CHPOINT ',IPCH1,1,IRET1) IF (IERR.NE.0) GOTO 10000 c creation d'un MPOVAL coherent avec ICPR c avec comme composantes COS et SIN N=NICPR NC=2 SEGINI,MPOVAL DO I=1,N VPOCHA(I,1)=1.D0 c VPOCHA(I,2)=0.D0 =sin(0) ENDDO c remplissage du MPOVAL depuis le MCHPO1 MCHPO1=IPCH1 SEGACT,MCHPO1 DO IP=1,MCHPO1.IPCHP(/1) MSOUP1=MCHPO1.IPCHP(IP) SEGACT,MSOUP1 IPT1=MSOUP1.IGEOC SEGACT,IPT1 MPOVA1=MSOUP1.IPOVAL SEGACT,MPOVA1 N1 = MPOVA1.VPOCHA(/1) NC1 = MPOVA1.VPOCHA(/2) IF (NC1.ne.1) THEN call ERREUR(180) GOTO 10000 ENDIF DO 303 J1=1,N1 I = ICPR(IPT1.NUM(1,J1)) IF(I.eq.0) GOTO 303 ANG=MPOVA1.VPOCHA(J1,1)/180.D0*XPI VPOCHA(I,1)=COS(ANG) VPOCHA(I,2)=SIN(ANG) 303 CONTINUE SEGDES,MPOVA1,IPT1 SEGDES,MSOUP1 ENDDO SEGDES,MCHPO1 ELSE ANG=ANG/180.D0*XPI CO=COS(ANG) SI=SIN(ANG) ENDIF c lecture des points centre et axe + fabrication repere de rotation CALL LIROBJ('POINT',IP1,1,IRETOU) IF (IDIM.EQ.3) CALL LIROBJ('POINT',IP2,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP1-1) XP(1)=XCOOR(IREF+1) XP(2)=XCOOR(IREF+2) IF (IDIM.LT.3) THEN XP(3)=0.D0 XW(1)=0.D0 XW(2)=0.D0 XW(3)=1.D0 ELSE XP(3)=XCOOR(IREF+3) IREF=idimp1*(IP2-1) XW(1)=XCOOR(IREF+1)-XP(1) XW(2)=XCOOR(IREF+2)-XP(2) XW(3)=XCOOR(IREF+3)-XP(3) XN=SQRT(XW(1)**2+XW(2)**2+XW(3)**2) IF (XN.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XW(1)=XW(1)/XN XW(2)=XW(2)/XN XW(3)=XW(3)/XN ENDIF XU(1)=-XW(2) XU(2)=XW(1) XN=SQRT(XU(1)**2+XU(2)**2) IF (XN.GE.0.1D0) THEN XU(3)=0.D0 XU(1)=XU(1)/XN XU(2)=XU(2)/XN ELSE XU(1)=0.D0 XU(2)=-XW(3) XU(3)=XW(2) XN=SQRT(XU(2)**2+XU(3)**2) XU(2)=XU(2)/XN XU(3)=XU(3)/XN ENDIF XV(1)=XW(2)*XU(3)-XW(3)*XU(2) XV(2)=XW(3)*XU(1)-XW(1)*XU(3) XV(3)=XW(1)*XU(2)-XW(2)*XU(1) c boucle sur les noeuds + rotation DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) XD=XCOOR(IREF+1)-XP(1) YD=XCOOR(IREF+2)-XP(2) ZD=0.D0 IF (IDIM.EQ.3) ZD=XCOOR(IREF+3)-XP(3) XE=XD*XU(1)+YD*XU(2)+ZD*XU(3) YE=XD*XV(1)+YD*XV(2)+ZD*XV(3) ZE=XD*XW(1)+YD*XW(2)+ZD*XW(3) IF(IPCH1.ne.0) THEN CO = VPOCHA(ICPR(i),1) SI = VPOCHA(ICPR(i),2) ENDIF XD=XE*CO-YE*SI YD=XE*SI+YE*CO ZD=ZE XCOOR(IREF+1)=XD*XU(1)+YD*XV(1)+ZD*XW(1)+XP(1) XCOOR(IREF+2)=XD*XU(2)+YD*XV(2)+ZD*XW(2)+XP(2) IF (IDIM.EQ.3) XCOOR(IREF+3)=XD*XU(3)+YD*XV(3)+ZD*XW(3)+XP(3) ENDIF ENDDO IF(IPCH1.ne.0) SEGSUP,MPOVAL GOTO 1000 C Option HOMO : C --------------- 400 CALL LIRREE(RAP,1,IRETOU) CALL LIROBJ('POINT',IP,1,IRETOU) IF (RAP.EQ.0.D0) CALL ERREUR(36) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP-1) DO j=1,IDIM XP(j)=XCOOR(IREF+j) ENDDO DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) DO j=1,IDIM XCOOR(IREF+j)=XP(j)+RAP*(XCOOR(IREF+j)-XP(j)) ENDDO ENDIF ENDDO GOTO 1000 C Option AFFI : C --------------- 500 CALL LIRREE(RAP,1,IRETOU) CALL LIROBJ('POINT',IPC,1,IRETOU) CALL LIROBJ('POINT',IPV,1,IRETOU) IF (RAP.EQ.0.D0) CALL ERREUR(36) IF (IERR.NE.0) GOTO 10000 RAP=RAP-1. IREF=idimp1*(IPC-1) XP(1)=XCOOR(IREF+1) XP(2)=XCOOR(IREF+2) XP(3)=0.D0 IF (IDIM.EQ.3) XP(3)=XCOOR(IREF+3) IREF=idimp1*(IPV-1) XV(1)=XCOOR(IREF+1)-XP(1) XV(2)=XCOOR(IREF+2)-XP(2) XV(3)=0.D0 IF (IDIM.EQ.3) XV(3)=XCOOR(IREF+3)-XP(3) XN=XV(1)**2+XV(2)**2+XV(3)**2 IF (XN.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XN=SQRT(XN) XV(1)=XV(1)/XN XV(2)=XV(2)/XN XV(3)=XV(3)/XN DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) XU(1)=XCOOR(IREF+1)-XP(1) XU(2)=XCOOR(IREF+2)-XP(2) XU(3)=0.D0 IF (IDIM.EQ.3) XU(3)=XCOOR(IREF+3)-XP(3) SCA=(XU(1)*XV(1)+XU(2)*XV(2)+XU(3)*XV(3))*RAP DO j=1,IDIM XCOOR(IREF+j)=XP(j)+XU(j)+SCA*XV(j) ENDDO ENDIF ENDDO GOTO 1000 C Option SYME : C --------------- 600 CALL LIRMOT(MCL2,3,JCLE,1) IF (IERR.NE.0) GOTO 10000 IF ((IDIM.EQ.1).AND.(JCLE.NE.1)) THEN MOTERR(1:4)=MCL2(JCLE) INTERR(1)=IDIM CALL ERREUR(971) GOTO 10000 ENDIF GOTO (610,620,630),JCLE C Option SYME POIN : C -------------------- 610 CALL LIROBJ('POINT',IP,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP-1) DO j=1,IDIM XP(j)=2*XCOOR(IREF+j) ENDDO DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) DO j=1,IDIM XCOOR(IREF+j)=XP(j)-XCOOR(IREF+j) ENDDO ENDIF ENDDO GOTO 1000 C Option SYME DROI : C -------------------- 620 CALL LIROBJ('POINT',IP1,1,IRETOU) CALL LIROBJ('POINT',IP2,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP1-1) XP(1)=XCOOR(IREF+1) XP(2)=XCOOR(IREF+2) XP(3)=0.D0 IF (IDIM.EQ.3) XP(3)=XCOOR(IREF+3) IREF=idimp1*(IP2-1) XU(1)=XCOOR(IREF+1)-XP(1) XU(2)=XCOOR(IREF+2)-XP(2) XU(3)=0.D0 IF (IDIM.EQ.3) XU(3)=XCOOR(IREF+3)-XP(3) XN=XU(1)**2+XU(2)**2+XU(3)**2 IF (XN.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XN=SQRT(XN) XU(1)=XU(1)/XN XU(2)=XU(2)/XN XU(3)=XU(3)/XN DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) XV(1)=XCOOR(IREF+1)-XP(1) XV(2)=XCOOR(IREF+2)-XP(2) XV(3)=0.D0 IF (IDIM.GE.3) XV(3)=XCOOR(IREF+3)-XP(3) SCA=2*(XU(1)*XV(1)+XU(2)*XV(2)+XU(3)*XV(3)) DO j=1,IDIM XCOOR(IREF+j)=XP(j)+SCA*XU(j)-XV(j) ENDDO ENDIF ENDDO GOTO 1000 C Option SYME PLAN : C -------------------- 630 IF (IDIM.EQ.2) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 CALL LIROBJ('POINT',IP1,1,IRETOU) CALL LIROBJ('POINT',IP2,1,IRETOU) CALL LIROBJ('POINT',IP3,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP1-1) XP(1)=XCOOR(IREF+1) XP(2)=XCOOR(IREF+2) XP(3)=XCOOR(IREF+3) IREF=idimp1*(IP2-1) XU(1)=XCOOR(IREF+1)-XP(1) XU(2)=XCOOR(IREF+2)-XP(2) XU(3)=XCOOR(IREF+3)-XP(3) IREF=idimp1*(IP3-1) XV(1)=XCOOR(IREF+1)-XP(1) XV(2)=XCOOR(IREF+2)-XP(2) XV(3)=XCOOR(IREF+3)-XP(3) XW(1)=XU(2)*XV(3)-XU(3)*XV(2) XW(2)=XU(3)*XV(1)-XU(1)*XV(3) XW(3)=XU(1)*XV(2)-XU(2)*XV(1) XN=XW(1)**2+XW(2)**2+XW(3)**2 IF (XN.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XN=SQRT(XN) XW(1)=XW(1)/XN XW(2)=XW(2)/XN XW(3)=XW(3)/XN DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=idimp1*(i-1) XU(1)=XCOOR(IREF+1)-XP(1) XU(2)=XCOOR(IREF+2)-XP(2) XU(3)=XCOOR(IREF+3)-XP(3) SCA=2*(XU(1)*XW(1)+XU(2)*XW(2)+XU(3)*XW(3)) XCOOR(IREF+1)=XP(1)+XU(1)-SCA*XW(1) XCOOR(IREF+2)=XP(2)+XU(2)-SCA*XW(2) XCOOR(IREF+3)=XP(3)+XU(3)-SCA*XW(3) ENDIF ENDDO GOTO 1000 C Option SYME PROJ : C -------------------- 700 CALL LIRMOT(MCL4,2,ITYYP,0) C Lecture du vecteur de projection ou du centre (projection conique) CALL LIROBJ('POINT',IP1,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP1-1) XVECS=XCOOR(IREF+1) YVECS=XCOOR(IREF+2) ZVECS=0.D0 IF (IDIM.GE.3) ZVECS=XCOOR(IREF+3) ICARR=1 IF (ITYYP.NE.0) ICARR=0 CALL LIRMOT(MCL3,7,JCLE,ICARR) IF (IERR.NE.0) GOTO 10000 IF (JCLE.EQ.0) THEN JCLE=ITYYP+2 ITYYP=0 ENDIF IF ((IDIM.EQ.2).AND.(JCLE.LE.5)) CALL ERREUR(34) IF ((IDIM.EQ.3).AND.(JCLE.GE.6)) CALL ERREUR(34) IF (IERR.NE.0) GOTO 10000 C Projection cylindrique normalisation du vecteur IF (ITYYP.NE.2) THEN SVECS=SQRT(XVECS*XVECS+YVECS*YVECS+ZVECS*ZVECS) IF (SVECS.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XVEC=XVECS/SVECS YVEC=YVECS/SVECS ZVEC=ZVECS/SVECS ENDIF GOTO (710,720,730,740,750,760,770),JCLE C Option PROJ PLAN : C -------------------- 710 CALL LIROBJ('POINT',IP1,1,IRETOU) CALL LIROBJ('POINT',IP2,1,IRETOU) CALL LIROBJ('POINT',IP3,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=idimp1*(IP1-1) XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=idimp1*(IP2-1) XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 IREF=idimp1*(IP3-1) XV3=XCOOR(IREF+1)-XPT1 YV3=XCOOR(IREF+2)-YPT1 ZV3=XCOOR(IREF+3)-ZPT1 XV1=YV2*ZV3-ZV2*YV3 YV1=ZV2*XV3-XV2*ZV3 ZV1=XV2*YV3-YV2*XV3 SV1=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (SV1.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 GOTO 780 C Option PROJ SPHE : C -------------------- 720 CALL LIROBJ('POINT',IPCEN,1,IRETOU) CALL LIROBJ('POINT',IP1,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=(IPCEN-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IP1-1)*idimp1 XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 ANGLE=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (ANGLE.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 GOTO 780 C Option PROJ CYLI : C -------------------- 730 CALL LIROBJ('POINT',IPOIN1,1,IRETOU) CALL LIROBJ('POINT',IPOIN2,1,IRETOU) CALL LIROBJ('POINT',IP1,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=(IPOIN1-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IPOIN2-1)*idimp1 XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 S=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (S.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV1=XV1/S YV1=YV1/S ZV1=ZV1/S IREF=(IP1-1)*idimp1 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*XV3+YV3*YV3+ZV3*ZV3) IF (ANGLE.EQ.0) CALL ERREUR(21) IF (IERR.NE.0) RETURN GOTO 780 C Option PROJ CONI : C -------------------- 740 CALL LIROBJ('POINT',IPT1,1,IRETOU) CALL LIROBJ('POINT',IP1,1,IRETOU) CALL LIROBJ('POINT',IP2,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=(IPT1-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IP1-1)*idimp1 XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 SV1=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (SV1.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 IREF=(IP2-1)*idimp1 XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=XCOOR(IREF+3)-ZPT1 SV2=SQRT(XV2*XV2+YV2*YV2+ZV2*ZV2) IF (SV2.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV2=XV2/SV2 YV2=YV2/SV2 ZV2=ZV2/SV2 ANGLE=(XV1*XV2+YV1*YV2+ZV1*ZV2)**2 GOTO 780 C Option PROJ TORI C ------------------ 750 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) GOTO 10000 IREF=(IPT1-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=XCOOR(IREF+3) IREF=(IPAX-1)*idimp1 XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=XCOOR(IREF+3)-ZPT1 SV1=XV1*XV1+YV1*YV1+ZV1*ZV1 IF (SV1.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 SV1=SQRT(SV1) XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 IREF=(IPCS-1)*idimp1 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*XV2+YV2*YV2+ZV2*ZV2 IF (GR2.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 IREF=(IP1-1)*idimp1 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*XC+YC*YC+ZC*ZC IF (SC.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 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.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV2=PR2 YV2=GR2 GOTO 780 C Option PROJ DROI : C -------------------- 760 CALL LIROBJ('POINT',IP1,1,IRETOU) CALL LIROBJ('POINT',IP2,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=(IP1-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=0.D0 IREF=(IP2-1)*idimp1 XV2=XCOOR(IREF+1)-XPT1 YV2=XCOOR(IREF+2)-YPT1 ZV2=0.D0 XV3=0.D0 YV3=0.D0 ZV3=1.D0 XV1=YV2*ZV3-ZV2*YV3 YV1=ZV2*XV3-XV2*ZV3 ZV1=XV2*YV3-YV2*XV3 SV1=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (SV1.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 XV1=XV1/SV1 YV1=YV1/SV1 ZV1=ZV1/SV1 GOTO 780 C Option PROJ CERC : C -------------------- 770 CALL LIROBJ('POINT',IPCEN,1,IRETOU) CALL LIROBJ('POINT',IP1,1,IRETOU) IF (IERR.NE.0) GOTO 10000 IREF=(IPCEN-1)*idimp1 XPT1=XCOOR(IREF+1) YPT1=XCOOR(IREF+2) ZPT1=0.D0 IREF=(IP1-1)*idimp1 XV1=XCOOR(IREF+1)-XPT1 YV1=XCOOR(IREF+2)-YPT1 ZV1=0.D0 ANGLE=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) IF (ANGLE.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 GOTO 780 C Traitement commun - Options PROJ 780 DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=(i-1)*idimp1 XPOIN=XCOOR(IREF+1) YPOIN=XCOOR(IREF+2) ZPOIN=0.D0 IF (IDIM.GE.3) ZPOIN=XCOOR(IREF+3) IF (ITYYP.EQ.2) THEN XVEC=XVECS-XPOIN YVEC=YVECS-YPOIN ZVEC=ZVECS-ZPOIN SVEC=SQRT(XVEC*XVEC+YVEC*YVEC+ZVEC*ZVEC) IF (SVEC.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) RETURN XVEC=XVEC/SVEC YVEC=YVEC/SVEC ZVEC=ZVEC/SVEC ENDIF GOTO (810,820,830,840,850,860,870),JCLE C Option PROJ PLAN 810 XV2=XPOIN-XPT1 YV2=YPOIN-YPT1 ZV2=ZPOIN-ZPT1 DENUM=XV2*XV1+YV2*YV1+ZV2*ZV1 DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 IF (DENOM.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 RAP=-DENUM/DENOM XPOIN=XPOIN+RAP*XVEC YPOIN=YPOIN+RAP*YVEC ZPOIN=ZPOIN+RAP*ZVEC GOTO 880 C Option PROJ SPHE 820 XV1=XPT1-XPOIN YV1=YPT1-YPOIN ZV1=ZPT1-ZPOIN SCA=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 XV1=XVEC*SCA YV1=YVEC*SCA ZV1=ZVEC*SCA S2=(XPOIN+XV1-XPT1)**2+(YPOIN+YV1-YPT1)**2+(ZPOIN+ZV1-ZPT1)**2 IF (S2.GT.(ANGLE*ANGLE)) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 C=SQRT(ANGLE*ANGLE-S2) IF (SCA.LT.0.D0) C=-C XPOIN=XPOIN+XV1-C*XVEC YPOIN=YPOIN+YV1-C*YVEC ZPOIN=ZPOIN+ZV1-C*ZVEC GOTO 880 C Option PROJ CYLI 830 XV2=YVEC*ZV1-ZVEC*YV1 YV2=ZVEC*XV1-XVEC*ZV1 ZV2=XVEC*YV1-YVEC*XV1 S2V2=XV2*XV2+YV2*YV2+ZV2*ZV2 IF (S2V2.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 C2=(XVEC*XV1+YVEC*YV1+ZVEC*ZV1)**2 IF (C2.EQ.1.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 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 (DIS2.GT.ANGLE**2) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 DMU=SQRT((ANGLE**2-DIS2)/(1.D0-C2)) DNUM=XV2*XV4+YV2*YV4+ZV2*ZV4 DLA=DNUM/S2V2 DMU=SIGN(DMU,DLA) XPOIN=XPOIN+XVEC*(DLA-DMU) YPOIN=YPOIN+YVEC*(DLA-DMU) ZPOIN=ZPOIN+ZVEC*(DLA-DMU) GOTO 880 C Option PROJ CONI 840 XV2=XPOIN-XPT1 YV2=YPOIN-YPT1 ZV2=ZPOIN-ZPT1 VECV1=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 VEC2=XVEC*XVEC+YVEC*YVEC+ZVEC*ZVEC V2V1=XV2*XV1+YV2*YV1+ZV2*ZV1 VECV2=XVEC*XV2+YVEC*YV2+ZVEC*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.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 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 880 C Option PROJ TORI 850 PR2=XV2 GR2=YV2 XOP=XPOIN-XPT1 YOP=YPOIN-YPT1 ZOP=ZPOIN-ZPT1 OPV=XOP*XVEC+YOP*YVEC+ZOP*ZVEC R2=XOP*XOP+YOP*YOP+ZOP*ZOP-GR2-PR2 VA=XVEC*XV1+YVEC*YV1+ZVEC*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.D0 C Resolution de l'equation du 4eme degre par iteration 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.D0) CALL ERREUR(40) IF (IERR.NE.0) GOTO 10000 CORR=FDLM/FPDLM RLMD=RLMD-CORR IF ((RLMD.NE.0.D0).AND.(ABS(CORR/RLMD).LT.1E-5)) GOTO 851 ENDDO CALL ERREUR(40) GOTO 10000 851 XPOIN=XPOIN+RLMD*XVEC YPOIN=YPOIN+RLMD*YVEC ZPOIN=ZPOIN+RLMD*ZVEC GOTO 880 C Option PROJ DROI 860 XV2=XPOIN-XPT1 YV2=YPOIN-YPT1 ZV2=0.D0 DENUM=XV2*XV1+YV2*YV1+ZV2*ZV1 DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 IF (DENOM.EQ.0.D0) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 RAP=-DENUM/DENOM XPOIN=XPOIN+RAP*XVEC YPOIN=YPOIN+RAP*YVEC ZPOIN=ZPOIN+RAP*ZVEC GOTO 880 C Option PROJ CERC 870 XV1=XPT1-XPOIN YV1=YPT1-YPOIN ZV1=0.D0 SCA=XVEC*XV1+YVEC*YV1+ZVEC*ZV1 XV1=XVEC*SCA YV1=YVEC*SCA ZV1=ZVEC*SCA S2=(XPOIN+XV1-XPT1)**2+(YPOIN+YV1-YPT1)**2+(ZPOIN+ZV1-ZPT1)**2 IF (S2.GT.ANGLE**2) CALL ERREUR(21) IF (IERR.NE.0) GOTO 10000 C=SQRT(ANGLE*ANGLE-S2) IF (SCA.LT.0.D0) C=-C XPOIN=XPOIN+XV1-C*XVEC YPOIN=YPOIN+YV1-C*YVEC ZPOIN=ZPOIN+ZV1-C*ZVEC GOTO 880 C Traitement commun des options 880 XCOOR(IREF+1)=XPOIN XCOOR(IREF+2)=YPOIN IF (IDIM.EQ.3) XCOOR(IREF+3)=ZPOIN ENDIF ENDDO GOTO 1000 C Option COOR : Changement de systeme de coordonnees C --------------- 1800 CALL LIRMOT(MCL5,2,ICYLI,1) IF (IERR.NE.0) GOTO 10000 IF (ICYLI.EQ.1) THEN CALL LIROBJ('POINT ',IP1,1,IRETOU) IF (IERR.NE.0) GO TO 10000 CALL LIROBJ('POINT ',IP2,1,IRETOU) IF (IERR.NE.0) GO TO 10000 IF (IDIM.EQ.3) CALL LIROBJ('POINT ',IP3,1,IRETOU) IF (IERR.NE.0) GO TO 10000 ENDIF C Transformation des coordonnees cartesiennes en cylindriques IF (ICYLI.EQ.1) THEN IREF=(IP1-1)*idimp1 XP11=XCOOR(IREF+1) XP12=XCOOR(IREF+2) XP13=0.D0 IF (IDIM.EQ.3) XP13=XCOOR(IREF+3) IREF=(IP2-1)*idimp1 XP21=XCOOR(IREF+1)-XP11 XP22=XCOOR(IREF+2)-XP12 XP23=0.D0 IF(IDIM.EQ.3) XP23 = XCOOR(IREF+3)-XP13 XL=SQRT(XP21*XP21+XP22*XP22+XP23*XP23) IF (XL.EQ.0.D0) THEN CALL ERREUR(21) GOTO 10000 ENDIF XP21=XP21/XL XP22=XP22/XL YP21=-XP22 YP22=XP21 YP23=0.D0 IF (IDIM.EQ.3) THEN XP23=XP23/XL XP31=XP21 XP32=XP22 XP33=XP23 IREF=(IP3-1)*idimp1 XP21=XCOOR(IREF+1)-XP11 XP22=XCOOR(IREF+2)-XP12 XP23=XCOOR(IREF+3)-XP13 XXX=XP21*XP31+XP22*XP32+XP23*XP33 XP21=XP21-XXX*XP31 XP22=XP22-XXX*XP32 XP23=XP23-XXX*XP33 XL=SQRT(XP21*XP21+XP22*XP22+XP23*XP23) IF (XL.EQ.0.D0) THEN CALL ERREUR(21) GOTO 10000 ENDIF XP21=XP21/XL XP22=XP22/XL XP23=XP23/XL YP21=XP32*XP23-XP33*XP22 YP22=XP33*XP21-XP31*XP23 YP23=XP31*XP22-XP32*XP21 ENDIF SCA=180.D0/XPI DO i=1,ICPR(/1) IF (ICPR(i).NE.0) THEN IREF=(i-1)*idimp1 XD=XCOOR(IREF+1)-XP11 YD=XCOOR(IREF+2)-XP12 ZD=0.D0 IF (IDIM.EQ.3) THEN ZD=XCOOR(IREF+3)-XP13 ZF=XD*XP31+YD*XP32+ZD*XP33 XD=XD-ZF*XP31 YD=YD-ZF*XP32 ZD=ZD-ZF*XP33 ENDIF XA=XD*XP21+YD*XP22+ZD*XP23 YA=XD*YP21+YD*YP22+ZD*YP23 RHO=SQRT(XA*XA+YA*YA) TETA=ATAN2(YA,XA) XCOOR(IREF+1)=RHO XCOOR(IREF+2)=TETA*SCA IF (IDIM.EQ.3) XCOOR(IREF+3)=ZF ENDIF ENDDO C Transformation des coordonnees cylindriques en cartesiennes C L'axe X est conserve, ainsi que l'axe Z. ELSE IF (ICYLI.EQ.2) THEN SCA=XPI/180.D0 DO i=1,ICPR(/1) IF (ICPR(I).NE.0) THEN IREF=(i-1)*idimp1 XD=XCOOR(IREF+1) YD=XCOOR(IREF+2)*SCA XCOOR(IREF+1)=XD*COS(YD) XCOOR(IREF+2)=XD*SIN(YD) ENDIF ENDDO ENDIF GOTO 1000 C Fin normale de l'operateur - Ecriture des objets (ICLE = 1 a 8) C ----------------------------------------------------------------- 1000 SEGSUP ICPR IF (ITYPLU.EQ.1) CALL ECROBJ('POINT',ITYPSA) IF (ITYPLU.EQ.2) CALL ECROBJ('MAILLAGE',ITYPSA) RETURN C Option MILI : C --------------- 1900 IF (IDIM.GE.2) THEN CALL ECROBJ('MAILLAGE',MELEME) CALL CHANLG CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU) IF (IERR.NE.0) RETURN ENDIF SEGACT MELEME IPT1=MELEME INSO=LISOUS(/1) DO j=1,MAX(1,INSO) IF (INSO.NE.0) THEN IPT1=LISOUS(j) SEGACT IPT1 ENDIF IF (IPT1.ITYPEL.EQ.3) THEN DO i=1,IPT1.NUM(/2) IRef1=idimp1*(IPT1.NUM(1,i)-1) IRef2=idimp1*(IPT1.NUM(3,i)-1) IRef3=idimp1*(IPT1.NUM(2,i)-1) X1=XCOOR(IRef1+1) X2=XCOOR(IRef2+1) X3=XCOOR(IRef3+1) X2X1=X2-X1 X3M=(X2+X1)*0.5-X3 XNUM=X2X1*X3M XDEN=X2X1*X2X1 IF (IDIM.GE.2) THEN Y1=XCOOR(IRef1+2) Y2=XCOOR(IRef2+2) Y3=XCOOR(IRef3+2) Y2Y1=Y2-Y1 Y3M=(Y2+Y1)*0.5-Y3 XNUM=XNUM+Y2Y1*Y3M XDEN=XDEN+Y2Y1*Y2Y1 IF (IDIM.GE.3) THEN Z1=XCOOR(IRef1+3) Z2=XCOOR(IRef2+3) Z3=XCOOR(IRef3+3) Z2Z1=Z2-Z1 Z3M=(Z2+Z1)*0.5-Z3 XNUM=XNUM+Z2Z1*Z3M XDEN=XDEN+Z2Z1*Z2Z1 ENDIF ENDIF IF (XDEN.LT.1.D-20) THEN wrITE (6,*) 'noeud double ' , IPT1.NUM(1,i) xnum = 0.D0 xden = 1.D0 c RETURN endif XLAMBD=XNUM/XDEN XCOOR(IRef3+1)=X3+XLAMBD*X2X1 IF (IDIM.GE.2) THEN XCOOR(IRef3+2)=Y3+XLAMBD*Y2Y1 IF (IDIM.GE.3) THEN XCOOR(IRef3+3)=Z3+XLAMBD*Z2Z1 ENDIF ENDIF ENDDO ENDIF SEGDES IPT1 ENDDO SEGDES MELEME RETURN C Traitement des erreurs (options ICLE = 1 a 8) C ----------------------------------------------- 10000 SEGSUP ICPR RETURN END