1 : C INTPSEG SOURCE SH236661 13/11/25 21:15:08 7869 2 : C 3 : SUBROUTINE INTPSEG(NRA,NRB,NRC,NRD,X1,Y1,Z1,X2,Y2,Z2, 4 : & IINT,XI1,YI1,ZI1,XI2,YI2,ZI2,XZERO) 5 : C----------------------------------------------------------------------C 6 : C Determination de l'intersection entre un segment et C 7 : C un demi espace definit par un plan et un point dans le demi espace C 8 : C C 9 : C NRA, NRB, NRC = entiers, numeros des noeuds du plan (ABC) C 10 : C NRD = entier, numero du noeud dans le demi espace C 11 : C X1,Y1,Z1 C 12 : C X2,Y2,Z2 = coordonnees des points P1 et P2 du segment [P1,P2] C 13 : C l'ordre P1,P2 est important (voir ci dessous) C 14 : C C 15 : C IINT = indicateur sur l'intersection : C 16 : C - IINT = 0 P1P2 est hors du demi plan C 17 : C --> les coordonnees de I1 et I2 sont nulles C 18 : C - IINT = 1 P1P2 rentre dans le demi plan C 19 : C --> I1 est le point d'intersection C 20 : C --> I2 = P2 C 21 : C - IINT = 2 P1P2 est dans le demi plan C 22 : C --> I1 = P1 C 23 : C --> I2 = P2 C 24 : C - IINT = 3 P1P2 sort du demi plan, cad : C 25 : C --> I1 = P1 C 26 : C --> I2 est le point d'intersection C 27 : C XI1,YI1,ZI1 C 28 : C XI2,YI2,ZI2 = coordonnees des point du segment d'intersection C 29 : C oriente comme le segment [P1,P2] C 30 : C----------------------------------------------------------------------C 31 : IMPLICIT INTEGER(I-N) 32 : IMPLICIT REAL*8 (A-H,O-Z) 33 : -INC CCOPTIO 34 : -INC SMCOORD 35 : c PARAMETER (ZERO=1D-14) 36 : ZERO = XZERO 37 : C 38 : SEGACT,MCOORD 39 : IDIMP1=IDIM+1 40 : C Extraction des coordonnees des points A, B et C 41 : XA=XCOOR((NRA-1)*IDIMP1+1) 42 : YA=XCOOR((NRA-1)*IDIMP1+2) 43 : ZA=XCOOR((NRA-1)*IDIMP1+3) 44 : XB=XCOOR((NRB-1)*IDIMP1+1) 45 : YB=XCOOR((NRB-1)*IDIMP1+2) 46 : ZB=XCOOR((NRB-1)*IDIMP1+3) 47 : XC=XCOOR((NRC-1)*IDIMP1+1) 48 : YC=XCOOR((NRC-1)*IDIMP1+2) 49 : ZC=XCOOR((NRC-1)*IDIMP1+3) 50 : C Vecteur normal au plan ABC : n=AB^AC 51 : XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA)) 52 : YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA)) 53 : ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA)) 54 : C Signe du produit scalaire n.AD 55 : XD=XCOOR((NRD-1)*IDIMP1+1) 56 : YD=XCOOR((NRD-1)*IDIMP1+2) 57 : ZD=XCOOR((NRD-1)*IDIMP1+3) 58 : PSD=(XN*(XD-XA))+(YN*(YD-YA))+(ZN*(ZD-ZA)) 59 : IF (ABS(PSD).LT.ZERO) THEN 60 : C Demi plan mal definit 60 : CALL ERREUR(363) 61 : ENDIF 62 : SIGD=SIGN(1D0,PSD) 63 : C Signe du produit scalaire n.AP1 64 : PS1=(XN*(X1-XA))+(YN*(Y1-YA))+(ZN*(Z1-ZA)) 65 : SIG1=SIGN(1D0,PS1) 66 : IF (ABS(PS1).LT.ZERO) SIG1=SIGD 67 : C Signe du produit scalaire n.AP2 68 : PS2=(XN*(X2-XA))+(YN*(Y2-YA))+(ZN*(Z2-ZA)) 69 : SIG2=SIGN(1D0,PS2) 70 : IF (ABS(PS2).LT.ZERO) SIG2=SIGD 71 : C Il y a 4 cas possibles 72 : IINT=-1 73 : XI1=0D0 74 : YI1=0D0 75 : ZI1=0D0 76 : XI2=0D0 77 : YI2=0D0 78 : ZI2=0D0 79 : C 1) [P1,P2] est hors du demi plan 80 : IF ((SIG1.NE.SIGD).AND.(SIG2.NE.SIGD)) THEN 81 : IINT=0 82 : GOTO 999 83 : ENDIF 84 : C 2) [P1,P2] rentre dans le demi plan 85 : IF ((SIG1.NE.SIGD).AND.(SIG2.EQ.SIGD)) THEN 86 : IINT=1 87 : XI2=X2 88 : YI2=Y2 89 : ZI2=Z2 90 : GOTO 100 91 : ENDIF 92 : C 3) [P1,P2] est entierment dans le demi plan 93 : IF ((SIG1.EQ.SIGD).AND.(SIG2.EQ.SIGD)) THEN 94 : IINT=2 95 : XI1=X1 96 : YI1=Y1 97 : ZI1=Z1 98 : XI2=X2 99 : YI2=Y2 100 : ZI2=Z2 101 : GOTO 999 102 : ENDIF 103 : C 4) [P1,P2] sort du demi plan 104 : IF ((SIG1.EQ.SIGD).AND.(SIG2.NE.SIGD)) THEN 105 : IINT=3 106 : XI1=X1 107 : YI1=Y1 108 : ZI1=Z1 109 : GOTO 100 110 : ENDIF 111 : C 112 : 100 CONTINUE 113 : C Calcul de l'intersection de la droite (P1,P2) avec le plan (ABC) 114 : C Le plan (ABC) est definit par l'equation catesienne 115 : C XN*x + YN*y + ZN*z + W = 0 116 : C calcul de W, sachant que le point A appartient au plan 117 : W=-1D0*((XN*XA)+(YN*YA)+(ZN*ZA)) 118 : C La droite (P1,P2) est definie par le systeme d'equations de 119 : C parametre t suivant 120 : C x = F1*t + F10 121 : C y = F2*t + F20 (F1,F2,F3) est le vecteur directeur de la droite 122 : C z = F3*t + F30 123 : C calcul d'un vecteur directeur de (P1,P2) 124 : F1=X2-X1 125 : F2=Y2-Y1 126 : F3=Z2-Z1 127 : C on peut prendre les coordonnees du point P1 pour les termes 128 : C F10,F20 et F30 puisque ce point est sur la droite 129 : F10=X1 130 : F20=Y1 131 : F30=Z1 132 : C Le point d'intersection de la droite et du plan est le point J 133 : C dont les coordonnees verifient les 4 equations precedentes 134 : C Calcul de la valeur de t (parametre de la droite) pour 135 : C correspondre avec le plan : 136 : C XN*(F1*t+F10) + YN*(F2*t+F20) + ZN*(F3*t+F30) + W = 0 137 : XDENOM=(XN*F1)+(YN*F2)+(ZN*F3) 138 : XNOM=-1D0*((XN*F10)+(YN*F20)+(ZN*F30)+W) 139 : TJ=XNOM/XDENOM 140 : XJ=(F1*TJ)+F10 141 : YJ=(F2*TJ)+F20 142 : ZJ=(F3*TJ)+F30 143 : C Ce point correspond a I1 ou I2 selon la valeur de IINT 144 : IF (IINT.EQ.1) THEN 145 : XI1=XJ 146 : YI1=YJ 147 : ZI1=ZJ 148 : ELSEIF (IINT.EQ.3) THEN 149 : XI2=XJ 150 : YI2=YJ 151 : ZI2=ZJ 152 : ENDIF 153 : C Fin du programme 154 : 999 RETURN 155 : END 156 : 157 : 158 :
© Cast3M 2003 - All rights reserved.
Disclaimer