Télécharger la plaquette Cast3M


fr en

Télécharger itrc.eso

Retour à la liste

Numérotation des lignes : 

itrc

  1. C ITRC      SOURCE    FANDEUR   22/01/03    21:15:22     11136          
  2.       SUBROUTINE ITRC
  3. CA   Ph. MIRANDA
  4. CD   30/11/2001
  5. CC--------------------------------------------------------------------------
  6. CC
  7. CC     Le sous-programme ITRC calcule par une methode analytique exacte
  8. CC     l'intersection entre des segments (connus par les coordonnees de
  9. CC     leurs extremites) et les facettes triangulaires d'un maillage
  10. CC     (connues par les coordonnees de leurs 3 extremites et leurs
  11. CC     normales).
  12. CC
  13. CC     ITRC reprend les fonctionnalites de la procedure GIBIANE @INTERC
  14. CC     concue par R.Mitteau (CEA/DSM/DRFC) et B.Riou (CSSI) en oct 1998.
  15. CC
  16. CC     version initiale : Ph.Miranda (CSSI)    en nov 2001
  17. CC     modification     : T.Charras  (CEA/DM2S) en oct 2003
  18. CC     modification     : A.Moal     (CSSI)    en jan 2004
  19. CC
  20. CC     Syntaxe generale :
  21. CC
  22. CC     CHDIST MINTER = ITRC CHOLD CHNEW TOL TAB1 ;
  23. CC
  24. CC--------------------------------------------------------------------------
  25. CE     En entree :
  26. CE
  27. CE     CHOLD : coordonnees des points extremites initiales des segments
  28. CE            (type CHPOINT a 3 composantes)
  29. CE     CHNEW : coordonnees des points extremites finales des segments
  30. CE            (type CHPOINT a 3 composantes appuye sur le meme support
  31. CE             que CH_OLD)
  32. CE     TOL   : tolerance (type FLOTTANT)
  33. CE     TAB1  : Table qui doit contenir les parametres suivants en entree
  34. CE
  35. CE      <chamx1 (type MCHAML) coordonnee x du premier noeud de l'element
  36. CE      <chamy1 (type MCHAML) coordonnee y du premier noeud de l'element
  37. CE      <chamz1 (type MCHAML) coordonnee z du premier noeud de l'element
  38. CE
  39. CE      <chamx2 (type MCHAML) coordonnee x du deuxieme noeud de l'element
  40. CE      <chamy2 (type MCHAML) coordonnee y du deuxieme noeud de l'element
  41. CE      <chamz2 (type MCHAML) coordonnee z du deuxieme noeud de l'element
  42. CE
  43. CE      <chamx3 (type MCHAML) coordonnee x du troisieme noeud de l'element
  44. CE      <chamy3 (type MCHAML) coordonnee y du troisieme noeud de l'element
  45. CE      <chamz3 (type MCHAML) coordonnee z du troisieme noeud de l'element
  46. CE
  47. CE      <cosx (type MCHAML) cosinus directeur selon x
  48. CE      <cosy (type MCHAML) cosinus directeur selon y
  49. CE      <cosz (type MCHAML) cosinus directeur selon z
  50. CC
  51. CS      En sortie :
  52. CS
  53. CS      CHDIST : CHPOINT appuye sur le support de CHOLD et contenant
  54. CS               le pas si pas d'intersection et la distance entre
  55. CS               point initial du segment et le point d'intersection
  56. CS               a la facette sinon (type CHPOINT)
  57. CS      MINTER : contient les noeuds intersectes du support de CHOLD.
  58. CS               (type MAILLAGE)
  59. CS
  60. CC      Remarque : il est necessaire de travailler avec les champs par
  61. CC                 element du maillage interceptant et non pas avec les
  62. CC                 coordonnées des points de l'objet maillage car ces
  63. CC                 champs peuvent avoir subi des changements de repere
  64. CC
  65. CC--------------------------------------------------------------------------
  66.       IMPLICIT INTEGER(I-N)
  67.       IMPLICIT REAL*8 (A-H,O-Z)
  68. C
  69.  
  70. -INC PPARAM
  71. -INC CCOPTIO
  72. -INC SMCHPOI
  73. -INC SMCHAML
  74. -INC SMELEME
  75. -INC SMCOORD
  76. -INC CCGEOME
  77. -INC CCREEL
  78. -INC SMTABLE
  79. C
  80.       SEGMENT MTRAV1
  81.         REAL*8 XPN(NZ1,NBP),YPN(NZ1,NBP),ZPN(NZ1,NBP)
  82.         REAL*8 XPN1(NZ1,NBP),YPN1(NZ1,NBP),ZPN1(NZ1,NBP)
  83.         REAL*8 CHPAS(NZ1,NBP),CHPAS0(NZ1,NBP)
  84.         LOGICAL INTER(NZ1,NBP)
  85.       ENDSEGMENT
  86. C
  87.       SEGMENT MTRAV2
  88.         REAL*8 COORFAC(NZ2,NEL,3,3)
  89.         REAL*8 TXN(NZ2,NEL),TYN(NZ2,NEL),TZN(NZ2,NEL)
  90.       ENDSEGMENT
  91. C
  92.       POINTEUR CHOLD.MCHPOI,CHNEW.MCHPOI
  93.       POINTEUR CHX1.MCHPOI,CHX2.MCHPOI,CHX3.MCHPOI
  94.       POINTEUR CHY1.MCHPOI,CHY2.MCHPOI,CHY3.MCHPOI
  95.       POINTEUR CHZ1.MCHPOI,CHZ2.MCHPOI,CHZ3.MCHPOI
  96.       POINTEUR TAB1.MTABLE
  97.       POINTEUR ICHOLD.MELEME
  98.       POINTEUR CHCOSX.MCHAML,CHCOSY.MCHAML,CHCOSZ.MCHAML
  99. C
  100.       CHARACTER*8 MTYPR
  101.       CHARACTER*4 REP
  102.       CHARACTER*9 MOTCLE(12)
  103.       logical login,logre
  104.       character*4 charre
  105. C
  106. C  LECTURE DES OBJETS
  107. C
  108.       CALL LIROBJ('CHPOINT',CHOLD,1,IRETOU)
  109.       IF (IERR .NE. 0) RETURN
  110.       CALL LIROBJ('CHPOINT',CHNEW,1,IRETOU)
  111.       IF (IERR .NE. 0) RETURN
  112.       CALL LIRREE(TOL1,1,IRETOU)
  113.       IF (IERR .NE. 0) RETURN
  114.       CALL LIROBJ('TABLE   ',TAB1,1,IRETOU)
  115.       IF (IERR .NE. 0) RETURN
  116. C
  117. C  Dimensionnement des tableaux du segment MTRAV1
  118. C
  119.       MCHPOI=CHOLD
  120.       SEGACT MCHPOI
  121.       NZ1=IPCHP(/1)
  122.       DO 50 I=1,IPCHP(/1)
  123.          MSOUPO=IPCHP(I)
  124.          SEGACT MSOUPO
  125.          MELEME=IGEOC
  126.          SEGACT MELEME
  127.          IF (I .EQ. 1) NBP=NUM(/2)
  128.          IF (NUM(/2) .GT. NBP) NBP=NUM(/2)
  129.          SEGDES MELEME
  130.          SEGDES MSOUPO
  131. 50    CONTINUE
  132.       SEGINI MTRAV1
  133. C
  134. C
  135. C   EXTRACTION DES COORDONNEES DES POINTS Pn et Pn+1 a partir des CHPOINT
  136. C       CHOLD et CHNEW
  137. C
  138.       DO 100 I=1,IPCHP(/1)
  139.          MSOUPO=IPCHP(I)
  140.          SEGACT MSOUPO
  141.          ICHOLD=IGEOC
  142.          MPOVAL=IPOVAL
  143.          SEGACT MPOVAL
  144.          DO 200 J=1,NOCOMP(/2)
  145.              IF (NOCOMP(J).EQ. 'X   ')then
  146.               do 331 k=1,vpocha(/1)
  147.                 XPN(I,K)=VPOCHA(K,J)
  148.   331         continue
  149.              elseif(NOCOMP(J).EQ. 'Y   ') then
  150.               do 332 k=1,vpocha(/1)
  151.                 YPN(I,K)=VPOCHA(K,J)
  152.   332         continue
  153.              elseif(NOCOMP(J).EQ. 'Z   ')then
  154.               do 333 k=1,vpocha(/1)
  155.                 ZPN(I,K)=VPOCHA(K,J)
  156.   333         continue
  157.              endif
  158. 200      CONTINUE
  159.          SEGDES MPOVAL
  160.          SEGDES MELEME
  161.          SEGDES MSOUPO
  162. 100   CONTINUE
  163.       SEGDES MCHPOI
  164. C
  165.       MCHPOI=CHNEW
  166.       SEGACT MCHPOI
  167.       DO 110 I=1,IPCHP(/1)
  168.          MSOUPO=IPCHP(I)
  169.          SEGACT MSOUPO
  170.          MPOVAL=IPOVAL
  171.          SEGACT MPOVAL
  172.          DO 210 J=1,NOCOMP(/2)
  173.                IF (NOCOMP(J).EQ. 'X   ')then
  174.                    DO 310 K=1,VPOCHA(/1)
  175.                     XPN1(I,K)=VPOCHA(K,J)
  176.   310              CONTINUE
  177.                ELSEIF (NOCOMP(J).EQ. 'Y   ')then
  178.                    DO 311 K=1,VPOCHA(/1)
  179.                     YPN1(I,K)=VPOCHA(K,J)
  180.   311              CONTINUE
  181.                ELSEIF (NOCOMP(J).EQ. 'Z   ')then
  182.                    DO 312 K=1,VPOCHA(/1)
  183.                     ZPN1(I,K)=VPOCHA(K,J)
  184.   312              CONTINUE
  185.                ENDIF
  186. 210      CONTINUE
  187.          SEGDES MPOVAL
  188.          SEGDES MSOUPO
  189. 110   CONTINUE
  190.       SEGDES MCHPOI
  191. C
  192.          MOTCLE(1) = '<COSX    '
  193.          MOTCLE(2) = '<COSY    '
  194.          MOTCLE(3) = '<COSZ    '
  195.          MOTCLE(4) = '<CHAMX1  '
  196.          MOTCLE(5) = '<CHAMY1  '
  197.          MOTCLE(6) = '<CHAMZ1  '
  198.          MOTCLE(7) = '<CHAMX2  '
  199.          MOTCLE(8) = '<CHAMY2  '
  200.          MOTCLE(9) = '<CHAMZ2  '
  201.          MOTCLE(10) = '<CHAMX3  '
  202.          MOTCLE(11) = '<CHAMY3  '
  203.          MOTCLE(12) = '<CHAMZ3  '
  204. C
  205.       SEGACT TAB1
  206.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(1),LOGIN,IOBIN,
  207.      .                 'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHCOSX)
  208.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(2),LOGIN,IOBIN,
  209.      .                 'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHCOSY)
  210.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(3),LOGIN,IOBIN,
  211.      .                 'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHCOSZ)
  212. C
  213.       MCHELM=CHCOSX
  214.       SEGACT MCHELM
  215. C
  216. C  Dimensionnement des tableaux du segment MTRAV2
  217. C
  218.       NZ2 = ICHAML(/1)
  219.       DO 60 I=1,ICHAML(/1)
  220.          MCHAML=ICHAML(I)
  221.          SEGACT MCHAML
  222.          DO 70 IC=1,IELVAL(/1)
  223.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  224.                 MELVAL=IELVAL(IC)
  225.                 SEGACT MELVAL
  226.                 IF (I .EQ. 1) NEL=VELCHE(/2)
  227.                 IF (VELCHE(/2) .GT. NEL) NEL=VELCHE(/2)
  228.                 SEGDES MELVAL
  229.              ENDIF
  230. 70      CONTINUE
  231.         SEGDES MCHAML
  232. 60    CONTINUE
  233.       SEGINI MTRAV2
  234. C
  235. C  Extraction des trois cosinus directeurs (CHAMELEM) a partir de TAB1
  236. C  Les elements du maillage associe sont ici des POI1
  237. C
  238.       DO 120 I=1,ICHAML(/1)
  239.          MCHAML=ICHAML(I)
  240.          SEGACT MCHAML
  241.          DO 220 IC=1,IELVAL(/1)
  242.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  243.                 MELVAL=IELVAL(IC)
  244.                 SEGACT MELVAL
  245.                 DO 320 IEL=1,VELCHE(/2)
  246.                    TXN(I,IEL) = VELCHE(1,IEL)
  247. 320             CONTINUE
  248.                 SEGDES MELVAL
  249.              ENDIF
  250. 220      CONTINUE
  251.          SEGDES MCHAML
  252. 120   CONTINUE
  253.       SEGDES MCHELM
  254. C
  255.       MCHELM=CHCOSY
  256.       SEGACT MCHELM
  257.       DO 130 I=1,ICHAML(/1)
  258.          MCHAML=ICHAML(I)
  259.          SEGACT MCHAML
  260.          DO 230 IC=1,IELVAL(/1)
  261.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  262.                 MELVAL=IELVAL(IC)
  263.                 SEGACT MELVAL
  264.                 DO 330 IEL=1,VELCHE(/2)
  265.                    TYN(I,IEL) = VELCHE(1,IEL)
  266. 330             CONTINUE
  267.                 SEGDES MELVAL
  268.              ENDIF
  269. 230      CONTINUE
  270.          SEGDES MCHAML
  271. 130   CONTINUE
  272.       SEGDES MCHELM
  273. C
  274.       MCHELM=CHCOSZ
  275.       SEGACT MCHELM
  276.       DO 140 I=1,ICHAML(/1)
  277.          MCHAML=ICHAML(I)
  278.          SEGACT MCHAML
  279.          DO 240 IC=1,IELVAL(/1)
  280.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  281.                 MELVAL=IELVAL(IC)
  282.                 SEGACT MELVAL
  283.                 DO 340 IEL=1,VELCHE(/2)
  284.                    TZN(I,IEL) = VELCHE(1,IEL)
  285. 340             CONTINUE
  286.                 SEGDES MELVAL
  287.              ENDIF
  288. 240      CONTINUE
  289.          SEGDES MCHAML
  290. 140   CONTINUE
  291.       SEGDES MCHELM
  292. C
  293. C   Extraction des coordonnees des 3 points de chaque facette
  294. C
  295.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(4),LOGIN,
  296.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHX1)
  297.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(5),LOGIN,
  298.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHY1)
  299.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(6),LOGIN,
  300.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHZ1)
  301. C
  302.       MCHELM=CHX1
  303.       SEGACT MCHELM
  304.       DO 141 I=1,ICHAML(/1)
  305.          MCHAML=ICHAML(I)
  306.          SEGACT MCHAML
  307.          DO 241 IC=1,IELVAL(/1)
  308.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  309.                 MELVAL=IELVAL(IC)
  310.                 SEGACT MELVAL
  311.                 DO 341 IEL=1,VELCHE(/2)
  312.                    COORFAC(I,IEL,1,1) = VELCHE(1,IEL)
  313. 341             CONTINUE
  314.                 SEGDES MELVAL
  315.              ENDIF
  316. 241      CONTINUE
  317.          SEGDES MCHAML
  318. 141   CONTINUE
  319.       SEGDES MCHELM
  320. C
  321.       MCHELM=CHY1
  322.       SEGACT MCHELM
  323.       DO 142 I=1,ICHAML(/1)
  324.          MCHAML=ICHAML(I)
  325.          SEGACT MCHAML
  326.          DO 242 IC=1,IELVAL(/1)
  327.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  328.                 MELVAL=IELVAL(IC)
  329.                 SEGACT MELVAL
  330.                 DO 342 IEL=1,VELCHE(/2)
  331.                    COORFAC(I,IEL,1,2) = VELCHE(1,IEL)
  332. 342             CONTINUE
  333.                 SEGDES MELVAL
  334.              ENDIF
  335. 242      CONTINUE
  336.          SEGDES MCHAML
  337. 142   CONTINUE
  338.       SEGDES MCHELM
  339. C
  340.       MCHELM=CHZ1
  341.       SEGACT MCHELM
  342.       DO 143 I=1,ICHAML(/1)
  343.          MCHAML=ICHAML(I)
  344.          SEGACT MCHAML
  345.          DO 243 IC=1,IELVAL(/1)
  346.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  347.                 MELVAL=IELVAL(IC)
  348.                 SEGACT MELVAL
  349.                 DO 343 IEL=1,VELCHE(/2)
  350.                    COORFAC(I,IEL,1,3) = VELCHE(1,IEL)
  351. 343             CONTINUE
  352.                 SEGDES MELVAL
  353.              ENDIF
  354. 243      CONTINUE
  355.          SEGDES MCHAML
  356. 143   CONTINUE
  357.       SEGDES MCHELM
  358. C
  359.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(7),LOGIN,
  360.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHX2)
  361.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(8),LOGIN,
  362.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHY2)
  363.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(9),LOGIN,
  364.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHZ2)
  365. C
  366.       MCHELM=CHX2
  367.       SEGACT MCHELM
  368.       DO 144 I=1,ICHAML(/1)
  369.          MCHAML=ICHAML(I)
  370.          SEGACT MCHAML
  371.          DO 244 IC=1,IELVAL(/1)
  372.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  373.                 MELVAL=IELVAL(IC)
  374.                 SEGACT MELVAL
  375.                 DO 344 IEL=1,VELCHE(/2)
  376.                     COORFAC(I,IEL,2,1) = VELCHE(1,IEL)
  377. 344             CONTINUE
  378.                 SEGDES MELVAL
  379.              ENDIF
  380. 244      CONTINUE
  381.          SEGDES MCHAML
  382. 144   CONTINUE
  383.       SEGDES MCHELM
  384. C
  385.       MCHELM=CHY2
  386.       SEGACT MCHELM
  387.       DO 145 I=1,ICHAML(/1)
  388.          MCHAML=ICHAML(I)
  389.          SEGACT MCHAML
  390.          DO 245 IC=1,IELVAL(/1)
  391.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  392.                 MELVAL=IELVAL(IC)
  393.                 SEGACT MELVAL
  394.                 DO 345 IEL=1,VELCHE(/2)
  395.                    COORFAC(I,IEL,2,2) = VELCHE(1,IEL)
  396. 345             CONTINUE
  397.                 SEGDES MELVAL
  398.              ENDIF
  399. 245      CONTINUE
  400.          SEGDES MCHAML
  401. 145   CONTINUE
  402.       SEGDES MCHELM
  403. C
  404.       MCHELM=CHZ2
  405.       SEGACT MCHELM
  406.       DO 146 I=1,ICHAML(/1)
  407.          MCHAML=ICHAML(I)
  408.          SEGACT MCHAML
  409.          DO 246 IC=1,IELVAL(/1)
  410.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  411.                 MELVAL=IELVAL(IC)
  412.                 SEGACT MELVAL
  413.                 DO 346 IEL=1,VELCHE(/2)
  414.                    COORFAC(I,IEL,2,3) = VELCHE(1,IEL)
  415. 346             CONTINUE
  416.                 SEGDES MELVAL
  417.              ENDIF
  418. 246      CONTINUE
  419.          SEGDES MCHAML
  420. 146   CONTINUE
  421.       SEGDES MCHELM
  422. C
  423.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(10),LOGIN,
  424.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHX3)
  425.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(11),LOGIN,
  426.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHY3)
  427.       CALL ACCTAB(TAB1,'MOT     ',IVALIN,XVALIN,MOTCLE(12),LOGIN,
  428.      .           IOBIN,'MCHAML  ',IVALRE,XVALRE,CHARRE,LOGRE,CHZ3)
  429. C
  430.       MCHELM=CHX3
  431.       SEGACT MCHELM
  432.       DO 147 I=1,ICHAML(/1)
  433.          MCHAML=ICHAML(I)
  434.          SEGACT MCHAML
  435.          DO 247 IC=1,IELVAL(/1)
  436.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  437.                 MELVAL=IELVAL(IC)
  438.                 SEGACT MELVAL
  439.                 DO 347 IEL=1,VELCHE(/2)
  440.                    COORFAC(I,IEL,3,1) = VELCHE(1,IEL)
  441. 347             CONTINUE
  442.                 SEGDES MELVAL
  443.              ENDIF
  444. 247      CONTINUE
  445.          SEGDES MCHAML
  446. 147   CONTINUE
  447.       SEGDES MCHELM
  448. C
  449.       MCHELM=CHY3
  450.       SEGACT MCHELM
  451.       DO 148 I=1,ICHAML(/1)
  452.          MCHAML=ICHAML(I)
  453.          SEGACT MCHAML
  454.          DO 248 IC=1,IELVAL(/1)
  455.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  456.                 MELVAL=IELVAL(IC)
  457.                 SEGACT MELVAL
  458.                 DO 348 IEL=1,VELCHE(/2)
  459.                    COORFAC(I,IEL,3,2) = VELCHE(1,IEL)
  460. 348             CONTINUE
  461.                 SEGDES MELVAL
  462.              ENDIF
  463. 248      CONTINUE
  464.          SEGDES MCHAML
  465. 148   CONTINUE
  466.       SEGDES MCHELM
  467. C
  468.       MCHELM=CHZ3
  469.       SEGACT MCHELM
  470.       DO 149 I=1,ICHAML(/1)
  471.          MCHAML=ICHAML(I)
  472.          SEGACT MCHAML
  473.          DO 249 IC=1,IELVAL(/1)
  474.              IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  475.                 MELVAL=IELVAL(IC)
  476.                 SEGACT MELVAL
  477.                 DO 349 IEL=1,VELCHE(/2)
  478.                    COORFAC(I,IEL,3,3) = VELCHE(1,IEL)
  479. 349             CONTINUE
  480.                 SEGDES MELVAL
  481.              ENDIF
  482. 249      CONTINUE
  483.          SEGDES MCHAML
  484. 149   CONTINUE
  485.       SEGDES MCHELM
  486. C
  487.       SEGDES TAB1
  488. C
  489. C  Initialisation du tableau des indicateurs d'interception
  490. C
  491.       DO 302 I=1,NZ1
  492.          DO 301 K=1,NBP
  493.             INTER(I,K)=.FALSE.
  494. 301      CONTINUE
  495. 302   CONTINUE
  496. C
  497. C   Boucle sur les facettes du maillage OMBRANT
  498. C
  499.       MCHELM=CHX1
  500.       SEGACT MCHELM
  501.       DO 400 I=1,ICHAML(/1)
  502.          MCHAML=ICHAML(I)
  503.          SEGACT MCHAML
  504.          DO 450 IC=1,IELVAL(/1)
  505.             IF (TYPCHE(IC)(1:6) .EQ. 'REAL*8') THEN
  506.                MELVAL=IELVAL(IC)
  507.                SEGACT MELVAL
  508.                DO 500 J=1,VELCHE(/2)
  509.                   XA1 =  COORFAC(I,J,1,1)
  510.                   YA1 =  COORFAC(I,J,1,2)
  511.                   ZA1 =  COORFAC(I,J,1,3)
  512. C
  513.                   XB1 =  COORFAC(I,J,2,1)
  514.                   YB1 =  COORFAC(I,J,2,2)
  515.                   ZB1 =  COORFAC(I,J,2,3)
  516. C
  517.                   XC1 =  COORFAC(I,J,3,1)
  518.                   YC1 =  COORFAC(I,J,3,2)
  519.                   ZC1 =  COORFAC(I,J,3,3)
  520. C
  521. C  Boucle sur les sous-zones du maillage OMBRE
  522. C
  523.                   DO 550 IZ=1,NZ1
  524. C
  525. C  Boucle sur les lignes de champ (ou sur les elements POI1 de OMBRE)
  526. C
  527.                      DO 600 K=1,NBP
  528.                         IF (.NOT. INTER(IZ,K)) THEN
  529. C
  530. C  Produits scalaires
  531. C
  532.                            PS1 = ((XPN(IZ,K)-XA1)*TXN(I,J)) +
  533.      .                           ((YPN(IZ,K)-YA1)*TYN(I,J)) +
  534.      .                           ((ZPN(IZ,K)-ZA1)*TZN(I,J))
  535.                            PS2 = ((XPN1(IZ,K)-XA1)*TXN(I,J)) +
  536.      .                           ((YPN1(IZ,K)-YA1)*TYN(I,J)) +
  537.      .                           ((ZPN1(IZ,K)-ZA1)*TZN(I,J))
  538. C             -> -->
  539. C   Calcul de PnPn+1 :
  540. C
  541.                            DEX = XPN1(IZ,K) - XPN(IZ,K)
  542.                            DEY = YPN1(IZ,K) - YPN(IZ,K)
  543.                            DEZ = ZPN1(IZ,K) - ZPN(IZ,K)
  544. C
  545. C   calcul de la distance de connection par defaut
  546. C
  547.              CHPAS0(IZ,K) = (DEX*DEX + DEY*DEY + DEZ*DEZ)**.5
  548. C
  549. C   Test d'intersection avec le plan de la facette
  550. C
  551.                   IF ( (PS1*PS2 .LT. (-1*TOL1*TOL1))
  552.      .             .OR. (PS2 .EQ. 0.) )               THEN
  553. C
  554. C                   -> -->   ->
  555. C   Calcul de AD1 = PnPn+1 . N  :
  556. C
  557.                 AD1 = DEX*TXN(I,J) + DEY*TYN(I,J) + DEZ*TZN(I,J)
  558. C
  559. C   Test de detection des lignes de champ paralleles a la facette
  560. C
  561.                            IF (ABS(AD1) .GT. TOL1) THEN
  562. C
  563. C
  564. C    Calcul des points d'intersection M entre les lignes de champ et
  565. C                                             le plan de la facette
  566. C
  567.                            XM = ( (XA1*DEX*TXN(I,J)) -
  568.      .                          ((((YPN(IZ,K)-YA1)*DEX)-
  569.      .                             (XPN(IZ,K)*DEY))*TYN(I,J)) -
  570.      .                          ((((ZPN(IZ,K)-ZA1)*DEX)-
  571.      .                       (XPN(IZ,K)*DEZ))*TZN(I,J)) ) / AD1
  572. C
  573.                            YM = ( (YA1*DEY*TYN(I,J)) -
  574.      .                          ((((ZPN(IZ,K)-ZA1)*DEY)-
  575.      .                             (YPN(IZ,K)*DEZ))*TZN(I,J)) -
  576.      .                          ((((XPN(IZ,K)-XA1)*DEY)-
  577.      .                       (YPN(IZ,K)*DEX))*TXN(I,J)) ) / AD1
  578. C
  579.                            ZM = ( (ZA1*DEZ*TZN(I,J)) -
  580.      .                          ((((XPN(IZ,K)-XA1)*DEZ)-
  581.      .                             (ZPN(IZ,K)*DEX))*TXN(I,J)) -
  582.      .                          ((((YPN(IZ,K)-YA1)*DEZ)-
  583.      .                       (ZPN(IZ,K)*DEY))*TYN(I,J)) ) / AD1
  584. C
  585. C   Calcul des vecteurs AM, BM, CM
  586. C
  587.                            XAM = XM - XA1
  588.                            XBM = XM - XB1
  589.                            XCM = XM - XC1
  590. C
  591. C
  592.                            YAM = YM - YA1
  593.                            YBM = YM - YB1
  594.                            YCM = YM - YC1
  595. C
  596.                            ZAM = ZM - ZA1
  597.                            ZBM = ZM - ZB1
  598.                            ZCM = ZM - ZC1
  599. C
  600. C   Denominateurs
  601. C
  602.                            DX = (YCM*ZBM)-(ZCM*YBM)+(ZAM*YBM)-
  603.      .                          (ZAM*YCM)-(YAM*ZBM)+(YAM*ZCM)
  604. C
  605.                            DY = (ZCM*XBM)-(XCM*ZBM)+(XAM*ZBM)-
  606.      .                          (XAM*ZCM)-(ZAM*XBM)+(ZAM*XCM)
  607. C
  608.                            DZ = (YCM*XBM)-(XCM*YBM)+(XAM*YBM)-
  609.      .                          (XAM*YCM)-(YAM*XBM)+(YAM*XCM)
  610. C
  611.                            IF (ABS(DZ) .GT. TOL1) THEN
  612.                                D1 = DZ
  613.                                A1 = YCM*XBM-XCM*YBM
  614.                                B1 = YAM*XCM-XAM*YCM
  615.                                C1 = YBM*XAM-XBM*YAM
  616.                            ELSEIF (ABS(DY) .GT. TOL1) THEN
  617.                                D1 = DY
  618.                                A1 = ZCM*XBM-XCM*ZBM
  619.                                B1 = ZAM*XCM-XAM*ZCM
  620.                                C1 = ZBM*XAM-XBM*ZAM
  621.                            ELSE
  622.                                D1 = DX
  623.                                A1 = YCM*ZBM-ZCM*YBM
  624.                                B1 = YAM*ZCM-ZAM*YCM
  625.                                C1 = YBM*ZAM-ZBM*YAM
  626.                            ENDIF
  627. C
  628. C    Coordonnees barycentriques
  629. C
  630.                            ALPHA = A1 / D1
  631.                            BETA =  B1 / D1
  632.                            GAMMA = C1 / D1
  633. C
  634.                            IF ( (ALPHA .GE. 0.) .AND.
  635.      .                          (BETA  .GE. 0.) .AND.
  636.      .                          (GAMMA .GE. 0.)      ) THEN
  637. C
  638.                                     INTER(IZ,K)=.TRUE.
  639. C
  640. C   Calcul de d(Pn,M) :
  641. C
  642.                          CHPAS(IZ,K) = ((XPN(IZ,K)-XM)**2. +
  643.      .          (YPN(IZ,K)-YM)**2. + (ZPN(IZ,K)-ZM)**2.)**.5
  644.                                  ENDIF
  645.                               ENDIF
  646.                            ENDIF
  647.                         ENDIF
  648. 600                  CONTINUE
  649. 550               CONTINUE
  650. 500            CONTINUE
  651.                SEGDES MELVAL
  652.             ENDIF
  653. 450      CONTINUE
  654.          SEGDES MCHAML
  655. 400   CONTINUE
  656.       SEGDES MCHELM
  657. C
  658. C--------------------------------------------------------
  659. C
  660. C  Creation du maillage MINTER des points d'intersection
  661. C
  662.       NBNN = 1
  663.       NBELEM = NBP
  664.       NBSOUS = 0
  665.       NBREF = 0
  666.       SEGINI MELEME
  667.       MINTER  = MELEME
  668.       ITYPEL = 1
  669. C
  670. C   Construction de MINTER
  671. C
  672.       IPT1 = ICHOLD
  673.       SEGACT IPT1
  674.       NINTER = 0
  675.       DO 701 I=1,NZ1
  676.          DO 700 K=1,NBP
  677.             IF (INTER(I,K)) THEN
  678.                NINTER = NINTER+1
  679.                NUM(1,NINTER) = IPT1.NUM(1,K)
  680.             ENDIF
  681. 700      CONTINUE
  682. 701   CONTINUE
  683.       NBELEM = NINTER
  684.       SEGADJ MELEME
  685.       SEGDES MELEME,IPT1
  686. C--------------------------------------------------------
  687. C
  688. C  Creation du CHPOINT CHDIST
  689. C
  690. C segment MCHPOI
  691.       NAT = 1
  692.       NSOUPO = NZ1
  693.       SEGINI MCHPOI
  694.       ICHDIS = MCHPOI
  695.       IFOPOI = 2
  696.       JATTRI(1) = 2
  697. C segment MSOUPO
  698.       NC = 1
  699.       DO 901 I=1,NZ1
  700.          SEGINI MSOUPO
  701.          IPCHP(I)=MSOUPO
  702.          NOCOMP(1)= 'SCAL'
  703. C segment MELEME : maillage associe a CHDIST
  704.          NBNN = 1
  705.          NBELEM = NBP
  706.          NBSOUS = 0
  707.          NBREF = 0
  708. C
  709. C   LE CHAMPOINT CHDIST S'APPUIE SUR LE MEME MAILLAGE QUE CHOLD
  710. C
  711.          IGEOC=ICHOLD
  712. C
  713. C   Construction de MPOVAL
  714. C
  715.          N = NBP
  716.          SEGINI MPOVAL
  717.          IPOVAL=MPOVAL
  718. C
  719. C   On met la valeur de la distance de connection dans le chpoint CHDIST
  720. C       Ce point n'existera plus dans le prochain maillage MINTER
  721. C
  722.          DO 900 K=1,NBP
  723.             IF (INTER(I,K)) THEN
  724.                 VPOCHA(K,1)=CHPAS(I,K)
  725.             ELSE
  726.                 VPOCHA(K,1)=CHPAS0(I,K)
  727.             ENDIF
  728. 900      CONTINUE
  729. C
  730.          SEGDES MPOVAL
  731.          SEGDES MSOUPO
  732. 901   CONTINUE
  733.       SEGDES MCHPOI
  734. C
  735. C--------------------------------------------------------
  736. C
  737. C   Ecriture des nouveaux objets en sortie
  738. C
  739.       CALL ECROBJ('MAILLAGE',MINTER)
  740. C
  741.       CALL ECROBJ('CHPOINT ',ICHDIS)
  742. C
  743.       SEGSUP MTRAV1
  744.       SEGSUP MTRAV2
  745. C
  746.       RETURN
  747.       END
  748.  
  749.  
  750.  

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