1 : C FACTR2 SOURCE CHAT 05/01/12 23:55:58 5004 2 : SUBROUTINE FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4,XDEP2,XARI2 3 : $ ,XN1,XN2,XN3,XINT,ITEST) 4 : *********************************************************************** 5 : *** SP 'FACTR2' : a partir d'une face definie par des pts, permet de 6 : *** dire si un segment (pt depart-pt arrivee) rencontre cette face. 7 : *** 8 : *** APPELES 1 = 'TRPOSE', 'MATMAT', 'SCVECT' (fonction) 9 : *** APPELES 2 = aucun 10 : *** 11 : *** E = 'NDIM' dimension de l'espace 12 : *** 'ITYG' entier caracterisant la geometrie de l'elemnt considere 13 : *** 'PT1', 'PT2', 'PT3', 'PT4' noeuds de la face 14 : *** 'XDEP2' coordonnees reelles de depart de la particule 15 : *** 'XARI2' coordonnees reelles d'arrivee de la particule 16 : *** 'XN1' vecteur unitaire normal à la face 17 : *** 'XN2' vecteur unitaire normal à trajectoire et 'XN1' 18 : *** 'XN3' vecteur unitaire normal à 'XN1' et 'XN2' 19 : *** 'XINT' pt intersection plan associé à face/droite associé à trajectoire 20 : *** 21 : *** S = 'ITEST' vaut 1 si face effectivement traversee, 0 sinon 22 : *** 23 :
Warning: preg_replace(): Compilation failed: unmatched parentheses at offset 6 in C:\xampp\htdocs\cast3m_kit_local_fr_en_2014\sources.php on line 158
24 : *** 25 : *** ORIGINE = CYRIL NOU 26 : *********************************************************************** 27 : IMPLICIT INTEGER(I-N) 28 : IMPLICIT REAL*8 (A-H,O-Z) 29 : -INC CCREEL 30 : DIMENSION XDEP2(3),XARI2(3) 31 : DIMENSION XINT(3),XN1(3),XN2(3),XN3(3) 32 : DIMENSION PT1(3),PT2(3),PT3(3),PT4(3) 33 : DIMENSION X1(3),X2(3),X3(3),X4(4),XPASS(3,3) 34 : DIMENSION PPT1(3),PPT2(3),PPT3(3),PPT4(3),XXINT(3) 35 : DIMENSION PPPT1(3),PPPT2(3),PPPT3(3),PPPT4(3),XXXINT(3) 36 : *** construction des vecteurs directeurs droite,plan, etc... 37 : DO 10 I=1,NDIM 38 : X1(I)=XARI2(I)-XDEP2(I) 39 : X2(I)=XINT(I)-XDEP2(I) 40 : X3(I)=PT2(I)-PT1(I) 41 : X4(I)=XINT(I)-PT1(I) 42 : 10 CONTINUE 43 : ITEST=0 44 : *** 'ITEST1' teste appartenance à trajectoire, 'ITEST2' appartenance à face 45 : ITEST1=0 46 : ITEST2=0 47 : 48 : ************************************************** 49 : *** PT INTERSECTION APPARTIENT À TRAJECTOIRE ? *** 50 : ************************************************** 51 : 52 : *** calcul de la distance entre les pts départ et arrivee 53 : X1NORM=0.D0 54 : DO 20 I=1,NDIM 55 : X1NORM=X1NORM+X1(I)**2 56 : 20 CONTINUE 57 : *** verification si 'XINT' appartient au segment depart-arrivee 58 : SCAL=SCVECT(X2,X1,NDIM) 59 : IF ((SCAL.GT.(-XZPREC*X1NORM)) 60 : $ .AND.(SCAL.LT.(X1NORM*(1+XZPREC)))) ITEST1=1 61 : 62 : ******************************************* 63 : *** PT INTERSECTION APPARTIENT À FACE ? *** 64 : ******************************************* 65 : 66 : ********** CAS 2D ********** 67 : IF (NDIM.EQ.2) THEN 68 : *** calcul de la norme de la face 69 : X3NORM=0.D0 70 : DO 30 I=1,NDIM 71 : X3NORM=X3NORM+X3(I)**2 72 : 30 CONTINUE 73 : *** verification si 'XINT' appartient au segment face 74 : SCAL=SCVECT(X4,X3,NDIM) 75 : IF ((SCAL.GT.(-XZPREC*X3NORM)) 76 : $ .AND.(SCAL.LT.(X3NORM*(1+XZPREC)))) ITEST2=1 77 : 78 : ********** CAS 3D ********** 79 : ELSEIF (NDIM.EQ.3) THEN 80 : *** matrice de passage du repere lié à la face vers repere absolu 81 : DO 40 I=1,NDIM 82 : XPASS(I,1)=XN1(I) 83 : XPASS(I,2)=XN2(I) 84 : XPASS(I,3)=XN3(I) 85 : 40 CONTINUE 86 : *** transpose de 'XPASS' (repere absolu vers repere face) 86 : CALL TRPOSE(XPASS) 87 : *** passage des coordonnees des pts du repere absolu au repere face 87 : CALL MATMAT(XPASS,PT1,NDIM,NDIM,1,PPT1) 87 : CALL MATMAT(XPASS,PT2,NDIM,NDIM,1,PPT2) 87 : CALL MATMAT(XPASS,PT3,NDIM,NDIM,1,PPT3) 87 : CALL MATMAT(XPASS,PT4,NDIM,NDIM,1,PPT4) 87 : CALL MATMAT(XPASS,XINT,NDIM,NDIM,1,XXINT) 88 : *** test par rapport à droite ppt1-ppt2, origine repere face=ppt3 89 : DO 50 I=1,NDIM 90 : PPPT1(I)=PPT1(I)-PPT3(I) 91 : PPPT2(I)=PPT2(I)-PPT3(I) 92 : XXXINT(I)=XXINT(I)-PPT3(I) 93 : 50 CONTINUE 94 : A=(XXXINT(2)-PPPT1(2))*(PPPT2(3)-PPPT1(3)) 95 : B=(XXXINT(3)-PPPT1(3))*(PPPT2(2)-PPPT1(2)) 96 : C=PPPT1(3)*(PPPT2(2)-PPPT1(2)) 97 : D=PPPT1(2)*(PPPT2(3)-PPPT1(3)) 98 : TEST=(A-B)/(C-D) 99 : IF (TEST.GT.(-XZPREC)) THEN 100 : ** test par rapport à droite ppt2-ppt3, origine repere face=ppt1 101 : DO 60 I=1,NDIM 102 : PPPT2(I)=PPT2(I)-PPT1(I) 103 : PPPT3(I)=PPT3(I)-PPT1(I) 104 : XXXINT(I)=XXINT(I)-PPT1(I) 105 : 60 CONTINUE 106 : A=(XXXINT(2)-PPPT2(2))*(PPPT3(3)-PPPT2(3)) 107 : B=(XXXINT(3)-PPPT2(3))*(PPPT3(2)-PPPT2(2)) 108 : C=PPPT2(3)*(PPPT3(2)-PPPT2(2)) 109 : D=PPPT2(2)*(PPPT3(3)-PPPT2(3)) 110 : TEST=(A-B)/(C-D) 111 : IF (TEST.GT.(-XZPREC)) THEN 112 : *** test par rapport à droite ppt3-ppt4, origine repere face=ppt2 113 : IF ((ABS(PPT3(2)-PPT4(2)).GT.(XZPREC*SQRT(X1NORM))) 114 : $ .OR.(ABS(PPT3(3)-PPT4(3)).GT.(XZPREC*SQRT(X1NORM)))) THEN 115 : DO 70 I=1,NDIM 116 : PPPT3(I)=PPT3(I)-PPT2(I) 117 : PPPT4(I)=PPT4(I)-PPT2(I) 118 : XXXINT(I)=XXINT(I)-PPT2(I) 119 : 70 CONTINUE 120 : A=(XXXINT(2)-PPPT3(2))*(PPPT4(3)-PPPT3(3)) 121 : B=(XXXINT(3)-PPPT3(3))*(PPPT4(2)-PPPT3(2)) 122 : C=PPPT3(3)*(PPPT4(2)-PPPT3(2)) 123 : D=PPPT3(2)*(PPPT4(3)-PPPT3(3)) 124 : TEST=(A-B)/(C-D) 125 : IF (TEST.GT.(-XZPREC)) THEN 126 : *** test par rapport à droite ppt4-ppt1, origine repere face=ppt3 127 : DO 80 I=1,NDIM 128 : PPPT4(I)=PPT4(I)-PPT3(I) 129 : PPPT1(I)=PPT1(I)-PPT3(I) 130 : XXXINT(I)=XXINT(I)-PPT3(I) 131 : 80 CONTINUE 132 : A=(XXXINT(2)-PPPT4(2))*(PPPT1(3)-PPPT4(3)) 133 : B=(XXXINT(3)-PPPT4(3))*(PPPT1(2)-PPPT4(2)) 134 : C=PPPT4(3)*(PPPT1(2)-PPPT4(2)) 135 : D=PPPT4(2)*(PPPT1(3)-PPPT4(3)) 136 : TEST=(A-B)/(C-D) 137 : IF (TEST.GT.(-XZPREC)) ITEST2=1 138 : ENDIF 139 : ELSE 140 : *** cas ou ppt3=ppt4 (3 pts seulement definissent la face) 141 : *** test par rapport à droite ppt3-ppt1, origine repere face=ppt2 142 : DO 90 I=1,NDIM 143 : PPPT3(I)=PPT3(I)-PPT2(I) 144 : PPPT1(I)=PPT1(I)-PPT2(I) 145 : XXXINT(I)=XXINT(I)-PPT2(I) 146 : 90 CONTINUE 147 : A=(XXXINT(2)-PPPT3(2))*(PPPT1(3)-PPPT3(3)) 148 : B=(XXXINT(3)-PPPT3(3))*(PPPT1(2)-PPPT3(2)) 149 : C=PPPT3(3)*(PPPT1(2)-PPPT3(2)) 150 : D=PPPT3(2)*(PPPT1(3)-PPPT3(3)) 151 : TEST=(A-B)/(C-D) 152 : IF (TEST.GT.(-XZPREC)) ITEST2=1 153 : ENDIF 154 : ENDIF 155 : ENDIF 156 : ENDIF 157 : *** verification des deux conditions à la fois 158 : IF ((ITEST1.EQ.1).AND.(ITEST2.EQ.1)) ITEST=1 159 : 160 : RETURN 161 : END 162 : 163 : 164 : 165 : 166 : 167 : 168 : 169 :
© Cast3M 2003 - All rights reserved.
Disclaimer