C FACTRA SOURCE CHAT 05/01/12 23:56:02 5004 SUBROUTINE FACTRA(NDIM,ITYG,XDEP2,XARI2,IZSH,IFACE,XINT,XN,NBFAC) ************************************************************************* *** SP 'FACTRA' : par convection ou diffusion explicite (où traj=droite), *** donne n° faces traversees, pts intersections associés, normales du *** repere local à la face associées. *** *** APPELES 1 = aucun *** APPELES 2 = 'FACTR1', 'FACTR2' *** *** E = 'NDIM' dimension de l'espace *** 'ITYG' entier caracterisant la geometrie de l'element *** 'XDEP2' coordonnees reelles de depart de la particule *** 'XARI2' coordonnees reelles d'arrivee de la particule *** 'IZSH' segment content coords reelles des noeuds de l'elemt considere *** *** S = 'IFACE' n° des faces traversees par particule *** 'XINT' pts intersection trajectoire particule-faces traversees *** 'XN' vecteurs unitaires normaux aux faces traversees *** 'NBFAC' nbre de faces de l'element considere traversees par particule *** *** Auteur Cyril Nou ************************************************************************* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) SEGMENT IZSH REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) ENDSEGMENT DIMENSION XDEP2(3),XARI2(3),IFACE(6),XINT(3,6),XN(3,6) DIMENSION XN1(3),XN2(3),XN3(3),XINTER(3) DIMENSION PT1(3),PT2(3),PT3(3),PT4(3) DIMENSION ITRIAN(4),ICARRE(5),ICUBE(4,6),IPRISM(4,5),ITETRA(4,4) *** donnees ordonnées specifiquement pour chaque type d'element afin de *** parcourir les faces dans l'ordre croissant dans les diverses boucles DATA ITRIAN/1,2,3,1/ DATA ICARRE/1,2,3,4,1/ 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/ DATA IPRISM/1,2,3,3, 4,5,6,6, 1,2,5,4, 2,3,6,5, 1,3,6,4/ DATA ITETRA/1,2,4,4, 1,2,3,3, 2,3,4,4, 1,3,4,4/ *** initialisation des arguments de sortie à 0 DO 10 I=1,6 IFACE(I)=0 DO 20 J=1,NDIM XINT(J,I)=0.D0 XN(J,I)=0.D0 20 CONTINUE 10 CONTINUE NBFAC=0 *** cas TRI3 (triangles), boucle sur les différentes faces IF (ITYG.EQ.4) THEN DO 30 I=1,3 *** recuperation des pts definissant la face DO 35 J=1,NDIM PT1(J)=XYZL(J,ITRIAN(I)) PT2(J)=XYZL(J,ITRIAN(I+1)) 35 CONTINUE *** recherche des normales et du pt d'intersection CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 $ ,XINTER,ITEST) *** 'ITEST'=1 s'il y a intersection avec plan associé à la face IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) *** 'ITEST'=1 s'il y a intersection avec la face IF (ITEST.EQ.1) THEN *** recuperation des infos de la face effectivement traversée NBFAC=NBFAC+1 IFACE(NBFAC)=I DO 40 K=1,NDIM XINT(K,NBFAC)=XINTER(K) XN(K,NBFAC)=XN1(K) 40 CONTINUE ENDIF 30 CONTINUE *** cas QUA4 (quadrangles) ELSEIF (ITYG.EQ.8) THEN DO 50 I=1,4 DO 55 J=1,NDIM PT1(J)=XYZL(J,ICARRE(I)) PT2(J)=XYZL(J,ICARRE(I+1)) 55 CONTINUE CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 $ ,XINTER,ITEST) IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) IF (ITEST.EQ.1) THEN NBFAC=NBFAC+1 IFACE(NBFAC)=I DO 60 K=1,NDIM XINT(K,NBFAC)=XINTER(K) XN(K,NBFAC)=XN1(K) 60 CONTINUE ENDIF 50 CONTINUE *** cas CUB8 (cubes) ELSEIF (ITYG.EQ.14) THEN DO 70 I=1,6 DO 75 J=1,NDIM PT1(J)=XYZL(J,ICUBE(1,I)) PT2(J)=XYZL(J,ICUBE(2,I)) PT3(J)=XYZL(J,ICUBE(3,I)) PT4(J)=XYZL(J,ICUBE(4,I)) 75 CONTINUE CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 $ ,XINTER,ITEST) IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) IF (ITEST.EQ.1) THEN NBFAC=NBFAC+1 IFACE(NBFAC)=I DO 80 K=1,NDIM XINT(K,NBFAC)=XINTER(K) XN(K,NBFAC)=XN1(K) 80 CONTINUE ENDIF 70 CONTINUE *** cas PRI6 (prismes) ELSEIF (ITYG.EQ.16) THEN DO 90 I=1,5 DO 95 J=1,NDIM PT1(J)=XYZL(J,IPRISM(1,I)) PT2(J)=XYZL(J,IPRISM(2,I)) PT3(J)=XYZL(J,IPRISM(3,I)) PT4(J)=XYZL(J,IPRISM(4,I)) 95 CONTINUE CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 $ ,XINTER,ITEST) IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) IF (ITEST.EQ.1) THEN NBFAC=NBFAC+1 IFACE(NBFAC)=I DO 100 K=1,NDIM XINT(K,NBFAC)=XINTER(K) XN(K,NBFAC)=XN1(K) 100 CONTINUE ENDIF 90 CONTINUE *** cas TET4 (tetraedres) ELSEIF (ITYG.EQ.23) THEN DO 110 I=1,4 DO 115 J=1,NDIM PT1(J)=XYZL(J,ITETRA(1,I)) PT2(J)=XYZL(J,ITETRA(2,I)) PT3(J)=XYZL(J,ITETRA(3,I)) PT4(J)=XYZL(J,ITETRA(4,I)) 115 CONTINUE CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3 $ ,XINTER,ITEST) IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4 $ ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST) IF (ITEST.EQ.1) THEN NBFAC=1 IFACE(NBFAC)=I DO 120 K=1,NDIM XINT(K,NBFAC)=XINTER(K) XN(K,NBFAC)=XN1(K) 120 CONTINUE ENDIF 110 CONTINUE ENDIF RETURN END