Télécharger la plaquette Cast3M


fr en

Télécharger factr2.eso

Retour à la liste

Numérotation des lignes : 

factr2

  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. *** Rq : 'XZPREC' (-INC CCREEL) erreur precision calcul machine
  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)
  87.          CALL TRPOSE(XPASS)
  88. ***      passage des coordonnees des pts du repere absolu au repere face
  89.          CALL MATMAT(XPASS,PT1,NDIM,NDIM,1,PPT1)
  90.          CALL MATMAT(XPASS,PT2,NDIM,NDIM,1,PPT2)
  91.          CALL MATMAT(XPASS,PT3,NDIM,NDIM,1,PPT3)
  92.          CALL MATMAT(XPASS,PT4,NDIM,NDIM,1,PPT4)
  93.          CALL MATMAT(XPASS,XINT,NDIM,NDIM,1,XXINT)
  94. ***      test par rapport à droite ppt1-ppt2, origine repere face=ppt3
  95.          DO 50 I=1,NDIM
  96.             PPPT1(I)=PPT1(I)-PPT3(I)
  97.             PPPT2(I)=PPT2(I)-PPT3(I)
  98.             XXXINT(I)=XXINT(I)-PPT3(I)
  99.  50      CONTINUE
  100.          A=(XXXINT(2)-PPPT1(2))*(PPPT2(3)-PPPT1(3))
  101.          B=(XXXINT(3)-PPPT1(3))*(PPPT2(2)-PPPT1(2))
  102.          C=PPPT1(3)*(PPPT2(2)-PPPT1(2))
  103.          D=PPPT1(2)*(PPPT2(3)-PPPT1(3))
  104.          TEST=(A-B)/(C-D)
  105.          IF (TEST.GT.(-XZPREC)) THEN
  106. **          test par rapport à droite ppt2-ppt3, origine repere face=ppt1
  107.             DO 60 I=1,NDIM
  108.                PPPT2(I)=PPT2(I)-PPT1(I)
  109.                PPPT3(I)=PPT3(I)-PPT1(I)
  110.                XXXINT(I)=XXINT(I)-PPT1(I)
  111.  60         CONTINUE
  112.             A=(XXXINT(2)-PPPT2(2))*(PPPT3(3)-PPPT2(3))
  113.             B=(XXXINT(3)-PPPT2(3))*(PPPT3(2)-PPPT2(2))
  114.             C=PPPT2(3)*(PPPT3(2)-PPPT2(2))
  115.             D=PPPT2(2)*(PPPT3(3)-PPPT2(3))
  116.             TEST=(A-B)/(C-D)
  117.             IF (TEST.GT.(-XZPREC)) THEN
  118. ***            test par rapport à droite ppt3-ppt4, origine repere face=ppt2
  119.                IF ((ABS(PPT3(2)-PPT4(2)).GT.(XZPREC*SQRT(X1NORM)))
  120.      $         .OR.(ABS(PPT3(3)-PPT4(3)).GT.(XZPREC*SQRT(X1NORM)))) THEN
  121.                   DO 70 I=1,NDIM
  122.                      PPPT3(I)=PPT3(I)-PPT2(I)
  123.                      PPPT4(I)=PPT4(I)-PPT2(I)
  124.                      XXXINT(I)=XXINT(I)-PPT2(I)
  125.  70               CONTINUE
  126.                   A=(XXXINT(2)-PPPT3(2))*(PPPT4(3)-PPPT3(3))
  127.                   B=(XXXINT(3)-PPPT3(3))*(PPPT4(2)-PPPT3(2))
  128.                   C=PPPT3(3)*(PPPT4(2)-PPPT3(2))
  129.                   D=PPPT3(2)*(PPPT4(3)-PPPT3(3))
  130.                   TEST=(A-B)/(C-D)
  131.                   IF (TEST.GT.(-XZPREC)) THEN
  132. ***                  test par rapport à droite ppt4-ppt1, origine repere face=ppt3
  133.                      DO 80 I=1,NDIM
  134.                         PPPT4(I)=PPT4(I)-PPT3(I)
  135.                         PPPT1(I)=PPT1(I)-PPT3(I)
  136.                         XXXINT(I)=XXINT(I)-PPT3(I)
  137.  80                  CONTINUE
  138.                      A=(XXXINT(2)-PPPT4(2))*(PPPT1(3)-PPPT4(3))
  139.                      B=(XXXINT(3)-PPPT4(3))*(PPPT1(2)-PPPT4(2))
  140.                      C=PPPT4(3)*(PPPT1(2)-PPPT4(2))
  141.                      D=PPPT4(2)*(PPPT1(3)-PPPT4(3))
  142.                      TEST=(A-B)/(C-D)
  143.                      IF (TEST.GT.(-XZPREC)) ITEST2=1
  144.                   ENDIF
  145.                ELSE
  146. ***               cas ou ppt3=ppt4 (3 pts seulement definissent la face)
  147. ***               test par rapport à droite ppt3-ppt1, origine repere face=ppt2
  148.                   DO 90 I=1,NDIM
  149.                      PPPT3(I)=PPT3(I)-PPT2(I)
  150.                      PPPT1(I)=PPT1(I)-PPT2(I)
  151.                      XXXINT(I)=XXINT(I)-PPT2(I)
  152.  90               CONTINUE
  153.                   A=(XXXINT(2)-PPPT3(2))*(PPPT1(3)-PPPT3(3))
  154.                   B=(XXXINT(3)-PPPT3(3))*(PPPT1(2)-PPPT3(2))
  155.                   C=PPPT3(3)*(PPPT1(2)-PPPT3(2))
  156.                   D=PPPT3(2)*(PPPT1(3)-PPPT3(3))
  157.                   TEST=(A-B)/(C-D)
  158.                   IF (TEST.GT.(-XZPREC)) ITEST2=1
  159.                ENDIF
  160.             ENDIF
  161.          ENDIF
  162.       ENDIF
  163. ***   verification des deux conditions à la fois
  164.       IF ((ITEST1.EQ.1).AND.(ITEST2.EQ.1)) ITEST=1
  165.  
  166.       RETURN
  167.       END
  168.  
  169.  
  170.  
  171.  
  172.  
  173.  
  174.  
  175.  
  176.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales