C J3LOAD SOURCE CHAT 05/01/13 00:46:59 5004 SUBROUTINE J3LOAD(NUM,WWORK,NPTO,TOL,IRET) C---------------------------------------------------- C TRANSFERT DES FACE 3D EN FACE 2D C C PP 9/97 C Pierre Pegon/JRC Ispra C---------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION NUM(2,NPTO),P1(3) C -INC PPARAM -INC CCOPTIO -INC SMCOORD C SEGMENT WWORK REAL*8 PORIG(3),VNORM(3),VI(3),VJ(3) INTEGER FWORK INTEGER TWORK(NTROU) ENDSEGMENT C SEGMENT WORK REAL*8 XYC(2,NPTO) INTEGER IST(3,NPTO) REAL*8 DENS(NPTO) INTEGER JUN ENDSEGMENT C IRET=0 IF(NPTO.LE.2)THEN WRITE(IOIMP,*)'J3LOAD: nb de point d"une face insuffisant' IRET=IRET+1 RETURN ENDIF WORK=FWORK C C LE PREMIER POINT DU CONTOURS DEFINIT L'ORIGINE C NUME=NUM(1,1) IREF=(NUME-1)*(IDIM+1) DO IE1=1,3 PORIG(IE1)=XCOOR(IREF+IE1) ENDDO XYC(1,1)=0.D0 XYC(2,1)=0.D0 DENS(1)=XCOOR(IREF+4) C C AVEC LE SECOND POINT, ON DEFINIT LE PREMIER VECTEUR C NUME=NUM(1,2) IREF=(NUME-1)*(IDIM+1) DO IE1=1,3 P1(IE1)=XCOOR(IREF+IE1) VI(IE1)=P1(IE1)-PORIG(IE1) ENDDO XNORM=SQRT(VI(1)**2+VI(2)**2+VI(3)**2) IF(XNORM.LT.TOL)THEN WRITE(IOIMP,*)'J3LOAD: 2 point d"une face sont trop proches' IRET=IRET+1 RETURN ENDIF DO IE1=1,3 VI(IE1)=VI(IE1)/XNORM ENDDO XYC(1,2)=XNORM XYC(2,2)=0.D0 DENS(2)=XCOOR(IREF+4) C C AVEC LE POINT SUIVANT, OU UN AUTRE POINT, ON DEFINIT LE VECTEUR C NORMAL PUIS LE SECOND VECTEUR C WARNING: LA + GD COMPOSANTE DE VNOR EST TJ POSITIVE! C DO IE1=3,NPTO NUME=NUM(1,IE1) IREF=(NUME-1)*(IDIM+1) DO IE2=1,3 VJ(IE2)=XCOOR(IREF+IE2)-P1(IE2) P1(IE2)=XCOOR(IREF+IE2) ENDDO XNORM=SQRT(VJ(1)**2+VJ(2)**2+VJ(3)**2) IF(XNORM.LT.TOL)THEN WRITE(IOIMP,*)'J3LOAD: 2 point d"une face sont trop proches' IRET=IRET+1 RETURN ENDIF DO IE2=1,3 VJ(IE2)=VJ(IE2)/XNORM ENDDO VNORM(1)=VI(2)*VJ(3)-VI(3)*VJ(2) VNORM(2)=VI(3)*VJ(1)-VI(1)*VJ(3) VNORM(3)=VI(1)*VJ(2)-VI(2)*VJ(1) XNORM=SQRT(VNORM(1)**2+VNORM(2)**2+VNORM(3)**2) IF(XNORM.GT.100*TOL)THEN XCMAX=VNORM(1) DO IE2=2,3 IF(ABS(VNORM(IE2)).GT.ABS(XCMAX))THEN XCMAX=VNORM(IE2) ENDIF ENDDO XNORM=SIGN(1.D0,XCMAX)*XNORM DO IE2=1,3 VNORM(IE2)=VNORM(IE2)/XNORM ENDDO VJ(1)=VNORM(2)*VI(3)-VNORM(3)*VI(2) VJ(2)=VNORM(3)*VI(1)-VNORM(1)*VI(3) VJ(3)=VNORM(1)*VI(2)-VNORM(2)*VI(1) GOTO 1 ENDIF ENDDO WRITE(IOIMP,*)'J3LOAD: On ne reussit pas a definir le 2nd vecteur' IRET=IRET+1 RETURN C C ON BOUCLE SUR LES POINTS ET ON LES PROJETTE SUR LE PLAN C ON TESTE LA DISTANCE DU POINT AU PLAN C 1 CONTINUE DO IE1=3,NPTO NUME=NUM(1,IE1) IREF=(NUME-1)*(IDIM+1) DO IE2=1,3 P1(IE2)=XCOOR(IREF+IE2) ENDDO XDIST=(P1(1)-PORIG(1))*VNORM(1)+(P1(2)-PORIG(2))*VNORM(2) > +(P1(3)-PORIG(3))*VNORM(3) IF(ABS(XDIST).GT.TOL)THEN WRITE(IOIMP,*)'J3LOAD: 1 point d"une face est hors plan' IRET=IRET+1 RETURN ENDIF XYC(1,IE1)=(P1(1)-PORIG(1))*VI(1)+(P1(2)-PORIG(2))*VI(2) > +(P1(3)-PORIG(3))*VI(3) XYC(2,IE1)=(P1(1)-PORIG(1))*VJ(1)+(P1(2)-PORIG(2))*VJ(2) > +(P1(3)-PORIG(3))*VJ(3) DENS(IE1)=XCOOR(IREF+4) ENDDO C JUN=0 C RETURN END