C IJET SOURCE PV 09/03/12 21:25:10 6325 C IJET SOURCE PM 01/12/11 21:44:48 4041 SUBROUTINE IJET CA Ph. MIRANDA CD 30/11/2001 CC-------------------------------------------------------------------------- CC CC Le sous-programme IJET calcule par une methode analytique exacte CC l'intersection entre des segments (connus par les coordonnees de CC leurs extremites) et les facettes triangulaires d'un maillage CC (connues par les coordonnees de leurs 3 extremites et leurs CC normales). De plus, ce programme permet de deduire les valeurs CC d'un champ appuye sur les noeuds du maillage aux points CC d'intersection des segments par interpolation. CC CC IJET reprend les fonctionnalites de la procedure GIBIANE @INTSEC CC inspiree de @INTERC concue par R.Mitteau (CEA/DSM/DRFC) et CC B.Riou (CSSI) en oct 1998. CC CC version initiale : Ph.Miranda (CSSI) en nov 2001 CC modification : T.Charras (CEA/DM2S) en oct 2003 CC modification : A.Moal (CSSI) en jan 2004 CC CC Syntaxe generale : CC CC CHDIST MINTER CHFN CHDEP = IJET CHOLD CHNEW TOL TAB1 ; CC CC-------------------------------------------------------------------------- CE En entree : CE CE CHOLD : coordonnees des points extremites initiales des segments CE (type CHPOINT a 3 composantes) CE CHNEW : coordonnees des points extremites finales des segments CE (type CHPOINT a 3 composantes appuye sur le meme support CE que CH_OLD) CE TOL : tolerance (type FLOTTANT) CE TAB1 : Table qui doit contenir les parametres suivants en entree CE CE --> C Calcul de PnPn+1 : C DEX = XPN1(IZ,K) - XPN(IZ,K) DEY = YPN1(IZ,K) - YPN(IZ,K) DEZ = ZPN1(IZ,K) - ZPN(IZ,K) C C calcul de la distance de connection par defaut C CHPAS0(IZ,K) = (DEX*DEX + DEY*DEY + DEZ*DEZ)**.5 C C Test d'intersection avec le plan de la facette C IF ( (PS1*PS2 .LT. (-1*TOL1*TOL1)) . .OR. (PS2 .EQ. 0.) ) THEN C C -> --> -> C Calcul de AD1 = PnPn+1 . N : C AD1 = DEX*TXN(I,J) + DEY*TYN(I,J) + DEZ*TZN(I,J) C C Test de detection des lignes de champ paralleles a la facette C IF (ABS(AD1) .GT. TOL1) THEN C C C Calcul des points d'intersection M entre les lignes de champ et C le plan de la facette C XM = ( (XA1*DEX*TXN(I,J)) - . ((((YPN(IZ,K)-YA1)*DEX)- . (XPN(IZ,K)*DEY))*TYN(I,J)) - . ((((ZPN(IZ,K)-ZA1)*DEX)- . (XPN(IZ,K)*DEZ))*TZN(I,J)) ) / AD1 C YM = ( (YA1*DEY*TYN(I,J)) - . ((((ZPN(IZ,K)-ZA1)*DEY)- . (YPN(IZ,K)*DEZ))*TZN(I,J)) - . ((((XPN(IZ,K)-XA1)*DEY)- . (YPN(IZ,K)*DEX))*TXN(I,J)) ) / AD1 C ZM = ( (ZA1*DEZ*TZN(I,J)) - . ((((XPN(IZ,K)-XA1)*DEZ)- . (ZPN(IZ,K)*DEX))*TXN(I,J)) - . ((((YPN(IZ,K)-YA1)*DEZ)- . (ZPN(IZ,K)*DEY))*TYN(I,J)) ) / AD1 C C Calcul des vecteurs AM, BM, CM C XAM = XM - XA1 XBM = XM - XB1 XCM = XM - XC1 C C YAM = YM - YA1 YBM = YM - YB1 YCM = YM - YC1 C ZAM = ZM - ZA1 ZBM = ZM - ZB1 ZCM = ZM - ZC1 C C Denominateurs C DX = (YCM*ZBM)-(ZCM*YBM)+(ZAM*YBM)- . (ZAM*YCM)-(YAM*ZBM)+(YAM*ZCM) C DY = (ZCM*XBM)-(XCM*ZBM)+(XAM*ZBM)- . (XAM*ZCM)-(ZAM*XBM)+(ZAM*XCM) C DZ = (YCM*XBM)-(XCM*YBM)+(XAM*YBM)- . (XAM*YCM)-(YAM*XBM)+(YAM*XCM) C IF (ABS(DZ) .GT. TOL1) THEN D1 = DZ A1 = YCM*XBM-XCM*YBM B1 = YAM*XCM-XAM*YCM C1 = YBM*XAM-XBM*YAM ELSEIF (ABS(DY) .GT. TOL1) THEN D1 = DY A1 = ZCM*XBM-XCM*ZBM B1 = ZAM*XCM-XAM*ZCM C1 = ZBM*XAM-XBM*ZAM ELSE D1 = DX A1 = YCM*ZBM-ZCM*YBM B1 = YAM*ZCM-ZAM*YCM C1 = YBM*ZAM-ZBM*YAM ENDIF C C Coordonnees barycentriques C ALPHA = A1 / D1 BETA = B1 / D1 GAMMA = C1 / D1 C IF ( (ALPHA .GE. 0.) .AND. . (BETA .GE. 0.) .AND. . (GAMMA .GE. 0.) ) THEN C INTER(IZ,K)=.TRUE. C C C Calcul du Flux Normalise C CHFN(IZ,K) = TFN(I,J,1)*ALPHA + . TFN(I,J,2)*BETA + . TFN(I,J,3)*GAMMA C C Calcul des champs de deplacement : C CHDEX(IZ,K) = XPN(IZ,K)-XM CHDEY(IZ,K) = YPN(IZ,K)-YM CHDEZ(IZ,K) = ZPN(IZ,K)-ZM C C Calcul de d(Pn,M) : C CHPAS(IZ,K) = ((XPN(IZ,K)-XM)**2. + . (YPN(IZ,K)-YM)**2. + (ZPN(IZ,K)-ZM)**2.)**.5 ENDIF ENDIF ENDIF ENDIF 600 CONTINUE 550 CONTINUE 500 CONTINUE SEGDES MELVAL ENDIF 450 CONTINUE SEGDES MCHAML 400 CONTINUE SEGDES MCHELM C C-------------------------------------------------------- C C Creation du maillage MINTER des points d'intersection C NBNN = 1 NBELEM = NBP NBSOUS = 0 NBREF = 0 SEGINI MELEME MINTER = MELEME ITYPEL = 1 C C Construction de MINTER C IPT1 = ICHOLD SEGACT IPT1 NINTER = 0 DO 701 I=1,NZ1 DO 700 K=1,NBP IF (INTER(I,K)) THEN NINTER = NINTER+1 NUM(1,NINTER) = IPT1.NUM(1,K) ENDIF 700 CONTINUE 701 CONTINUE NBELEM = NINTER SEGADJ MELEME SEGDES MELEME,IPT1 C-------------------------------------------------------- C C Creation du CHPOINT CHDIST C C segment MCHPOI NAT = 1 NSOUPO = NZ1 SEGINI MCHPOI ICHDIS = MCHPOI IFOPOI = 2 JATTRI(1) = 2 C segment MSOUPO NC = 1 DO 901 I=1,NZ1 SEGINI MSOUPO IPCHP(I)=MSOUPO NOCOMP(1)= 'SCAL' C segment MELEME : maillage associe a CHDIST NBNN = 1 NBELEM = NBP NBSOUS = 0 NBREF = 0 C C LE CHAMPOINT CHDIST S'APPUIE SUR LE MEME MAILLAGE QUE CHOLD C IGEOC=ICHOLD C C Construction de MPOVAL C N = NBP SEGINI MPOVAL IPOVAL=MPOVAL C C On met la valeur de la distance de connection dans le chpoint CHDIST C Ce point n'existera plus dans le prochain maillage MINTER C DO 900 K=1,NBP IF (INTER(I,K)) THEN VPOCHA(K,1)=CHPAS(I,K) ELSE VPOCHA(K,1)=CHPAS0(I,K) ENDIF 900 CONTINUE C SEGDES MPOVAL SEGDES MSOUPO 901 CONTINUE SEGDES MCHPOI C C-------------------------------------------------------- C C Creation du CHPOINT CHFN C C segment MCHPOI NAT = 1 NSOUPO = NZ1 SEGINI MCHPOI ICHFN = MCHPOI IFOPOI = 2 JATTRI(1)=2 C segment MSOUPO NC = 1 DO 911 I=1,NZ1 SEGINI MSOUPO IPCHP(I)=MSOUPO NOCOMP(1)= 'SCAL' C segment MELEME : maillage associe a CHDIST NBNN = 1 NBELEM = NBP NBSOUS = 0 NBREF = 0 C C LE CHAMPOINT CHFN S'APPUIE SUR LE MEME MAILLAGE QUE CHOLD C IGEOC=ICHOLD C C Construction de MPOVAL C N = NBP SEGINI MPOVAL IPOVAL=MPOVAL C C On met la valeur de la distance de connection dans le chpoint CHDIST C Ce point n'existera plus dans le prochain maillage MINTER C DO 910 K=1,NBP IF (INTER(I,K)) THEN VPOCHA(K,1)=CHFN(I,K) ELSE VPOCHA(K,1)=0. ENDIF 910 CONTINUE C C C Desactivation des segments C SEGDES MPOVAL SEGDES MSOUPO 911 CONTINUE SEGDES MCHPOI C-------------------------------------------------------- C C Creation du CHPOINT CHDEP C C segment MCHPOI NAT = 1 NSOUPO = 1 SEGINI MCHPOI ICHDEP = MCHPOI IFOPOI = 2 JATTRI(1) = 1 C segment MSOUPO NC = 3 SEGINI MSOUPO IPCHP(1)=MSOUPO NOCOMP(1)= 'UX ' NOCOMP(2)= 'UY ' NOCOMP(3)= 'UZ ' C segment MELEME : maillage associe a CHDIST NBNN = 1 NBELEM = NBP NBSOUS = 0 NBREF = 0 C C LE CHAMPOINT CHDEP S'APPUIE SUR LE MEME MAILLAGE QUE MINTER C IGEOC=MINTER C C Construction de MPOVAL C N = NINTER SEGINI MPOVAL IPOVAL=MPOVAL C NINTER = 0 C C On met la valeur du champ de deplacement dans le chpoint CHDEP C DO 921 I=1,NZ1 DO 920 K=1,NBP IF (INTER(I,K)) THEN NINTER = NINTER+1 VPOCHA(NINTER,1)= 0. - CHDEX(I,K) VPOCHA(NINTER,2)= 0. - CHDEY(I,K) VPOCHA(NINTER,3)= 0. - CHDEZ(I,K) ENDIF 920 CONTINUE 921 CONTINUE C C C Desactivation des segments C SEGDES MPOVAL SEGDES MSOUPO SEGDES MCHPOI C-------------------------------------------------------- C C Ecriture des nouveaux objets en sortie C CALL ECROBJ('CHPOINT ',ICHDEP) CALL ECROBJ('CHPOINT ',ICHFN) C CALL ECROBJ('MAILLAGE',MINTER) C CALL ECROBJ('CHPOINT ',ICHDIS) C SEGSUP MTRAV1 SEGSUP MTRAV2 C RETURN END