1 : C FACTRA SOURCE CHAT 05/01/12 23:56:02 5004 2 : SUBROUTINE FACTRA(NDIM,ITYG,XDEP2,XARI2,IZSH,IFACE,XINT,XN,NBFAC) 3 : ************************************************************************* 4 : *** SP 'FACTRA' : par convection ou diffusion explicite (où traj=droite), 5 : *** donne n° faces traversees, pts intersections associés, normales du 6 : *** repere local à la face associées. 7 : *** 8 : *** APPELES 1 = aucun 9 : *** APPELES 2 = 'FACTR1', 'FACTR2' 10 : *** 11 : *** E = 'NDIM' dimension de l'espace 12 : *** 'ITYG' entier caracterisant la geometrie de l'element 13 : *** 'XDEP2' coordonnees reelles de depart de la particule 14 : *** 'XARI2' coordonnees reelles d'arrivee de la particule 15 : *** 'IZSH' segment content coords reelles des noeuds de l'elemt considere 16 : *** 17 : *** S = 'IFACE' n° des faces traversees par particule 18 : *** 'XINT' pts intersection trajectoire particule-faces traversees 19 : *** 'XN' vecteurs unitaires normaux aux faces traversees 20 : *** 'NBFAC' nbre de faces de l'element considere traversees par particule 21 : *** 22 : *** Auteur Cyril Nou 23 : ************************************************************************* 24 : 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8 (A-H,O-Z) 27 : SEGMENT IZSH 28 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 29 : ENDSEGMENT 30 : DIMENSION XDEP2(3),XARI2(3),IFACE(6),XINT(3,6),XN(3,6) 31 : DIMENSION XN1(3),XN2(3),XN3(3),XINTER(3) 32 : DIMENSION PT1(3),PT2(3),PT3(3),PT4(3) 33 : DIMENSION ITRIAN(4),ICARRE(5),ICUBE(4,6),IPRISM(4,5),ITETRA(4,4) 34 : *** donnees ordonnées specifiquement pour chaque type d'element afin de 35 : *** parcourir les faces dans l'ordre croissant dans les diverses boucles 36 : DATA ITRIAN/1,2,3,1/ 37 : DATA ICARRE/1,2,3,4,1/ 38 : DATA ICUBE/1,2,3,4, 5,6,7,8, 1,2,6,5, 2,3,7,6, 3,4,8,7, 1,4,8,5/ 39 : DATA IPRISM/1,2,3,3, 4,5,6,6, 1,2,5,4, 2,3,6,5, 1,3,6,4/ 40 : DATA ITETRA/1,2,4,4, 1,2,3,3, 2,3,4,4, 1,3,4,4/ 41 : *** initialisation des arguments de sortie à 0 42 : DO 10 I=1,6 43 : IFACE(I)=0 44 : DO 20 J=1,NDIM 45 : XINT(J,I)=0.D0 46 : XN(J,I)=0.D0 47 : 20 CONTINUE 48 : 10 CONTINUE 49 : NBFAC=0 50 : 51 : *** cas TRI3 (triangles), boucle sur les différentes faces 52 : IF (ITYG.EQ.4) THEN 53 : DO 30 I=1,3 54 : *** recuperation des pts definissant la face 55 : DO 35 J=1,NDIM 56 : PT1(J)=XYZL(J,ITRIAN(I)) 57 : PT2(J)=XYZL(J,ITRIAN(I+1)) 58 : 35 CONTINUE 59 : *** recherche des normales et du pt d'intersection 59 : CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 60 : $ ,XINTER,ITEST) 61 : *** 'ITEST'=1 s'il y a intersection avec plan associé à la face 61 : IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 62 : $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) 63 : *** 'ITEST'=1 s'il y a intersection avec la face 64 : IF (ITEST.EQ.1) THEN 65 : *** recuperation des infos de la face effectivement traversée 66 : NBFAC=NBFAC+1 67 : IFACE(NBFAC)=I 68 : DO 40 K=1,NDIM 69 : XINT(K,NBFAC)=XINTER(K) 70 : XN(K,NBFAC)=XN1(K) 71 : 40 CONTINUE 72 : ENDIF 73 : 30 CONTINUE 74 : *** cas QUA4 (quadrangles) 75 : ELSEIF (ITYG.EQ.8) THEN 76 : DO 50 I=1,4 77 : DO 55 J=1,NDIM 78 : PT1(J)=XYZL(J,ICARRE(I)) 79 : PT2(J)=XYZL(J,ICARRE(I+1)) 80 : 55 CONTINUE 80 : CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 81 : $ ,XINTER,ITEST) 81 : IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 82 : $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) 83 : IF (ITEST.EQ.1) THEN 84 : NBFAC=NBFAC+1 85 : IFACE(NBFAC)=I 86 : DO 60 K=1,NDIM 87 : XINT(K,NBFAC)=XINTER(K) 88 : XN(K,NBFAC)=XN1(K) 89 : 60 CONTINUE 90 : ENDIF 91 : 50 CONTINUE 92 : *** cas CUB8 (cubes) 93 : ELSEIF (ITYG.EQ.14) THEN 94 : DO 70 I=1,6 95 : DO 75 J=1,NDIM 96 : PT1(J)=XYZL(J,ICUBE(1,I)) 97 : PT2(J)=XYZL(J,ICUBE(2,I)) 98 : PT3(J)=XYZL(J,ICUBE(3,I)) 99 : PT4(J)=XYZL(J,ICUBE(4,I)) 100 : 75 CONTINUE 100 : CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 101 : $ ,XINTER,ITEST) 101 : IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 102 : $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) 103 : IF (ITEST.EQ.1) THEN 104 : NBFAC=NBFAC+1 105 : IFACE(NBFAC)=I 106 : DO 80 K=1,NDIM 107 : XINT(K,NBFAC)=XINTER(K) 108 : XN(K,NBFAC)=XN1(K) 109 : 80 CONTINUE 110 : ENDIF 111 : 70 CONTINUE 112 : *** cas PRI6 (prismes) 113 : ELSEIF (ITYG.EQ.16) THEN 114 : DO 90 I=1,5 115 : DO 95 J=1,NDIM 116 : PT1(J)=XYZL(J,IPRISM(1,I)) 117 : PT2(J)=XYZL(J,IPRISM(2,I)) 118 : PT3(J)=XYZL(J,IPRISM(3,I)) 119 : PT4(J)=XYZL(J,IPRISM(4,I)) 120 : 95 CONTINUE 120 : CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 121 : $ ,XINTER,ITEST) 121 : IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 122 : $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) 123 : IF (ITEST.EQ.1) THEN 124 : NBFAC=NBFAC+1 125 : IFACE(NBFAC)=I 126 : DO 100 K=1,NDIM 127 : XINT(K,NBFAC)=XINTER(K) 128 : XN(K,NBFAC)=XN1(K) 129 : 100 CONTINUE 130 : ENDIF 131 : 90 CONTINUE 132 : *** cas TET4 (tetraedres) 133 : ELSEIF (ITYG.EQ.23) THEN 134 : DO 110 I=1,4 135 : DO 115 J=1,NDIM 136 : PT1(J)=XYZL(J,ITETRA(1,I)) 137 : PT2(J)=XYZL(J,ITETRA(2,I)) 138 : PT3(J)=XYZL(J,ITETRA(3,I)) 139 : PT4(J)=XYZL(J,ITETRA(4,I)) 140 : 115 CONTINUE 140 : CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 141 : $ ,XINTER,ITEST) 141 : IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 142 : $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) 143 : IF (ITEST.EQ.1) THEN 144 : NBFAC=1 145 : IFACE(NBFAC)=I 146 : DO 120 K=1,NDIM 147 : XINT(K,NBFAC)=XINTER(K) 148 : XN(K,NBFAC)=XN1(K) 149 : 120 CONTINUE 150 : ENDIF 151 : 110 CONTINUE 152 : ENDIF 153 : RETURN 154 : END 155 : 156 : 157 : 158 : 159 : 160 : 161 : 162 :
© Cast3M 2003 - All rights reserved.
Disclaimer